{"cells":[{"cell_type":"markdown","metadata":{"_cell_guid":"419ad209-85bf-c4ef-d3d7-d8fa2b6be0ea"},"source":"# Cervix EDA\n\nIn this competition we have a multi-class classification problem with **three** classes. We are asked, given an image, to identify the cervix type.\n\nFrom the data description:\n\n*In this competition, you will develop algorithms to correctly classify cervix types based on cervical images. These different types of cervix in our data set are all considered normal (not cancerous), but since the transformation zones aren't always visible, some of the patients require further testing while some don't. This decision is very important for the healthcare provider and critical for the patient. Identifying the transformation zones is not an easy task for the healthcare providers, therefore, an algorithm-aided decision will significantly improve the quality and efficiency of cervical cancer screening for these patients.*\n\nThe submission format is asking for a probability for each of the three different cervix types.\n\nI'll show some of the basic properties of this dataset and a few sample images from the training set."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"6f6e8717-2620-3fbc-dff4-7093bde04be9"},"outputs":[],"source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom skimage.io import imread, imshow\nimport cv2\n\n%matplotlib inline\nimport plotly.offline as py\npy.init_notebook_mode(connected=True)\nimport plotly.graph_objs as go\nimport plotly.tools as tls\n\nfrom subprocess import check_output\nprint(check_output([\"ls\", \"../input/train\"]).decode(\"utf8\"))"},{"cell_type":"markdown","metadata":{"_cell_guid":"a640fa62-669c-a08c-aa2b-20c4c713bcdd"},"source":"We are given training images for each of cervix types. Lets first count them for each class."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"df96bdda-f3bf-84f9-0abd-ee6161018846"},"outputs":[],"source":"from glob import glob\nbasepath = '../input/train/'\n\nall_cervix_images = []\n\nfor path in sorted(glob(basepath + \"*\")):\n    cervix_type = path.split(\"/\")[-1]\n    cervix_images = sorted(glob(basepath + cervix_type + \"/*\"))\n    all_cervix_images = all_cervix_images + cervix_images\n\nall_cervix_images = pd.DataFrame({'imagepath': all_cervix_images})\nall_cervix_images['filetype'] = all_cervix_images.apply(lambda row: row.imagepath.split(\".\")[-1], axis=1)\nall_cervix_images['type'] = all_cervix_images.apply(lambda row: row.imagepath.split(\"/\")[-2], axis=1)\nall_cervix_images.head()"},{"cell_type":"markdown","metadata":{"_cell_guid":"172accdb-dc4d-17da-fa39-0740c32a54f9"},"source":"## Image types\n\nNow that we have the data in a handy dataframe we can do a few aggregations on the data. Let us first see how many images there are for each cervix type and which file types they have.\n\nAll files are in JPG format and Type 2 is the most common one with a little bit more than 50% in the training data in total, Type 1 on the other hand has a little bit less than 20% in the training data."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"1c42f54f-45fb-7a2e-6c89-ba7c7a6a48a3"},"outputs":[],"source":"print('We have a total of {} images in the whole dataset'.format(all_cervix_images.shape[0]))\ntype_aggregation = all_cervix_images.groupby(['type', 'filetype']).agg('count')\ntype_aggregation_p = type_aggregation.apply(lambda row: 1.0*row['imagepath']/all_cervix_images.shape[0], axis=1)\n\nfig, axes = plt.subplots(nrows=2, ncols=1, figsize=(10, 8))\n\ntype_aggregation.plot.barh(ax=axes[0])\naxes[0].set_xlabel(\"image count\")\ntype_aggregation_p.plot.barh(ax=axes[1])\naxes[1].set_xlabel(\"training size fraction\")"},{"cell_type":"markdown","metadata":{"_cell_guid":"4b3581d5-c9a3-49e4-427c-08ef8ede1053"},"source":"Now, lets read the files for each type to get an idea about how the images look like.\n\nThe images seem to vary alot in they formats, the first two samples have only a circular area with the actual image, the last sample has the image in a rectangle."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"0b500f12-8d54-4dab-a9d1-8bc2fb2de755"},"outputs":[],"source":"fig = plt.figure(figsize=(12,8))\n\ni = 1\nfor t in all_cervix_images['type'].unique():\n    ax = fig.add_subplot(1,3,i)\n    i+=1\n    f = all_cervix_images[all_cervix_images['type'] == t]['imagepath'].values[0]\n    plt.imshow(plt.imread(f))\n    plt.title('sample for cervix {}'.format(t))"},{"cell_type":"markdown","metadata":{"_cell_guid":"4db0231d-6d7c-5780-400d-d27cbea25434"},"source":"## Image dimensions\n\nNow, in order to get an idea of how many different shapes of images by class there are, lets have a look at. To reduce runtime, take only a subsample per class."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"24dca5d4-1055-00b3-d5e1-a9eb8de6cb64"},"outputs":[],"source":"from collections import defaultdict\n\nimages = defaultdict(list)\n\nfor t in all_cervix_images['type'].unique():\n    sample_counter = 0\n    for _, row in all_cervix_images[all_cervix_images['type'] == t].iterrows():\n        #print('reading image {}'.format(row.imagepath))\n        try:\n            img = imread(row.imagepath)\n            sample_counter +=1\n            images[t].append(img)\n        except:\n            print('image read failed for {}'.format(row.imagepath))\n        if sample_counter > 50:\n            break"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"c3aa56f1-e78e-55e4-a961-52a189c18403"},"outputs":[],"source":"dfs = []\nfor t in all_cervix_images['type'].unique():\n    t_ = pd.DataFrame(\n        {\n            'nrows': list(map(lambda i: i.shape[0], images[t])),\n            'ncols': list(map(lambda i: i.shape[1], images[t])),\n            'nchans': list(map(lambda i: i.shape[2], images[t])),\n            'type': t\n        }\n    )\n    dfs.append(t_)\n\nshapes_df = pd.concat(dfs, axis=0)\nshapes_df_grouped = shapes_df.groupby(by=['nchans', 'ncols', 'nrows', 'type']).size().reset_index().sort_values(['type', 0], ascending=False)\nshapes_df_grouped"},{"cell_type":"markdown","metadata":{"_cell_guid":"dabf9549-9b4c-11b8-6182-ec71de2aab5d"},"source":"All of the images in our sample have three channels, we can ignore this information for now. Now lets build a barplot to get an idea of the distribution of image dimensions by cervix type."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"336a5362-4eee-090d-5ac5-3a84e6970b00"},"outputs":[],"source":"shapes_df_grouped['size_with_type'] = shapes_df_grouped.apply(lambda row: '{}-{}-{}'.format(row.ncols, row.nrows, row.type), axis=1)\nshapes_df_grouped = shapes_df_grouped.set_index(shapes_df_grouped['size_with_type'].values)\nshapes_df_grouped['count'] = shapes_df_grouped[[0]]\n\nplt.figure(figsize=(10,8))\n#shapes_df_grouped['count'].plot.barh(figsize=(10,8))\nsns.barplot(x=\"count\", y=\"size_with_type\", data=shapes_df_grouped)"},{"cell_type":"markdown","metadata":{"_cell_guid":"1008461b-6a10-64a2-1716-4ba7d2d2dce4"},"source":"# TSNE embedding\n\nWe will now take all of the sample images, rescale them & convert them to grayscale. This will result in a matrix, where each row are all flattened pixel for the grayscale images.\n\nThe original images have, as we have seen earlier, quite a high resolution, so scaling them down to **100 x 100** is resulting in a great loss of information, so the embedding to two dimensions is likely not going to have a good structure where we can separate visually by cervical cancer types. Also, we are giving only very few images per class that TSNE can work with to find a good, distance preserving, embedding.\n\nI've added the option to choose to convert them to grayscale before passing them to TSNE. Per default each image will now keep its RGB information and therefore will be transformed to a 100*100*3 = 30.000 dimensional vector."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"f83a4a94-b331-a1a5-def9-68678758c30b"},"outputs":[],"source":"def transform_image(img, rescaled_dim, to_gray=False):\n    resized = cv2.resize(img, (rescaled_dim, rescaled_dim), cv2.INTER_LINEAR)\n\n    if to_gray:\n        resized = cv2.cvtColor(resized, cv2.COLOR_RGB2GRAY).astype('float')\n    else:\n        resized = resized.astype('float')\n\n    normalized = cv2.normalize(resized, None, 0.0, 1.0, cv2.NORM_MINMAX)\n    timg = normalized.reshape(1, np.prod(normalized.shape))\n\n    return timg/np.linalg.norm(timg)\n\n#rescaled_dim = 100\nrescaled_dim = 2322\n\nall_images = []\nall_image_types = []\n\nfor t in all_cervix_images['type'].unique():\n    all_images = all_images + images[t]\n    all_image_types = all_image_types + len(images[t])*[t]\n\nprint('Resized to 2322x2322')\n# - normalize each uint8 image to the value interval [0, 1] as float image\n# - rgb to gray\n# - downsample image to rescaled_dim X rescaled_dim\n# - L2 norm of each sample = 1\ngray_all_images_as_vecs = [transform_image(img, rescaled_dim) for img in all_images]\n\ngray_imgs_mat = np.array(gray_all_images_as_vecs).squeeze()\nall_image_types = np.array(all_image_types)\ngray_imgs_mat.shape, all_image_types.shape"},{"cell_type":"markdown","metadata":{"_cell_guid":"ec6a4789-f799-8ce7-7854-0a1fe0a88baf"},"source":"### 3D t-SNE with cervix indicators\n\nNow let's project the 100x100x3 images to three dimensions to check for low dimensional patterns."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"e7fd77ee-8449-2ca7-2691-1906ee10286b"},"outputs":[],"source":"from sklearn.manifold import TSNE\ntsne = TSNE(\n    n_components=3,\n    init='random', # pca\n    random_state=101,\n    method='barnes_hut',\n    n_iter=500,\n    verbose=2\n).fit_transform(gray_imgs_mat)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"508824d4-7475-5722-bb85-73c1bd79c1ee"},"outputs":[],"source":"from sklearn import preprocessing\n\ntrace1 = go.Scatter3d(\n    x=tsne[:,0],\n    y=tsne[:,1],\n    z=tsne[:,2],\n    mode='markers',\n    marker=dict(\n        sizemode='diameter',\n        color = preprocessing.LabelEncoder().fit_transform(all_image_types),\n        colorscale = 'Portland',\n        colorbar = dict(title = 'cervix types'),\n        line=dict(color='rgb(255, 255, 255)'),\n        opacity=0.9\n    )\n)\n\ndata=[trace1]\nlayout=dict(height=800, width=800, title='3D embedding of images')\nfig=dict(data=data, layout=layout)\npy.iplot(fig, filename='3DBubble')"},{"cell_type":"markdown","metadata":{"_cell_guid":"5f8e4fcc-7686-ed91-bcfb-78a7384cf911"},"source":"We can clearly see that there is a big heterogeneous cluster and a few quite distant outliers."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"c9a14afc-2aab-8519-b4d2-53c9eb36337c"},"outputs":[],"source":"for t in all_cervix_images['type'].unique():\n    tsne_t = tsne[np.where(all_image_types == t), :][0]\n    plt.scatter(tsne_t[:, 0], tsne_t[:, 1])\nplt.legend(all_cervix_images['type'].unique())"},{"cell_type":"markdown","metadata":{"_cell_guid":"dba5ac7f-3868-04c6-79f6-3236bfea95a1"},"source":"# Image clustering with TSNE embedding\n\nLet us now use the compressed images and the distance-preserving transformation of TSNE to get an understanding of what image clusters there are.\n\nFor that we will use matplotlib to replace the previously red, green and blue points with the actual images that were input to the TSNE transformation."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"72cd02f8-1bb5-77d6-3029-6b78111abe6c"},"outputs":[],"source":"from matplotlib.offsetbox import OffsetImage, AnnotationBbox\ndef imscatter(x, y, images, ax=None, zoom=0.01):\n    ax = plt.gca()\n    images = [OffsetImage(image, zoom=zoom) for image in images]\n    artists = []\n    for x0, y0, im0 in zip(x, y, images):\n        ab = AnnotationBbox(im0, (x0, y0), xycoords='data', frameon=False)\n        artists.append(ax.add_artist(ab))\n    ax.update_datalim(np.column_stack([x, y]))\n    ax.autoscale()\n    #return artists\n\nnimgs = 60\nplt.figure(figsize=(10,8))\nimscatter(tsne[0:nimgs,0], tsne[0:nimgs,1], all_images[0:nimgs])"},{"cell_type":"markdown","metadata":{"_cell_guid":"660e028d-fba5-e767-850e-6eaa2c463f84"},"source":"# Clustering of pairwise image distances\n\nTo get a different view of how images relate to each other from a purely numerical point of view, lets now look at pairwise distances. For that we'll use scipy's pdist method.\n\nThe yellow somewhat clustered area tells us there are a few images that have relatively high distance to all other images in the sample of our training images we read.\nOn the left and top of the clustermap we find one of three colors for each row and column, this color indicates the type of cervix.\n\n* Type 1: Red\n* Type 2: Green\n* Type 3: Blue"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"fe9f7daa-880d-b944-9705-7bec5a550e32"},"outputs":[],"source":"pal = sns.color_palette(\"hls\", 3)\nsns.palplot(pal)"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"f54fc443-44d5-69a9-0d2a-1d506ebd467f"},"outputs":[],"source":"from scipy.spatial.distance import pdist, squareform\n\nsq_dists = squareform(pdist(gray_imgs_mat))\n\nall_image_types = list(all_image_types)\n\nd = {\n    'Type_1': pal[0],\n    'Type_2': pal[1],\n    'Type_3': pal[2]\n}\n\n# translate each sample to its color\ncolors = list(map(lambda t: d[t], all_image_types))\n\nsns.clustermap(\n    sq_dists,\n    figsize=(12,12),\n    row_colors=colors, col_colors=colors,\n    cmap=plt.get_cmap('viridis')\n)"},{"cell_type":"markdown","metadata":{"_cell_guid":"3c09b509-3249-d7e0-2c07-f9df3c81cb72"},"source":"Here is the unclustered distance matrix."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"1becfcaa-0503-2ee6-dc38-8d29ac8e463a"},"outputs":[],"source":"\nmask = np.zeros_like(sq_dists, dtype=np.bool)\nmask[np.triu_indices_from(mask)] = True\n\nplt.figure(figsize=(12,12))\nsns.heatmap(sq_dists, cmap=plt.get_cmap('viridis'), square=True, mask=mask)"},{"cell_type":"markdown","metadata":{"_cell_guid":"763f409d-3e8d-2045-8014-cbd2d4e48001"},"source":"----------\n\n# Image neighbourhood\n\nGiven the pool of cervix images, and the similarity matrix above, we will now plot the extrem cases of the similarity matrix.\n\nTo do that we find the image that has:\n\n* maximal distance on average compared to the all other images\n* minimal distance on average compared to the all other images\n\nand compare both candidates to the mean image from out set of images."},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"32ac6a8f-4636-1f6a-b762-1d4992704a1d"},"outputs":[],"source":"# upper triangle of matrix set to np.nan\nsq_dists[np.triu_indices_from(mask)] = np.nan\nsq_dists[0, 0] = np.nan\n\nfig = plt.figure(figsize=(12,8))\n# maximally dissimilar image\nax = fig.add_subplot(1,3,1)\nmaximally_dissimilar_image_idx = np.nanargmax(np.nanmean(sq_dists, axis=1))\nplt.imshow(all_images[maximally_dissimilar_image_idx])\nplt.title('maximally dissimilar')\n\n# maximally similar image\nax = fig.add_subplot(1,3,2)\nmaximally_similar_image_idx = np.nanargmin(np.nanmean(sq_dists, axis=1))\nplt.imshow(all_images[maximally_similar_image_idx])\nplt.title('maximally similar')\n\n# now compute the mean image\nax = fig.add_subplot(1,3,3)\nmean_img = gray_imgs_mat.mean(axis=0).reshape(rescaled_dim, rescaled_dim, 3)\nplt.imshow(cv2.normalize(mean_img, None, 0.0, 1.0, cv2.NORM_MINMAX))\nplt.title('mean image')"},{"cell_type":"markdown","metadata":{"_cell_guid":"940d21df-393e-99e2-8a17-0b314167fc75"},"source":"The image on the left has quite a lot of blue in it which I would not expect to be present in the majority of samples. I think it is a good representation of maximally dissimilar."},{"cell_type":"markdown","metadata":{"_cell_guid":"af9d24c3-b3e0-6fc5-eca6-0eb5f672c793"},"source":"----------"},{"cell_type":"code","execution_count":null,"metadata":{"_cell_guid":"a7b22b2d-86d2-1567-ba25-57d734402371"},"outputs":[],"source":""}],"metadata":{"_change_revision":0,"_is_fork":false,"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.6.0"}},"nbformat":4,"nbformat_minor":0}