title: Obtaining orientation map of fingerprint image using OpenCV
body: <p>I'm trying to implement the method of improving fingerprint images by <a href="http://www.math.tau.ac.il/~turkel/imagepapers/fingerprint.pdf" rel="nofollow noreferrer">Anil Jain</a>. As a starter, I encountered some difficulties while extracting the orientation image, and am strictly following those steps described in Section 2.4 of that paper. </p>
  
  <p>So, this is the input image:</p>
  
  <p><img src="https://i.stack.imgur.com/nwmPa.png" alt="enter image description here"></p>
  
  <p>And this is after normalization using exactly the same method as in that paper:</p>
  
  <p><img src="https://i.stack.imgur.com/vd5d9.png" alt="enter image description here"></p>
  
  <p>I'm expecting to see something like this (an example from the internet):</p>
  
  <p><img src="https://i.stack.imgur.com/k4yEd.png" alt="enter image description here"></p>
  
  <p>However, this is what I got for displaying obtained orientation matrix:</p>
  
  <p><img src="https://i.stack.imgur.com/YhN0T.png" alt="enter image description here"></p>
  
  <p>Obviously this is wrong, and it also gives non-zero values for those zero points in the original input image.</p>
  
  <p>This is the code I wrote:</p>
  
  <pre><code>cv::Mat orientation(cv::Mat inputImage)
  {
      cv::Mat orientationMat = cv::Mat::zeros(inputImage.size(), CV_8UC1);
  
      // compute gradients at each pixel
      cv::Mat grad_x, grad_y;
      cv::Sobel(inputImage, grad_x, CV_16SC1, 1, 0, 3, 1, 0, cv::BORDER_DEFAULT);
      cv::Sobel(inputImage, grad_y, CV_16SC1, 0, 1, 3, 1, 0, cv::BORDER_DEFAULT);
  
      cv::Mat Vx, Vy, theta, lowPassX, lowPassY;
      cv::Mat lowPassX2, lowPassY2;
  
      Vx = cv::Mat::zeros(inputImage.size(), inputImage.type());
      Vx.copyTo(Vy);
      Vx.copyTo(theta);
      Vx.copyTo(lowPassX);
      Vx.copyTo(lowPassY);
      Vx.copyTo(lowPassX2);
      Vx.copyTo(lowPassY2);
  
      // estimate the local orientation of each block
      int blockSize = 16;
  
      for(int i = blockSize/2; i &lt; inputImage.rows - blockSize/2; i+=blockSize)
      {    
          for(int j = blockSize / 2; j &lt; inputImage.cols - blockSize/2; j+= blockSize)
          {
              float sum1 = 0.0;
              float sum2 = 0.0;
  
              for ( int u = i - blockSize/2; u &lt; i + blockSize/2; u++)
              {
                  for( int v = j - blockSize/2; v &lt; j+blockSize/2; v++)
                  {
                      sum1 += grad_x.at&lt;float&gt;(u,v) * grad_y.at&lt;float&gt;(u,v);
                      sum2 += (grad_x.at&lt;float&gt;(u,v)*grad_x.at&lt;float&gt;(u,v)) * (grad_y.at&lt;float&gt;(u,v)*grad_y.at&lt;float&gt;(u,v));
                  }
              }
  
              Vx.at&lt;float&gt;(i,j) = sum1;
              Vy.at&lt;float&gt;(i,j) = sum2;
  
              double calc = 0.0;
  
              if(sum1 != 0 &amp;&amp; sum2 != 0)
              {
                  calc = 0.5 * atan(Vy.at&lt;float&gt;(i,j) / Vx.at&lt;float&gt;(i,j));
              }
  
              theta.at&lt;float&gt;(i,j) = calc;
  
              // Perform low-pass filtering
              float angle = 2 * calc;
              lowPassX.at&lt;float&gt;(i,j) = cos(angle * pi / 180);
              lowPassY.at&lt;float&gt;(i,j) = sin(angle * pi / 180);
  
              float sum3 = 0.0;
              float sum4 = 0.0;
  
              for(int u = -lowPassSize / 2; u &lt; lowPassSize / 2; u++)
              {
                 for(int v = -lowPassSize / 2; v &lt; lowPassSize / 2; v++)
                 {
                    sum3 += inputImage.at&lt;float&gt;(u,v) * lowPassX.at&lt;float&gt;(i - u*lowPassSize, j - v * lowPassSize);
                    sum4 += inputImage.at&lt;float&gt;(u, v) * lowPassY.at&lt;float&gt;(i - u*lowPassSize, j - v * lowPassSize);
                 }
              }
          lowPassX2.at&lt;float&gt;(i,j) = sum3;
          lowPassY2.at&lt;float&gt;(i,j) = sum4;
  
          float calc2 = 0.0;
  
          if(sum3 != 0 &amp;&amp; sum4 != 0)
          {
             calc2 = 0.5 * atan(lowPassY2.at&lt;float&gt;(i, j) / lowPassX2.at&lt;float&gt;(i, j)) * 180 / pi;
          }
          orientationMat.at&lt;float&gt;(i,j) = calc2;
  
          }
  
      }
  return orientationMat;
  
  }
  </code></pre>
  
  <p>I've already searched a lot on the web, but almost all of them are in Matlab. And there exist very few ones using OpenCV, but they didn't help me either. I sincerely hope someone could go through my code and point out any error to help. Thank you in advance. </p>
  
  <p><strong>Update</strong></p>
  
  <p>Here are the steps that I followed according to the paper:</p>
  
  <ol>
  <li>Obtain normalized image G.</li>
  <li>Divide G into blocks of size wxw (16x16).</li>
  <li>Compute the x and y gradients at each pixel (i,j).</li>
  <li><p>Estimate the local orientation of each block centered at pixel (i,j) using equations:
  <img src="https://i.stack.imgur.com/lKfbR.png" alt="enter image description here"></p></li>
  <li><p>Perform low-pass filtering to remove noise. For that, convert the orientation image into a continuous vector field defined as:</p></li>
  </ol>
  
  <p><img src="https://i.stack.imgur.com/5QR4O.png" alt="enter image description here"></p>
  
  <p><img src="https://i.stack.imgur.com/Y2ydY.png" alt="enter image description here"></p>
  
  <p>where W is a two-dimensional low-pass filter, and w(phi) x w(phi) is its size, which equals to 5. </p>
  
  <ol start="6">
  <li>Finally, compute the local ridge orientation at (i,j) using:</li>
  </ol>
  
  <p><img src="https://i.stack.imgur.com/Ol9I6.png" alt="enter image description here"></p>
  
  <p><strong>Update2</strong></p>
  
  <p>This is the output of orientationMat after changing the mat type to CV_16SC1 in Sobel operation as Micka suggested:</p>
  
  <p><img src="https://i.stack.imgur.com/A5dYx.png" alt="enter image description here"></p>