title: Canny Edge Detection- Non Maximum Surpression Implementation
body: <p>I am trying to implement the Canny edge detection algorithm from scratch with the help of OpenCV. I am facing a problem implementing the Non-Maximum Suppression step which helps to thin the edges. </p>
  
  <p>My logic is to first compute the intensity gradient vector, then group it in either 0,45,90,135 degrees direction and then try to find local maxima. The method to find this local maxima is by making sure that current pixel is greater than succeeding and preceeding pixel in same direction. If not, I assign value of zero to this pixel. Using this logic, I'm still not able to thin the edges. I feel the error is when I'm computing the intensity gradient vector for each pixel.</p>
  
  <p>Here is my code-</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;bits/stdc++.h&gt;
  #include "opencv2/core/core.hpp"
  #include "opencv2/highgui/highgui.hpp"
  #include "opencv2/opencv.hpp"
  #include "opencv2/imgproc/imgproc.hpp"
  #include &lt;opencv2/objdetect/objdetect.hpp&gt;
  #include &lt;math.h&gt;
  
  using namespace cv; 
  using namespace std; 
  
  int main()
  {
  // Reading image
  Mat img = imread("1.jpg");
  
  // Displaying image
  //imshow("Original Image",img);
  //waitKey(0);
  
  // Converting to grayscale
  Mat img_gray,image_gray;
  cvtColor(img,image_gray,CV_RGB2GRAY);
  GaussianBlur( image_gray, img_gray, Size(15,15), 3, 3);
  
  
  // Displaying grayscale image
  imshow("Original Image",img_gray);
  waitKey(0);
  
  
  int cols = img_gray.cols;
  int rows = img_gray.rows;
  
  // Creating sobel operator in x direction
  int sobel_x[3][3] = {-1,0,1,-2,0,2,-1,0,1};
  // Creating sobel operator in y direction
  int sobel_y[3][3] = {1,2,1,0,0,0,-1,-2,-1};
  
  
  int radius = 1;
  
  // Handle border issues
  Mat _src;
  copyMakeBorder(img_gray, _src, radius, radius, radius, radius, BORDER_REFLECT101);
  
  // Create output matrix
  Mat gradient_x = img_gray.clone();
  Mat gradient_y = img_gray.clone();
  Mat gradient_f = img_gray.clone();
  Mat gradient_mag = img_gray.clone();
  
  // Conrrelation loop in x direction 
  
  // Iterate on image 
  for (int r = radius; r &lt; _src.rows - radius; ++r)
  {
      for (int c = radius; c &lt; _src.cols - radius; ++c)
      {
          int s = 0;
  
          // Iterate on kernel
          for (int i = -radius; i &lt;= radius; ++i)
          {
              for (int j = -radius; j &lt;= radius; ++j)
              {
                  s += _src.at&lt;uchar&gt;(r + i, c + j) * sobel_x[i + radius][j + radius];
              }
          }
          gradient_x.at&lt;uchar&gt;(r - radius, c - radius) = s/8;
  
          /*if(s&gt;200)
              gradient.at&lt;uchar&gt;(r - radius, c - radius) = 255;
          else
              gradient.at&lt;uchar&gt;(r - radius, c - radius) = 0;
          */    
      }
  }
  
  // Conrrelation loop in y direction 
  
  // Iterate on image 
  for (int r = radius; r &lt; _src.rows - radius; ++r)
  {
      for (int c = radius; c &lt; _src.cols - radius; ++c)
      {
          int s = 0;
  
          // Iterate on kernel
          for (int i = -radius; i &lt;= radius; ++i)
          {
              for (int j = -radius; j &lt;= radius; ++j)
              {
                  s += _src.at&lt;uchar&gt;(r + i, c + j) * sobel_y[i + radius][j + radius];
              }
          }
          gradient_y.at&lt;uchar&gt;(r - radius, c - radius) = s/8;
  
          /*if(s&gt;200)
              gradient.at&lt;uchar&gt;(r - radius, c - radius) = 255;
          else
              gradient.at&lt;uchar&gt;(r - radius, c - radius) = 0;
          */    
      }
  }
  
  ///cout&lt;&lt;endl&lt;&lt;"max:"&lt;&lt;max;
  //cout&lt;&lt;img_gray.rows;
  //cout&lt;&lt;endl&lt;&lt;_src.rows;
  cout&lt;&lt;endl&lt;&lt;gradient_x.rows;
  cout&lt;&lt;endl&lt;&lt;gradient_y.rows;
  cout&lt;&lt;endl&lt;&lt;gradient_f.rows&lt;&lt;gradient_f.cols;   
  
  
  //Calculating gradient magnitude
  for(int i=0; i&lt;gradient_mag.rows; i++)
  {
      for(int j=0; j&lt;gradient_mag.cols; j++)
      {
          gradient_mag.at&lt;uchar&gt;(i,j) = sqrt( pow(gradient_x.at&lt;uchar&gt;(i,j),2) + pow(gradient_y.at&lt;uchar&gt;(i,j),2) );  
  
           if(gradient_mag.at&lt;uchar&gt;(i,j) &gt;250)
              gradient_f.at&lt;uchar&gt;(i,j) = 255;
          else
              gradient_f.at&lt;uchar&gt;(i,j) = 0;
      }
  }
  
  
  /*
  imshow("grad x",gradient_x);
  waitKey(0);
  
  imshow("grad y",gradient_y);
  waitKey(0);
  */
  
  imshow("grad magnitude",gradient_f);
  waitKey(0); 
  
  int max=0;
  
  
  // Performing Non-Maximum Surpression
  float theta; // Calculate intensity gradient vector theta=atan2(Gy,Gx);
  Mat nonMaxSupp= Mat(gradient_mag.rows-2, gradient_mag.cols-2, CV_8UC1); //CV_8UC1 is 8-bit single channel image i.e grayscale
  
  for(int i=1; i&lt;gradient_x.rows-1; i++)
  {
      for(int j=1; j&lt;gradient_x.cols-1; j++)
      {
          //if(gradient_x.at&lt;uchar&gt;(i,j) ==0) //Arctan Fix
           //   theta = 90;
          //else
              theta = atan2(gradient_y.at&lt;uchar&gt;(i,j),gradient_x.at&lt;uchar&gt;(i,j))*(180/3.14);
              //theta = atan(gradient_y.at&lt;uchar&gt;(i,j)/gradient_x.at&lt;uchar&gt;(i,j))*(180/3.14);
          //cout&lt;&lt;theta&lt;&lt;endl;
          //if(theta&gt;max)
          //    max=theta;
          nonMaxSupp.at&lt;uchar&gt;(i-1, j-1) = gradient_mag.at&lt;uchar&gt;(i,j);
  
          // For horizontal edge
          if(((-22.5 &lt; theta) &amp;&amp; (theta &lt;= 22.5)) || ((157.5 &lt; theta) &amp;&amp; (theta &lt;= -157.5)))
          {
              if ((gradient_mag.at&lt;uchar&gt;(i,j) &lt; gradient_mag.at&lt;uchar&gt;(i,j+1)) || (gradient_mag.at&lt;uchar&gt;(i,j) &lt; gradient_mag.at&lt;uchar&gt;(i,j-1)))
                  nonMaxSupp.at&lt;uchar&gt;(i-1, j-1) = 0;
          }
  
          //For vertical edge
          if (((-112.5 &lt; theta) &amp;&amp; (theta &lt;= -67.5)) || ((67.5 &lt; theta) &amp;&amp; (theta &lt;= 112.5)))
          {
              if ((gradient_mag.at&lt;uchar&gt;(i,j) &lt; gradient_mag.at&lt;uchar&gt;(i+1,j)) || (gradient_mag.at&lt;uchar&gt;(i,j) &lt; gradient_mag.at&lt;uchar&gt;(i-1,j)))
                  nonMaxSupp.at&lt;uchar&gt;(i-1, j-1) = 0;
          }
  
          // For 135 degree or -45 degree edge
          if (((-67.5 &lt; theta) &amp;&amp; (theta &lt;= -22.5)) || ((112.5 &lt; theta) &amp;&amp; (theta &lt;= 157.5)))
          {
              if ((gradient_mag.at&lt;uchar&gt;(i,j) &lt; gradient_mag.at&lt;uchar&gt;(i-1,j+1)) || (gradient_mag.at&lt;uchar&gt;(i,j) &lt; gradient_mag.at&lt;uchar&gt;(i+1,j-1)))
                  nonMaxSupp.at&lt;uchar&gt;(i-1, j-1) = 0;
          }
  
          // For 45 Degree Edge
          if (((-157.5 &lt; theta) &amp;&amp; (theta &lt;= -112.5)) || ((22.5 &lt; theta) &amp;&amp; (theta &lt;= 67.5)))
          {
              if ((gradient_mag.at&lt;uchar&gt;(i,j) &lt; gradient_mag.at&lt;uchar&gt;(i+1,j+1)) || (gradient_mag.at&lt;uchar&gt;(i,j) &lt; gradient_mag.at&lt;uchar&gt;(i-1,j-1)))
                  nonMaxSupp.at&lt;uchar&gt;(i-1, j-1) = 0;
          }
  
      }
  
  }
  //cout&lt;&lt;endl&lt;&lt;"max"&lt;&lt;max;
  imshow("Non-Maximum Surpression",nonMaxSupp);
  waitKey(0);
  
  return 0;
  </code></pre>
  
  <p>}</p>