{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Clouds - A Topological Approach\n\nRecently I've had a strong interest in topological approaches to data analysis (typically abbreviated TDA). To those unfamiliar, the key idea in topology is to throw away coordinate systems and analyze systems and structures based soley on _connectedness_. For phenomena with complicated, intricate structure (i.e. clouds), TDA is good at extracting global qualities. For things like facial recognition, one is better off using conventional CNN techniques, which rely heavily on layers of feature detection.\n\nSince this competition is about recognizing a complicated pattern, it seems appropriate to use TDA. An example process might use TDA to extract local features, then use a CNN to wrap these up into a segmentation map. In this kernel, I mostly investigate topological features qualitatively to see if it is a viable option. A rough segmentation is attempted at the end to validate it, but it is so slow as to be practically unusable.\n\nFor more information, see I. Obayashi and Y. Hiraoka, “Persistence Diagrams with Linear Machine Learning Models.”\n\n## Setup"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"%matplotlib inline\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport matplotlib\nmatplotlib.rcParams['figure.figsize'] = [15,15]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The first step is to load the training data. I do a couple of transformations to clean it up and make it easier to work with."},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_train():\n    df = pd.read_csv('../input/understanding_cloud_organization/train.csv')\n    df = df.dropna()\n    df[['fname','label']] = df.Image_Label.str.split('_',expand=True)\n    df.EncodedPixels = df.EncodedPixels.str.split().apply(np.array,dtype=int)\n    df = df.set_index(['fname','label']).EncodedPixels\n    df = df.loc[~df.index.duplicated(keep='first')]\n    df = df.unstack('label')\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = read_train()\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The next step is to define an easy function to load images. I use PIL (like everybody else) for this."},{"metadata":{"trusted":true},"cell_type":"code","source":"def load_image(fname,test=False):\n    from PIL import Image\n    path = '../input/understanding_cloud_organization/{}_images/{}'.format('train' if not test else 'test', fname)\n    return np.asarray(Image.open(path))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fname = df.index[0]\nimg = load_image(fname)\nplt.imshow(img);","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Wow, this has extremely intricate structure! To transform it into something we can apply persistent homology to, we need to apply the following steps:\n\n1. Shrink the image (optional, but speeds it up a lot).\n2. Threshold the image. This effective separates the clouds from the ocean.\n3. Apply the exact euclidean distance transform to the thresholded image and it's inverse. Subtract the two."},{"metadata":{"trusted":true},"cell_type":"code","source":"def transform_image(x,thresh,rescale,blur):\n    import numpy as np\n    from scipy.ndimage import zoom,gaussian_filter,distance_transform_edt\n    # normalize between 0 and 1\n    img = np.asarray(x,dtype=float)/255\n\n    # take average of RGB channels\n    img = img.mean(axis=2)\n\n    # Shrink image\n    img = zoom(img,rescale)\n\n    # Remove artifacts from shrinking\n    img = gaussian_filter(img,blur)\n\n    # threshold image\n    img = (img < thresh).astype(int)\n\n    # Compute distance from black pixels to nearest white pixels (and normalize distances)\n    img = distance_transform_edt(img) - (distance_transform_edt(1-img)-1).clip(0,None)\n    img = img/rescale\n\n    return img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"timg = transform_image(img,thresh=0.5,rescale=0.2,blur=0.1)\nplt.imshow(timg,cmap='coolwarm')\nplt.clim([-50,50])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The goal of the transformation is to highlight the intricate structure of the cloud formations by creating a filtration. Specifically, I do a lower-star filtration (i.e. the birth points of connected components occur within clouds). Furthermore, since I'm using the Ripser library, I only compute the 0-dimensional homology groups (connected components). Looking at the picture above, I can visibly pick out a cycle, so there is the possibility that the 1-dimensional groups may also be of value.<sup>1</sup> To visulize what the persistence code is doing, let's look at what happens when we threshold the transformed image at different values (i.e. distances).\n\n<sup>1</sup>As a side note, it looks like it would be pretty minor to make that change to the Ripser python code."},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.imshow(np.block([[timg < t for t in ts] for ts in np.linspace(-20,20,16).reshape(4,4)]), cmap='gray');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Ok, so this makes sense. As the threshold increases, the sublevel sets appear, merge, and then fill up the image. The final step is to create a persistence image."},{"metadata":{"trusted":true},"cell_type":"code","source":"from ripser import lower_star_img as lower_star\n\ndef pers_image(dgm,res,rng,spread):\n    from scipy.ndimage import gaussian_filter,zoom\n    \n    # birth and death\n    idx = np.where(~np.isinf(dgm).any(axis=1))\n    b,d = dgm[idx].T\n    \n    # compute histogram at 3x resolution\n    img = np.histogram2d(b,d,bins=res*3,range=[[-rng,rng],[-rng,rng]])[0]\n    \n    # apply blurring to histogram (relative to total shape)\n    img = gaussian_filter(img, np.array(img.shape)*spread)\n    \n    # decimate histogram (hopefully the blurring prevents aliasing)\n    img = img[::3,::3]\n    \n    return img","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"dgm = lower_star(timg)\npim = pers_image(dgm,res=500,rng=20,spread=0.02)\nplt.imshow(pim.T,origin='lower',cmap='jet');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Huh, looks pretty concentrated in the middle. Roughly, these correspond to small connected components (maybe a pixel or two wide). They seem to be fairly uniformly separated and die at about the same time.\n\nFinally, we need some utility functions to mask the images and convert the persistence images to vectors."},{"metadata":{"trusted":true},"cell_type":"code","source":"def apply_mask(img,mask):\n    px1 = mask[0::2]\n    px2 = px1 + mask[1::2]\n    mask = np.zeros(img.shape[:2][::-1],dtype=int)\n    mask.flat[px1] = 1\n    mask.flat[px2[px2<len(mask.flat)]] = -1\n    mask.flat = np.cumsum(mask.flat)\n    mask = mask.T\n    return img*(mask>0)[:,:,None]\n\ndef training_vectors(df,thresh,rescale,blur,res,spread,rng):\n    from tqdm import tqdm_notebook\n    todo = [(label,col,fname,mask) for label,col in df.items() for fname,mask in col.dropna().items()]\n    vecs,labels = [],[]\n    for label,col,fname,mask in tqdm_notebook(todo):\n        img = load_image(fname)\n        img = apply_mask(img,mask)\n        img = transform_image(img,thresh,rescale,blur)\n        dgm = lower_star(img)\n        img = pers_image(dgm,res,rng,spread)\n        vecs += [img]\n        labels += [label]\n    return np.array(vecs),np.array(labels)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fname2 = df.index[3]\nimg2 = apply_mask(load_image(fname2),df.Flower.iloc[3])\nplt.imshow(img2); plt.show();","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Training\n\nNow we can do some regular old-fashioned machine learning. The goal right now is just to classify the masks. It is especially useful if the algorithm has a coefficient or weight vector since we can visualize it directly as a persistence image (logistic regression, SVM, etc). Since this is a kernel I chose to scale things down as much as possible, probably at the cost of accuracy."},{"metadata":{"trusted":true},"cell_type":"code","source":"X,y = training_vectors(df.sample(1000),thresh=0.5,rescale=0.2,blur=0.1,res=20,rng=20,spread=0.02)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.linear_model import LogisticRegression\nfrom scipy.ndimage import gaussian_filter\nlg = LogisticRegression(C=1e-1,solver='lbfgs',multi_class='multinomial',max_iter=300)\nXi = gaussian_filter(X,(0,1,1)).reshape(-1,20**2)\nidx = np.arange(len(Xi))\nnp.random.shuffle(idx)\ni1,i2 = idx[:-300],idx[-300:]\nlg.fit(Xi[i1],y[i1])\nlg.score(Xi[i2],y[i2])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's look at some of the coefficient vectors."},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy.ndimage import zoom\nfig,axes = plt.subplots(2,2,figsize=[15,15])\nfor ci,ax in zip(lg.coef_,axes.flat):\n    ax.imshow(zoom(ci.reshape(20,20).T,30),origin='lower',cmap='RdBu');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.svm import LinearSVC\n\nsvm = LinearSVC(max_iter=500,C=1e-1)\nsvm.fit(Xi[i1],y[i1])\nsvm.score(Xi[i2],y[i2])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fig,axes = plt.subplots(2,2,figsize=[15,15])\nfor ci,ax in zip(svm.coef_,axes.flat):\n    ax.imshow(zoom(ci.reshape(20,20).T,30),origin='lower',cmap='RdBu');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Segmentation\n\nFor a mostly topological method, this performs pretty well! My theory is that the approach can be used to generate topological features that can be used by a CNN. For now, let's see if we can do segmentation with this method.[](http://)"},{"metadata":{"trusted":true},"cell_type":"code","source":"def segment(img,size,stride,thresh,rescale,blur,res,rng,spread):\n    from tqdm import tqdm_notebook\n    from itertools import product\n    trans = transform_image(img,thresh,rescale,blur)\n    preds = np.empty(np.array(img.shape[:2])//stride+1,dtype='U10')\n    for i,j in tqdm_notebook(list(product(range(0,img.shape[0],stride),range(0,img.shape[1],stride)))):\n        i1,i2,j1,j2 = ((np.array([[i],[j]])+[-size/2,size/2])*rescale).astype(int).flat\n        small_img = trans[i1:i2,j1:j2]\n        dgm = lower_star(small_img)\n        pim = pers_image(dgm,res,rng,spread)\n        preds[i//stride,j//stride] = svm.predict(pim.reshape(1,-1))[0]\n    return preds","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"si = 121\nsimg = load_image(df.index[si])\nseg = segment(simg,size=500,stride=50,thresh=0.5,rescale=0.2,blur=0.1,res=20,rng=20,spread=0.02)\n\nfig,axes = plt.subplots(4,2,figsize=[15,20])\ndice = 0\nfor ax,label in zip(axes,'Flower Fish Gravel Sugar'.split()):\n    segl = seg == label\n    try:\n        ax[0].imshow(apply_mask(img,df[label].iloc[si]));\n        segY = apply_mask(np.ones_like(img,dtype=int),df[label].iloc[si])[:,:,0] == 1\n    except TypeError:\n        segY = np.zeros(img.shape[:2],dtype=bool)\n        ax[0].imshow(np.zeros_like(img))\n    segX = (zoom(segl.astype(float),50)[:img.shape[0],:img.shape[1]] > 0.5) == 1\n    if (segX|segY).sum() == 0:\n        dice += 0.25\n    else:\n        dice += (2*(segX&segY).sum()/(segX|segY).sum())/4\n    ax[1].imshow(segX,cmap='Blues');\n    ax[0].set_title(label)\nprint(dice)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Conclusion\n\nThe segmentation is so-so, but this is likely due to all the rough approximations I made to keep the computation time manageable for playing around. My goal is to refine the approach in another notebook, so stay tuned. At the very least, this has been a learning experience for me in TDA approaches. Feel free to leave comments, especially if you notice any problems or errors in my code."}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}