{"cells":[{"metadata":{},"cell_type":"markdown","source":"Since I am not at all competitive in this competition, I share my code in case it is useful to someone else.\n\nhttps://github.com/junkoda/junkoda_cellularlib\n\nThis is not a model and does not contribute directory to the score, but maybe useful indirectly.\n\nThe code finds nuclei in channel 1.\n\n* My original idea was to use image clips centered on nuclei, rotating them along the major axis; not sure how useful this is.\n* The pixel values of nuclei may be useful as a characteristic brightness of the image, which depends on the observational condition, to calibrate the overall pixel values.\n\n\n"},{"metadata":{},"cell_type":"markdown","source":"## Install\n\n### Requirements:\n\n* C++11 compiler\n* Eigen linear algebra library "},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"# Install Eigen\n# You need to turn on \"Internet\" from Settings on the right panel\n! apt install libeigen3-dev\n\n# Download code from github\n%cd /tmp\n! git clone https://github.com/junkoda/junkoda_cellularlib.git\n\n%cd /tmp/junkoda_cellularlib/py\n\n# Specify the version in case backward compatibility is lost\n! git checkout bc97d6ab419d8e9e1c295a7662d94cfd1f5b3501\n\n# Compile the library\n! make","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Load modules"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"%matplotlib inline\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport pandas as pd\nimport junkoda_cellularlib as cellularlib\n\nfrom scipy.signal import find_peaks\nfrom sklearn.neighbors import KernelDensity","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"cellularlib.data.set_data_dir('/kaggle/input/recursion-cellular-image-classification')\nimg = cellularlib.data.load('train', 'RPE-01_1_B03', site=1)\n\nprint('img shape', img.shape)\n\nplt.title('Channel 1')\nplt.imshow(img[0, ])  # channel 1\nplt.clim(0, 0.1)\nplt.show()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Find the brightness distribution of nuclei pixels\n\nThis may be useful to determine the overall normalisation of the image brightness.\n\n### Definition\n\n- *cluster* is a connected component of pixels above some threshold; img >= threshold;\n- 4 adjacent pixels are connected;\n- *size* of a cluster is the number of pixels that belongs to that cluster;\n\n\n### Algorithm\n\nPixels are explored in a \"watershed\" way, from highest peak to lower. A pixel belongs to nuclei if it belongs to some cluster of size between `size_min` and `size_max` for some threshold. The code is in `watershed_nuclei.cpp`."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Binary mask: 1 for nuclei pixel in the algorithm above, 0 otherwise\nnuclei = cellularlib.threshold.obtain_nuclei_pixels(img[0,], size_min=50, size_max=1000)\n# Good size_min may depend on the cell_line, HUVEC2, RPE, ...\n\nprint('shape', nuclei.shape)\n\nplt.figure(figsize=(8, 4))\nplt.subplot(1, 2, 1)\nplt.title('Channel 1')\nplt.imshow(img[0, ])\nplt.clim(0, 0.1)\n\nplt.subplot(1, 2, 2)\nplt.title('Nuclei pixels')\nplt.imshow(nuclei, cmap='Greys')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Pixel values of nuclei pixels (converted to range [0, 255])\nvalues = img[0,].flatten()[nuclei.flatten()]*255\n\n# Kernel density estimation\nkde = KernelDensity(kernel='gaussian', bandwidth=1.0)\n# Extra smoothness, e.g., bandwidth = 2 may be better for stability.\nkde.fit(values.reshape(-1, 1))\nx = np.arange(40)\nf = np.exp(kde.score_samples(x.reshape(-1, 1)))*len(values)\n\n\n# Find peaks\nheight_threshold = 0.1*np.max(f)\npeaks = find_peaks(f, height_threshold)[0]\nx_peaks = x[peaks]\nnuclei_brightness = x_peaks[-1]\n\nbins = 0.5 + np.arange(40)\nplt.figure()\nplt.title('Nuclei brightness')\nplt.xlim(0, 25)\nplt.ylim(0, 20000)\nplt.xlabel('pixel value [0, 255]')\nplt.hist(img[0, ].flatten()*255, bins, alpha=0.5, label='All pixels')\nplt.hist(values, bins, alpha=0.5, label='Nuclei pixels')\nplt.plot(x, f, label='KDE')\nplt.axvline(nuclei_brightness, ls='--', color='gray', label='nuclei brightness')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Fitting nuclei with ellipses\n\nThe mean and covariance matrix of each cluster, or in physics, center of mass and moment of inertia, are computed. Diagonalizing the matrix gives the major and minor axes of an ellipse."},{"metadata":{"trusted":true},"cell_type":"code","source":"threshold = 0.6*nuclei_brightness/255  # convert back [0, 255] to [0, 1.0]\nes = cellularlib.ellipses.obtain(img[0,], threshold, 50)  # 50 is the minimum size\n\n# size x y semi-major axis semi-minor axis theta\n# theta is the angle between the major axis and x axis in radian\nes[:3, ]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(6, 6))\ncellularlib.ellipses.plot(es, img[0,])\nplt.clim(0, 0.05)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Number of clusters - threshold relation"},{"metadata":{"trusted":true},"cell_type":"code","source":"thresholds, nclusters = cellularlib.compute_nclusters(img[0,], size_threshold=5)\n\nplt.xlabel('threshold')\nplt.ylabel('number of clusters')\nplt.xlim(0, 0.2)\nplt.plot(thresholds, nclusters)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"From right to left, as the threshold goes down, the number of clusters first increses, then as clusters merge, the number decreases, and then there is a second peak on the left due to noise."},{"metadata":{},"cell_type":"markdown","source":""}],"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}