{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"

Introduction

\n","metadata":{}},{"cell_type":"markdown","source":"

Greeting fellow Kagglers,

\n

Most of us are familiar with graph theory; some understand it better, and some have a general idea of this interesting field of mathematics.

\n

Looking at the structure of the chemical compounds in the images, it is hard not to see a graph theory graph; a question arising goes as following: given the same formula represented as a graph, can the extracted graph attributes provide features that can improve the model that solves this problem, e.g. if I have an image of a chemical compound and I know that the graph the represent the compound has a maximum degree of some value and a minimum degree of another value, I know the chromatic number of the graph and the size of the largest clique or independent set, can such features included in a model improve performance or accuracy?

\n


\n

Well, it is a hard question to answer, and even before we look into such a question, how can we even create a graph out of an image? That's exactly what I was thinking in the days working on this notebook. 

\n

First, we need to detect the vertices, the points between which we will connect our edges, and the following section dedicated to creating a pipeline that will get us to that first step, knowing how much and where the graph edges are located.

","metadata":{}},{"cell_type":"markdown","source":"

Note: for those of you not familiar with graph theory an introduction to the basics of graph theory can be reviewed in this notebook

","metadata":{}},{"cell_type":"markdown","source":"

Libraries And Utilities

\n","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom collections import Counter, defaultdict\nimport cv2, os\nimport skimage.io as io\nfrom tqdm.auto import tqdm\ntqdm.pandas()\nimport matplotlib.pyplot as plt\nimport matplotlib as mpl\nimport seaborn as sns\nfrom scipy.spatial import distance_matrix\nimport itertools\nplt.rc('figure',figsize=(18,11))\nplt.rc('image',cmap='Blues')\nsns.set_context('paper',font_scale=2)\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Loading The Data

\n","metadata":{}},{"cell_type":"code","source":"labels = pd.read_csv('../input/bms-molecular-translation/train_labels.csv')\nss = pd.read_csv('../input/bms-molecular-translation/sample_submission.csv', index_col = 0)\n\nlabels['path'] = labels['image_id'].progress_apply(\n lambda x: \"../input/bms-molecular-translation/train/{}/{}/{}/{}.png\".format(\n x[0], x[1], x[2], x))\nlabels.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Going Through The Logic

\n","metadata":{}},{"cell_type":"markdown","source":"

Understanding Harris Corrner

\n","metadata":{}},{"cell_type":"markdown","source":"

We will use the Harris Corner Detection algorithms to extract the corners i.e "the significant" pixel's which in our case is the intersection or the begining of every straight line.

\n


\n

Harris Corners Detector Steps

\n

1. Compute image gradients: Gx, Gy

\n

2. Compute products: Gx*Gx, Gx*Gy, Gy*Gy

\n

3. Filter products with a Gaussian window

\n

4. For each pixel (i,j) define the matrix M

\n

5. For each pixel compute the score R

\n

6. Threshold R, and perform non-maxima suppression

\n


\n


","metadata":{}},{"cell_type":"markdown","source":"\n

1) Compute image gradients: Gx, Gy

\n","metadata":{}},{"cell_type":"code","source":"def calculate_image_gradients(img):\n #Derivative Kernels\n Kx = np.array([[-1,0,1],[-2,0,2],[-1,0,1]]) \n Ky = Kx.T\n #Calculate The Derivative with respect to the x axis\n Gx = cv2.filter2D(img,-1,Kx)\n #Calculate The Derivative with respect to the y axis\n Gy = cv2.filter2D(img,-1,Ky)\n \n return Gx,Gy","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"img = np.zeros((21,21), dtype=np.float32)\nimg[5:-5,5:-5]=200\n\nGx,Gy = calculate_image_gradients(img)\n\nplt.subplot(2,3,3)\nplt.title(\"Gx\",fontsize=16,fontweight='bold')\nplt.imshow(Gx,cmap='jet')\nplt.colorbar()\nplt.subplot(2,3,2)\nplt.title(\"Gy\",fontsize=16,fontweight='bold')\nplt.imshow(Gy,cmap='jet')\nplt.subplot(2,3,1)\nplt.title(\"Original Image\",fontsize=16,fontweight='bold')\nplt.imshow(img,cmap='jet')\nplt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n

2) Compute the second order moments: ($Gx\\cdot Gx$, $Gx\\cdot Gy$, $Gy\\cdot Gy$)

\n","metadata":{}},{"cell_type":"code","source":"titles=[r\"$Gx^{2}$\",r\"$Gy^{2}$\",r\"$Gy\\cdot Gx}$\"]\n\nimg = np.zeros((21,21), dtype=np.float32)\nimg[5:-5,5:-5]=200\nGx,Gy = calculate_image_gradients(img)\nmoments = [Gx**2,Gy**2,Gx*Gy]\n\nimg = np.zeros((21,21), dtype=np.float32)\nimg[5:-5,5:-5]=200\n\nfor i in range(0,3):\n plt.subplot(2,3,i+1)\n plt.title(titles[i],fontsize=16,fontweight='bold')\n plt.imshow(moments[i],cmap='jet')\n \nplt.colorbar()\nplt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n

3) Filter products with a Gaussian window

\n","metadata":{}},{"cell_type":"code","source":"titles=[r\"$Gx^{2}$\",r\"$Gy^{2}$\",r\"$Gy\\cdot Gx}$\"]\n\nimg = np.zeros((21,21), dtype=np.float32)\nimg[5:-5,5:-5]=200\nGx,Gy = calculate_image_gradients(img)\nmoments = [Gx**2,Gx*Gy,Gy**2]\nmoments = [cv2.GaussianBlur(moments[i],(3,3),2) for i in range(0,3)]\n\nimg = np.zeros((21,21), dtype=np.float32)\nimg[5:-5,5:-5]=200\n\nfor i in range(0,3):\n plt.subplot(2,3,i+1)\n plt.title('Filtered ' +titles[i],fontsize=16,fontweight='bold')\n plt.imshow(moments[i],cmap='jet')\n\nplt.colorbar()\nplt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n

4) For each pixel $(i,j)$ define the matrix $M$

\n","metadata":{}},{"cell_type":"markdown","source":"For each (i,j) define a 2x2 moments-matrix: $M = \\begin{bmatrix}\n \\sum{G_x^2}&\\sum{G_xG_y} \\\\ \n \\sum{G_xG_y}&\\sum{G_y^2} \n\\end{bmatrix}$\n\nWe will define the Harris Score as : $R = det(M) - \\alpha{[tr(M)]^2}$\nWhere det is the determinante operator and tr is the trace operator","metadata":{}},{"cell_type":"code","source":"img = np.zeros((21,21), dtype=np.float32)\nimg[5:-5,5:-5]=200\n\nplt.subplot(2,2,1)\nplt.title('Original Image',weight='bold')\nplt.imshow(img,cmap='jet')\nplt.subplot(2,2,2)\nGx,Gy =calculate_image_gradients(img)\nM = [Gx**2,Gx*Gy,Gy**2]\nM = [cv2.GaussianBlur(Prod,(3,3),2) for Prod in M]\ndetM = M[0]*M[2]-M[1]**2\ntraceM = M[0]+M[2]\nR_scores = detM-0.06*(traceM**2)\nplt.title('R Score Matrix',weight='bold')\nplt.imshow(R_scores,cmap='jet')\nplt.colorbar()\nplt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n

6) Threshold $R$, and perform non-maxima suppression

\n","metadata":{}},{"cell_type":"code","source":"def get_R_scores(Gx,Gy,Sigma=2,alpha=0.06,tsh=0.35):\n M = [Gx**2,Gx*Gy,Gy**2]\n M = [cv2.GaussianBlur(Prod,(3,3),Sigma) for Prod in M]\n detM = M[0]*M[2]-M[1]**2\n traceM = M[0]+M[2]\n R_scores = detM-alpha*(traceM**2)\n cross_kernel = np.ones((3,3),np.uint8)\n #cross_kernel[1,:]=1\n #cross_kernel[:,1]=1\n R_dilate = cv2.dilate(R_scores,cross_kernel)\n R_th = R_scores > R_scores.max()*tsh; # threshold\n R_nms = R_scores >= R_dilate;# NMS\n R_final = R_th * R_nms # threshold and NMS\n #plt.imshow(R_scores,cmap='jet')\n return R_final.astype(np.int)\n\nimg = np.zeros((21,21), dtype=np.float32)\nimg[5:-5,5:-5]=200\n\nplt.subplot(2,2,1)\nplt.title('Original Image',weight='bold')\nplt.imshow(img,cmap='jet')\nplt.subplot(2,2,2)\nplt.title('R Score Matrix after NMS',weight='bold')\nGx,Gy =calculate_image_gradients(img)\nplt.imshow(get_R_scores(Gx,Gy),cmap='jet')\nplt.colorbar()\nplt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Note that what we are left with the corners of our shape, the vertices,\nunfortunately, Harris corner by itself won't do us much, taking into account the amount of instability, variability, and noise in our images.\nBut now that we are familiar with the Harris corner detection algorithms, we can extend it and fit it to our needs; that is exactly what is about to be explained in the following section.

\n","metadata":{}},{"cell_type":"markdown","source":"

Understanding Our Pipeline

\n","metadata":{}},{"cell_type":"markdown","source":"

Our extension to the Harris corner algorithm requires us to find corners that may have lower R scores due to noise and the overall quality of quite many images.

\n

If we use more strict parameters in the Harris detector, i.e., the aperture and window size and the alpha value, we will get more quality and stable vertices, representing the intersection between straight lines; but we will miss a lot of vertices due to the already mentioned noise and quality.

\n

The solution implemented in the following pipeline goes as follows:

\n
    \n
  1.  Use Large Harris Detector Parameters To Get A Noisy Set Of Potential Vertecis, A Set Which Definitely Covers All Vertices.
    2. \n
  2. Preform A Modified Version Of KNN, Where Depending On A Threshold, All Vertecis That Grouped Within A Certiean Distance Leave Only The First Detected Vertex.
    3. \n
","metadata":{}},{"cell_type":"markdown","source":"

1- Use Harris Corenr with Robust Parameters

\n","metadata":{}},{"cell_type":"code","source":"def clean_noise(image,stepSize=2,windowSize=(3,3)):\n\n for y in range(0, image.shape[0], stepSize):\n for x in range(0, image.shape[1], stepSize):\n if (windowSize[0]**2 - image[y:y + windowSize[1], x:x + windowSize[0]].flatten().sum()) == 1:\n image[y:y + windowSize[1], x:x + windowSize[0]] = 1\n return image\n\nimg = plt.imread(labels.path[55])\nimg = clean_noise(img)\npimg = cv2.erode(img,np.ones((3,3)),iterations=2)\ncimg = cv2.cvtColor(img,cv2.COLOR_GRAY2RGB)\n\nplt.subplot(2,2,1)\nplt.title('Eroded Example')\nplt.imshow(pimg)\nplt.subplot(2,2,2)\n\nR = cv2.cornerHarris(np.float32(pimg),7,3,0.04)\n\nR_dilate = cv2.dilate(R, np.ones((3,3)))\nR_nms = R >= R_dilate;# NMS\n\nR_th = R > R.max()*0.1; # threshold\n\nR_final = R_th * R_nms # threshold and NMS\n[y,x] = np.nonzero(R_final)\n\nfor x,y in zip(x,y):\n cimg = cv2.circle(cimg,(x,y),3,(255,0,0))\n\nplt.title('Detected Vertices')\nplt.imshow(cimg)\nplt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

2- Filter Point Using Modified KN Clustering

\n","metadata":{}},{"cell_type":"code","source":"\ndef threshold_neighbours(distance_matrix,index,threshold):\n return np.where(distance_matrix[index,:] = R_dilate;# NMS\n R_th = R > R.max()*nms_threshold; # threshold\n R_final = R_th * R_nms # threshold and NMS\n [y,x] = np.nonzero(R_final)\n xx,yy = cluster_vertices(x,y,cluster_threshold)\n return xx,yy\n\nplt.subplot(2,2,1)\nplt.title('Detected Vertices Before Filtration')\nplt.imshow(cimg)\nplt.subplot(2,2,2)\nplt.title('Detected Vertices After Filtration')\nx,y = tag_vertices(img,apertureSize=3,cluster_threshold=8)\nccimg = cv2.cvtColor(img,cv2.COLOR_GRAY2RGB)\nfor x,y in zip(x,y):\n ccimg = cv2.circle(ccimg,(x,y),4,(255,0,0),-1)\nplt.imshow(ccimg)\n\nplt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Vertex Detection Pipeline

\n","metadata":{}},{"cell_type":"code","source":"def sliding_window(image, stepSize, windowSize):\n for y in range(0, image.shape[0], stepSize):\n for x in range(0, image.shape[1], stepSize):\n yield (x, y, image[y:y + windowSize[1], x:x + windowSize[0]])\n\n \ndef clean_noise(image,stepSize=2,windowSize=(3,3)):\n \n \"\"\"\n \n Parameters\n ----------\n image : np.array\n an image of a chemical compound\n stepSize : int\n The number of pixels the sliding window moves each step\n windowSize : tuple-(int,int)\n The size of the sliding window \n\n \"\"\"\n for y in range(0, image.shape[0], stepSize):\n for x in range(0, image.shape[1], stepSize):\n if (windowSize[0]**2 - image[y:y + windowSize[1], x:x + windowSize[0]].flatten().sum()) == 1:\n image[y:y + windowSize[1], x:x + windowSize[0]] = 1\n return image\n\n\ndef threshold_neighbours(distance_matrix,index,threshold):\n \"\"\"\n \n Parameters\n ----------\n distance_matrix : np.array\n The distnace matrix of each detected vertex with every other\n index : int\n The index of our target vertex\n threshold : float\n The maximum distance allowed before considering a vertex to close to our target i.e a neighbour\n\n \"\"\"\n return np.where(distance_matrix[index,:] = R_dilate;# NMS\n R_th = R > R.max()*nms_threshold; # threshold\n R_final = R_th * R_nms # threshold and NMS\n \n \n [y,x] = np.nonzero(R_final)\n xx,yy = cluster_vertices(x,y,cluster_threshold,minkowski_p)\n return xx,yy\n \n ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Examples and A Brief Analysis

\n","metadata":{}},{"cell_type":"code","source":"imgs = [230,5,55]\nfig,axs = plt.subplots(2,3)\nfor ix,img in enumerate(imgs):\n img = plt.imread(labels.path[img])\n img = clean_noise(img)\n pimg = cv2.erode(img,np.ones((3,3)),iterations=2)\n\n x,y = tag_vertices(pimg,apertureSize=3,blockSize=7,cluster_threshold=12,minkowski_p=2)\n ccimg = cv2.cvtColor(img,cv2.COLOR_GRAY2RGB)\n for x,y in zip(x,y):\n ccimg = cv2.circle(ccimg,(x,y),3,(255,0,0),-1)\n axs[1,ix].imshow(ccimg)\n axs[0,ix].imshow(img)\n ","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def count_nodes(img_path):\n #Read An Image\n img = plt.imread(img_path)\n #Remove Salt And Paper Noise - Not Using MedianBlur\n img = clean_noise(img)\n #Erode To Close Gaps\n img = cv2.erode(img,np.ones((3,3)),iterations=2)\n #Run Our Pipeline\n x,y = tag_vertices(img)\n return len(x)\n\nSCAN_N = 500\nvert_df = pd.DataFrame({'path':labels.loc[:SCAN_N,'path']})\nvert_df['Num_Of_Vertices'] = labels.loc[0:SCAN_N,'path'].progress_apply(count_nodes)\nvert_df['Img_Height'] = labels.loc[0:SCAN_N,'path'].progress_apply(lambda x: plt.imread(x).shape[0])\nvert_df['Img_Width'] = labels.loc[0:SCAN_N,'path'].progress_apply(lambda x: plt.imread(x).shape[1])","metadata":{"_kg_hide-output":false,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"vert_df","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.title(f'Distribution of Vertices Count in First {SCAN_N} Images - With Noise')\nsns.histplot(vert_df.Num_Of_Vertices)\nplt.grid()\nplt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy import stats\n\n\nplt.suptitle(\"Linear Relationship Between Image Scale And Number of Vertices\", fontsize=16,fontweight='bold')\nplt.subplot(2,1,1)\nslope, intercept, _,_,_ = stats.linregress(vert_df['Num_Of_Vertices'],vert_df['Img_Width'])\nplt.title(r'$Width ={} \\cdot X+{}$'.format(np.round(slope,2),np.round(intercept,2)))\nsns.regplot(x=vert_df.Num_Of_Vertices,y=vert_df.Img_Width,line_kws=dict(color='red',lw='2',ls='-.'))\nplt.grid()\nplt.subplot(2,1,2)\nslope, intercept, _,_,_ = stats.linregress(vert_df['Num_Of_Vertices'],vert_df['Img_Height'])\nplt.title(r'$Height ={} \\cdot X+{}$'.format(np.round(slope,2),np.round(intercept,2)))\nsns.regplot(x=vert_df.Num_Of_Vertices,y=vert_df.Img_Height,line_kws=dict(color='red',lw='2',ls='-.'))\nplt.grid()\nplt.tight_layout()\nplt.subplots_adjust(top=0.88)\nplt.show()","metadata":{"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"

Conclusions and Further Directions

\n","metadata":{}},{"cell_type":"markdown","source":"\n


\n\n


\n","metadata":{}}]}