{"cells":[{"metadata":{"trusted":false,"_uuid":"e2553cdaab6b228fbee45fd7ff86735b37abd94d"},"cell_type":"code","source":"# load python libraries\n%matplotlib inline\n\nimport matplotlib.image as img\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport seaborn as sns\nimport csv\nimport sys\nimport math\nfrom sklearn.decomposition import PCA\nfrom os import listdir","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"8bf947edef11d395011988c2ff66089828eefa72"},"cell_type":"markdown","source":"The following image will be taken as an example"},{"metadata":{"scrolled":true,"trusted":false,"_uuid":"f86f6005d9ee38e1866ca04f4c89bfa22394a76d"},"cell_type":"code","source":"picture = img.imread('../input/stage1_test/df40099c6306ca1f47fcc8a62e2fa39486d4e223177afdc51b2ad189691802d8/images/df40099c6306ca1f47fcc8a62e2fa39486d4e223177afdc51b2ad189691802d8.png')\nnx = len(picture) # number of pixels for width\nny = len(picture[0]) # number of pixels for height","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"65cf456543f72eee7e6f7016ff2e1a1ad01ef10c"},"cell_type":"markdown","source":""},{"metadata":{"trusted":false,"_uuid":"d4151559dfe49cf4ac652d4748194eb98397b98c"},"cell_type":"code","source":"'''\n# For colorful images, PCA is tried to reduce RGB images to white&black images\n# Colorful images actually need more careful and complicated treatment\npicture = picture[:,:,0:3]\npca = PCA(n_components = 3)\npicture = picture.reshape(nx*ny,3)\npicture = pca.fit_transform(picture)\npicture = picture.reshape(nx,ny,3)\n'''\npicture = picture[:,:,0] # In this case the image is just white&black","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"8f4b46bd5978747c18552ebe0d2bd698f81fe50a"},"cell_type":"code","source":"# Normalize images\npicture_stat = np.copy(picture)\npicture_stat = picture_stat.reshape(1,nx*ny)\npicture_stat.sort()\n# Pixels are normalized according to the 100th largest and smallest numbers instead of the largest and smallest to avoid outlier pixels\npicture = (picture-picture_stat[0,100])/(picture_stat[0,nx*ny-100]-picture_stat[0,100])\npicture_stat = picture.reshape(1,nx*ny)\n# Images can have bright background from fluorescent light or dark background\n# Pixels are reversed if background is bright\naverage = np.mean(picture_stat)\nif average > 0.5:\n picture = 1-picture","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b5c594a6aca27de34d7ab32b425856ea348fd770"},"cell_type":"markdown","source":"### 1. Apply adaptive threshold"},{"metadata":{"_uuid":"5d9bd0c388c8b5349a05f773630aa0d778fd8fd1"},"cell_type":"markdown","source":"Images can have uneven background. Background in the green triangle in the following image is brighter than other parts. "},{"metadata":{"_uuid":"715ee83591c58de7a55abb478fb9948a72aecb64"},"cell_type":"markdown","source":""},{"metadata":{"_uuid":"dc117daa75c645ee17f4c7463459b19279dce0b0"},"cell_type":"markdown","source":"Then if one tries to keep pixels brighter than a threshold and abandon pixels darker than the threshold, some noisy pixels in the bottom right corner will be kept and some useful pixels in the top left part will be abandoned."},{"metadata":{"trusted":false,"_uuid":"e2500de0966d2b792f8c586dc27a15c25caa219f"},"cell_type":"code","source":"picture_copy = np.copy(picture)\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n if picture_copy[x,y] > 0.11:\n picture_copy[x,y] = 1\n else:\n picture_copy[x,y] = 0\nplt.imshow(picture_copy)\n_ = plt.title('Constant Threshold')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"2830d24a79b32701cc37a6860b3b2d85cbb9c98f"},"cell_type":"code","source":"# This function smoothens the picture by assigning a pixel average brightness within a square of size s around it.\n# Alternatively, a gaussian kernel can be used\ndef blurPicture(picture_copy,s):\n dp = np.zeros((nx+1,ny+1)) # numpy array dp[x,y] stores the summation of brightness of all pixels in the image to the left and top of point (x,y)\n for x in np.arange(0,nx):\n dp[x,0] = 0\n for y in np.arange(1,ny):\n dp[0,y] = 0\n for x in np.arange(1,nx+1):\n for y in np.arange(1,ny+1):\n dp[x,y] = dp[x-1,y]+dp[x,y-1]-dp[x-1,y-1]+picture_copy[x-1,y-1]\n for x in np.arange(0,nx):\n for y in np.arange(0,ny):\n left = max(x-s,0)\n right = min(x+s+1,nx)\n upper = max(y-s,0)\n lower = min(y+s+1,ny)\n picture_copy[x,y] = (dp[right,lower]-dp[left,lower]-dp[right,upper]+dp[left,upper])/(right-left)/(lower-upper)\n \n\n# Images can have uneven illumination.\n# Therefore, when certain threshold is applied to decide whether a pixel should be kept or not, the threshold should be adapted \n# according to local brightness.\n# Function blurPicture can be used to calculate local brightness.\nlocal_brightness = np.copy(picture)\nblurPicture(local_brightness,math.floor(min(nx,ny)/2)) # Local brightness is defined as the average brightness of pixels in a square of size min((nx,ny)/2) around the point","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"f93b74710a43cb713ae62eedbcf6509359da8c9c"},"cell_type":"code","source":"picture_copy = np.copy(picture)\naverage = np.mean(picture_copy.reshape(1,nx*ny))\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n if picture_copy[x,y] > 0.11*local_brightness[x,y]/average:\n picture_copy[x,y] = 1\n else:\n picture_copy[x,y] = 0\nplt.imshow(picture_copy)\n_ = plt.title('Adaptive Threshold')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"181e7457733fb45c4e384ebbe7a4343841366377"},"cell_type":"markdown","source":"## 2. Choose threshold according to brightness histogram"},{"metadata":{"trusted":false,"_uuid":"eab00524930818c5be28d6ee8efa2a745805fe68"},"cell_type":"code","source":"picture_stat = np.copy(picture)\npicture_stat = picture_stat.reshape(1,nx*ny)\npicture_stat.sort()\nedges = []\nfor x in np.arange(0,100):\n edges.append(x/100)\ndistribution = np.histogram(picture_stat,bins = edges)\nsns.regplot(distribution[1][:-1],distribution[0],fit_reg=False)\n_ = plt.title('brightness histogram')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"a0bea8ad00e2c6b1b66960365f5f0d7950372bdd"},"cell_type":"markdown","source":"Graph above shows the distribution of pixel brightness. The seemingly gaussian distribution at low brightness is background noise. So the threshold can be set according to this distribution"},{"metadata":{"trusted":false,"_uuid":"59b3b6600607118968558139c5a1c43f8a6659f3"},"cell_type":"code","source":"origin_picture = np.copy(picture)\nlargest_value = 0\nsecond_value = 0\nlargest_pos = 0\nsecond_pos = 0\nfor x in np.arange(0,len(distribution[0])):\n if distribution[0][x] > largest_value:\n second_value = largest_value\n second_pos = largest_pos\n largest_value = distribution[0][x]\n largest_pos = distribution[1][x]\n elif distribution[0][x] > second_value:\n second_value = distribution[0][x]\n second_pos = distribution[1][x]\nhalf_height = 0\nfor x in np.arange(len(distribution[0])-1,-1,-1):\n if distribution[0][x] > 0.5*largest_value:\n half_height = distribution[1][x]\n break\n# Threshold is set to the larger one of 0.1 or half height position.\nvalve = max(half_height,0.25)\naverage = np.mean(picture_stat)\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n if picture[x,y] > max(valve*local_brightness[x,y]/average,valve):\n picture[x,y] = 1\n else:\n picture[x,y] = 0\n \npicture_helper = np.copy(picture) # picture_helper is used to deal with the problem caused by nuclei with large contrast described in part 4. \n_ = plt.imshow(picture)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"090f104d95cf66662fae181f83c3ad6a7c19126d"},"cell_type":"markdown","source":"## 3. Identify boundaries between adjacent nuclei using gradient"},{"metadata":{"_uuid":"328b3447c776dafdff16b86dd74dfa2bbc855599"},"cell_type":"markdown","source":"Different approaches can be used to identify boundaries between adjacent nuclei. For example, when two nuclei overlap, sharp turns are highly possible to exist when one walks along the overall boundary of the two. However, an easy method is to utilize the gradient of the image. "},{"metadata":{"trusted":false,"_uuid":"0795f67fa74739a093e5dc7855b17e329725c931"},"cell_type":"code","source":"blurPicture(origin_picture,1)\n# Learnt from some experience, extracting gradient from the original image works better than from the smoothened one \n# px stores gradient in horizontal direction. py stores gradient in vertical direction.\n# pz is the magnitude of gradient\npx = np.zeros((nx,ny))\nfor x in np.arange(0,nx-1):\n for y in np.arange(0,ny):\n px[x,y] = origin_picture[x+1,y]-origin_picture[x,y]\n\npy = np.zeros((nx,ny))\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny-1):\n py[x,y] = origin_picture[x,y+1]-origin_picture[x,y]\n\npz = np.zeros((nx,ny))\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n pz[x,y] = np.sqrt(np.power(px[x,y],2)+np.power(py[x,y],2))\n \nplt.imshow(pz)\n_ = plt.title('Gradient Image')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b0771888e996ed82744cab81c04337dcded7e5bc"},"cell_type":"markdown","source":"From the graph above, boundaries between nuclei are clearly seen."},{"metadata":{"trusted":false,"_uuid":"9325edaead0c8809a3a268c9ab10d6b3707673a3"},"cell_type":"code","source":"# Normalize the gradient image and reverse it if necessary.\npz_stat = np.copy(pz)\npz_stat = pz_stat.reshape(1,nx*ny)\npz_stat.sort()\npz = (pz-pz_stat[0,100])/(pz_stat[0,nx*ny-100]-pz_stat[0,100])\npz_stat = pz.reshape(1,nx*ny)\npz_average = np.mean(pz_stat)\nif pz_average > 0.5:\n pz = 1-pz\npz_stat = np.copy(pz)\npz_stat = pz_stat.reshape(1,nx*ny)\npz_stat.sort()\n\n# We define threshold for the gradient image to identify boundaries between nuclei.\npz_valve = pz_stat[0,0]+0.07*(pz_stat[0,nx*ny-1]-pz_stat[0,0])\npz_local_brightness = np.copy(pz)\n# Suppose thicknesses of boundaries are close, brighter nuclei should have boundaries with larger gradient.\n# Therefore, an adaptive threshold changing with brightness of the original image should be used.\n# In order to handle boundaries between a bright and a dark nuclei, a smoothened original image is used.\nblurPicture(pz_local_brightness,5)\naverage = np.mean(pz_local_brightness.reshape(1,nx*ny))\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n if picture[x,y] == 1 and pz[x,y] > pz_valve*pz_local_brightness[x,y]/average:\n pz[x,y] = 1\n else:\n pz[x,y] = 0\nplt.imshow(pz)\n_ = plt.title('Boundaries between Nuclei')","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"ba7f0df96873640e13578e24f00b1ca47b35cdad"},"cell_type":"code","source":"for x in np.arange(0,nx):\n for y in np.arange(0,ny):\n picture[x,y] = max(picture[x,y]-pz[x,y],0)\n\n_ = plt.imshow(picture)","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"1014c54363047777c0e5b64629e2a66265bd3da8"},"cell_type":"markdown","source":"After boundaries are subtracted, patches of nuclei are separated to individual nuclei. Some disadvantages of this method are also exposed. For example, small nuclei may totally disappear after boundaries are subtracted."},{"metadata":{"_uuid":"35f06d195d3a9bf001695b08fe66f7d57326873b"},"cell_type":"markdown","source":""},{"metadata":{"_uuid":"7c151b1e6a6009657dedc79196cade6591232b5f"},"cell_type":"markdown","source":"### 4. Handle nuclei with large contrast"},{"metadata":{"_uuid":"3c780cd2ce50a97992d698e45108d1822c94d037"},"cell_type":"markdown","source":"Another drawback of this method is that when a nucleus itself has very large contrast. After boundaries are subtracted, the nucleus may be splitted to several parts as shown in the following graph. The nucleus is a single piece before boundaries are subtracted."},{"metadata":{"_uuid":"eae032ec5e7a2a28ef316f36e838dcc5228f09db"},"cell_type":"markdown","source":""},{"metadata":{"_uuid":"0eb02b954c808bb649e51a2b1488ebaa3ba1f80f"},"cell_type":"markdown","source":"One observation is that after boundaries are subtracted, these nuclei with large contrast will have lots of alternate 0,1 between adjacent pixels. If real boundaries are subtracted, the resulting image will be less scattered as shown in the following graph. This can be used as a criteria to destinguish real or fake boundaries."},{"metadata":{"_uuid":"923ab4e58f68d0196bf15747c4ff5175049e2d7f"},"cell_type":"markdown","source":""},{"metadata":{"trusted":false,"_uuid":"20e28999ab80c64003d329a57e29b0a87748b991"},"cell_type":"code","source":"# This function find connected areas in an image\ndef labelNuclei(picture_copy,x,y):\n if x < 0 or x > nx-1 or y < 0 or y > ny-1:\n return\n if picture_copy[x,y] != 1:\n return\n picture_copy[x,y] = num\n # Depth first search is not allowed to be performed along diagonal directions since sometimes after boundaries are substracted, some adjacent nuclei are still\n # connected through one pixel along diagonal directions.\n labelNuclei(picture_copy,x+1,y)\n labelNuclei(picture_copy,x-1,y)\n labelNuclei(picture_copy,x,y+1)\n labelNuclei(picture_copy,x,y-1)\n\n# Label connected areas in an image with numbers\nnum = 1\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n if picture_helper[x,y] == 1:\n num = num+1\n labelNuclei(picture_helper,x,y)","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"ec665196e4f248171649ae2712a909b91dc19617"},"cell_type":"code","source":"directions = [[1,1],[1,0],[1,-1],[0,1],[0,0],[0,-1],[-1,1],[-1,0],[-1,-1]]\n\n# This function returns numbers of pixels with the other value(0 vs 1) around point (x,y)\ndef getScatter(picture_helper_copy,picture_copy,x,y):\n result = 0\n for z in np.arange(0,9):\n x_new = x+directions[z][0]\n y_new = y+ directions[z][1]\n if x_new >= 0 and x_new < nx and y_new >= 0 and y_new < ny:\n # If (x,y) is on the boundary of a connected area, 0 is returned. This is to reduce the effect of the size of a nucleus.\n if picture_helper_copy[x_new,y_new] == 0:\n return 0\n if picture_copy[x,y] != picture_copy[x_new,y_new]:\n result = result+1\n return result\n\n# Check whether some connected areas in picture_helper(before boundaries are subtracted) become very scattered in picture(after boundaries are subtracted).\n# If yes, recombine the area to a single piece in picture.\nscatter_valve = 1.7\nneed_union = set()\nscatter = [0]*(num+1) # stores numbers of alternative 0,1 in picture within each connected area in picture_helper.\nnuclei_sizes = [0]*(num+1) # stores sizes of connected areas in picture_helper.\n# Sometimes a nucleus with large contrast will just become very sparse instead of scattered after boundaries are subtracted.\nzero_cnt = [0]*(num+1) # countes numbers of 0 in picture within each connected area in picture_helper.\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n label = int(picture_helper[x,y])\n if label != 0:\n nuclei_sizes[label] = nuclei_sizes[label]+1\n scatter[label] = scatter[label]+getScatter(picture_helper,picture,x,y)\n if picture[x,y] == 0:\n zero_cnt[label] = zero_cnt[label]+1\nmax_size = np.max(nuclei_sizes)\nfor x in np.arange(2,num+1):\n if scatter[x]/nuclei_sizes[x] > scatter_valve and nuclei_sizes[x] > 0.1*max_size: # If the size of a nucleus is too small, it can be background noise.\n need_union.add(x)\n elif zero_cnt[x]/nuclei_sizes[x] > 0.7:\n need_union.add(x)\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n if picture_helper[x,y] in need_union:\n picture[x,y] = 1","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"dc57e8fcf12ff9efd83fe0a0eb84b3411ce5119d"},"cell_type":"code","source":"# It's now time to label individual nucleus.\nnum = 1\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n if picture[x,y] == 1:\n num = num+1\n labelNuclei(picture,x,y)\n \nplt.figure(figsize=(20,7),dpi=80) \nplt.imshow(picture)\n_ = plt.title('Nuclei Tagging')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"c8829525bf8cf6f588ca52eb6871daa805752708"},"cell_type":"markdown","source":"## 5. Post processing "},{"metadata":{"_uuid":"a730f393664af7705fa9d40652f3639e3a21ee18"},"cell_type":"markdown","source":"### a. Grow nuclei"},{"metadata":{"_uuid":"0c57e49214f5107d006fb565f9ed827432fb08ee"},"cell_type":"markdown","source":"Since after boundaries are subtracted, nuclei may become too small, to allow nuclei to grow, each pixel is assigned the majority of nuclei label around it."},{"metadata":{"trusted":false,"_uuid":"cfd8cf835257cfec7963cd3e5b57cd41a93cf9c1"},"cell_type":"code","source":"# This function finds the majority vote of a pixels' neighbors.\ndef major(final,x,y):\n cnt = dict()\n for z in np.arange(0,9):\n x_1 = x+directions[z][0]\n y_1 = y+directions[z][1]\n if x_1 > 0 and x_1 < nx and y_1 > 0 and y_1 < ny:\n if final[x_1,y_1] not in cnt:\n cnt.update({final[x_1,y_1]:1})\n else:\n cnt[final[x_1,y_1]] = cnt[final[x_1,y_1]]+1\n majority = 0\n major_cnt = 0\n for x in cnt.items():\n if x[0] != 0:\n if x[1] > major_cnt:\n majority = x[0]\n major_cnt = x[1]\n if major_cnt <= 2:\n return 0\n return majority\n\n# Nuclei are allowed to grow in the mentioned manner for two iterations.\npicture_helper = np.zeros((nx,ny))\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n picture_helper[x,y] = major(picture,x,y)\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n picture[x,y] = major(picture_helper,x,y)\n\nplt.figure(figsize=(20,7),dpi=80) \nplt.imshow(picture)\n_ = plt.title('After Growth')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"e670b70ebcd9bfae2d80b5f59dc24b232116b457"},"cell_type":"markdown","source":"### b. Adjust nuclei"},{"metadata":{"trusted":false,"_uuid":"a07abbc6cf6ac7a598b17f17251b7ab411a9c3bf"},"cell_type":"code","source":"# Calculate average brightness of each nucleus\nnuclei_brightness = dict()\nnuclei_sizes = dict()\nfor x in range(0,nx):\n for y in range(0,ny):\n if picture[x,y] not in nuclei_sizes:\n nuclei_sizes.update({picture[x,y]:1})\n nuclei_brightness.update({picture[x,y]:origin_picture[x,y]})\n else:\n nuclei_sizes[picture[x,y]] = nuclei_sizes[picture[x,y]]+1\n nuclei_brightness[picture[x,y]] = nuclei_brightness[picture[x,y]]+origin_picture[x,y]\nfor z in nuclei_brightness.items():\n if z[0] != 0:\n nuclei_brightness[z[0]] = nuclei_brightness[z[0]]/nuclei_sizes[z[0]]\n\n# According to how the score is calculated by Kaggle, adjustment is needed to make sure masks covers exactly nuclei.\n# After the above step, pixels are adjusted according to following criteria.\n# 1. If a pixel is darker than the brightness of all surrounding nuclei times a factor, the pixel is removed.\n# 2. Otherwise a pixel is assigned to the nucleus with the closest brightness around it.\ndef adjustNuclei(picture_copy,origin_pic,x,y,local_brightness,average,nuclei_brightness):\n factor = 0.7\n cnt = dict()\n for x_1 in range(x-1,x+2):\n for y_1 in range(y-1,y+2):\n if x_1 > 0 and x_1 < nx and y_1 > 0 and y_1 < ny:\n if picture_copy[x_1,y_1] not in cnt:\n cnt.update({picture_copy[x_1,y_1]:1})\n else:\n cnt[picture_copy[x_1,y_1]] = cnt[picture_copy[x_1,y_1]]+1\n darkest_nucleus = 5\n majority = 0\n major_cnt = 0\n for z in cnt.items():\n if z[0] != 0:\n if nuclei_brightness[z[0]] < darkest_nucleus:\n darkest_nucleus = nuclei_brightness[z[0]]\n if z[1] > major_cnt:\n major_cnt = z[1]\n majority = z[0]\n if darkest_nucleus == 5 or origin_pic[x,y] <= min(factor*darkest_nucleus*local_brightness[x,y]/average,factor*darkest_nucleus):\n return 0\n for z in cnt.items():\n if z[0] != 0 and z[0]!=majority:\n if abs(nuclei_brightness[z[0]]-nuclei_brightness[majority]) > 0.2:\n possibility_1 = abs(origin_pic[x,y]-nuclei_brightness[majority])\n possibility_2 = abs(origin_pic[x,y]-nuclei_brightness[z[0]])\n if possibility_1 > possibility_2:\n return z[0]\n return majority\n\npicture_helper = np.zeros((nx,ny))\nfor n in range(0,2):\n for x in range(0,nx):\n for y in range(0,ny):\n picture_helper[x,y] = adjustNuclei(picture,origin_picture,x,y,local_brightness,average,nuclei_brightness)\n for x in range(0,nx):\n for y in range(0,ny):\n picture[x,y] = picture_helper[x,y]\n \nplt.figure(figsize=(20,7),dpi=80) \nplt.imshow(picture)\n_ = plt.title('After Adjustment')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"daff1b6799e7b10af8f6f965030d36ff48cd0da0"},"cell_type":"markdown","source":"### c. Filtering"},{"metadata":{"_uuid":"cdc93da6626d088bf8f718c2d78528f5dc69a2c4"},"cell_type":"markdown","source":"Two criteria are applied to filter the final result.\n2. A nucleus is abandoned if it's both darker than a factor times the brightness of the brightest nucleus and smaller than a factor times the size of the largest nucleus.\n3. A nucleus is abandoned if it's not close to a circle. Because of surface tension, a nucleus should be close to a circle. Correlation of a nucleus' x,y coordinates is used as an easy evaluation to decide whether it's close to circle."},{"metadata":{"trusted":false,"_uuid":"924697bf95b50b6f4dfb704a49385ee40fa381b1"},"cell_type":"code","source":"nuclei_brightness = dict() # stores average brightness of each nucleus.\nnuclei_sizes = dict() # stores sizes of each nucleus.\nnuclei_pixels_x = dict() # stores an array of x coordinates of each nucleus\nnuclei_pixels_y = dict() # stores an array of y coordinates of each nucleus\nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n if picture[x,y] not in nuclei_sizes:\n nuclei_sizes.update({picture[x,y]:1})\n nuclei_brightness.update({picture[x,y]:origin_picture[x,y]})\n nuclei_pixels_x.update({picture[x,y]:[x]})\n nuclei_pixels_y.update({picture[x,y]:[y]})\n else:\n nuclei_sizes[picture[x,y]] = nuclei_sizes[picture[x,y]]+1\n nuclei_brightness[picture[x,y]] = nuclei_brightness[picture[x,y]]+origin_picture[x,y]\n nuclei_pixels_x[picture[x,y]].append(x)\n nuclei_pixels_y[picture[x,y]].append(y)\n\nfor z in nuclei_brightness.items():\n if z[0] != 0:\n nuclei_brightness[z[0]] = nuclei_brightness[z[0]]/nuclei_sizes[z[0]]","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"70229f12e2ef9a1cf7b2bf16dd28474c4f50f6e0"},"cell_type":"code","source":"need_remove = set()\n# Remove all nuclei smaller than 5 pixels.\nfor m in nuclei_sizes.items():\n if m[0] != 0 and m[1] < 5:\n need_remove.add(m[0])\n \n# Remove all nuclei darker than 0.2 times the brightness of the brightest nucleus and smaller than 0.2 times the size of the largest nucleus.\nlargest = 0\nbrightest = 0\nfor m in nuclei_sizes.items():\n if m[0] != 0:\n if m[1] > largest:\n largest = m[1]\n if nuclei_brightness[m[0]] > brightest:\n brightest = nuclei_brightness[m[0]]\n\nfor m in nuclei_sizes.items():\n if m[0] != 0 and m[1] < 0.2*largest and nuclei_brightness[m[0]] < 0.2*brightest:\n need_remove.add(m[0])\n \n# Remove all nuclei with a shape not close to a circle.\nfor m in nuclei_pixels_x.items():\n if m[0] != 0:\n covariance = np.cov(m[1],nuclei_pixels_y[m[0]])\n if abs(covariance[0,1]/covariance[0,0]/covariance[1,1]) > 0.5:\n need_remove.add(m[0])\n \nfor x in np.arange(0,nx):\n for y in np.arange(0,ny):\n if picture[x,y] in need_remove:\n picture[x,y] = 0","execution_count":null,"outputs":[]},{"metadata":{"trusted":false,"_uuid":"3e91368f3e2379ec56238769ccd86f7684999a81"},"cell_type":"code","source":"plt.figure(figsize=(20,7),dpi=80) \nplt.imshow(picture)\n_ = plt.title('Final Result')","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"b347e573ba2d065ed3b2f5d0c4720a4596b4e922"},"cell_type":"markdown","source":"And the raw image for comparation."},{"metadata":{"_uuid":"01973487f4696384b89056163f14f38ea16dc0aa"},"cell_type":"markdown","source":""},{"metadata":{"trusted":true,"_uuid":"b2964524c4b873a1857e3351b1a0f69405151322"},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}