Weird orientation map using sobel. What's wrong?

Question

I tried to follow this method of drawing orientation map http://answers.opencv.org/question/9493/fingerprint-orientation-map-through-gradient/

And i used a block size of 5x5 on my 480x320 image. The gradients i got was from 0-270 degrees. And there are constant values that keep on repeating like 44.7623 and 224.762. I wonder if my gradients are wrong.

After that, i add all the gradients in the 5x5 block and divide them by 25 (averaging) like what the link said. I divided the degrees into 8 sections of 45degree intervals and plotted them out. But it looks nothing like my original image. Can anyone tell me what's wrong? I just want to detect the core(circle-like) feature of the image.

My original image is this _________________________________________________ But my orientation map is this:

This is what I'm doing

What's wrong ? =(

I got the gradients from this method :

/// Gradient X
 cv::Sobel(original_Mat, grad_x, CV_32FC1, 1, 0, 3);

 /// Gradient Y
 cv::Sobel(original_Mat, grad_y, CV_32FC1, 0, 1, 3);

  Mat orientation = Mat(grad_x.rows, grad_y.cols, CV_32F);

 for(int i = 0; i < grad_x.rows; i++){
    for(int j = 0; j < grad_x.cols; j++){
        // Retrieve a single value
        float valueX = grad_x.at<float>(i,j);
        float valueY = grad_x.at<float>(i,j);
        // Calculate the corresponding single direction, done by applying the arctangens function
        float result = fastAtan2(valueX,valueY);
        // Store in orientation matrix element
        orientation.at<float>(i,j) = result;

    }
    }

Here's the full code.

 int main()
 {

cv::Mat original_Mat=cv::imread("Source.bmp", 1);

cv::Mat grad = cv::Mat::zeros(original_Mat.size(),CV_64F);

 /// Generate grad_x and grad_y
 cv::Mat grad_x = cv::Mat::zeros(original_Mat.size(), CV_64F); 
 cv::Mat grad_y = cv::Mat::zeros(original_Mat.size(), CV_64F);
 cv::Mat grad_angle = cv::Mat::zeros(original_Mat.size(), CV_64F);

 /// Gradient X
 cv::Sobel(original_Mat, grad_x, CV_32FC1, 1, 0, 3);

 /// Gradient Y
 cv::Sobel(original_Mat, grad_y, CV_32FC1, 0, 1, 3);

  Mat orientation = Mat(grad_x.rows, grad_y.cols, CV_32F); //to store the gradients
  Mat img=Mat(grad_x.rows, grad_y.cols, CV_32F);//to draw out the map
  img = cv::Scalar(255,255,255);//all white

// Calculate orientations of gradients --> in degrees
// Loop over all matrix values and calculate the accompanied orientation

 for(int i = 0; i < grad_x.rows; i++){
    for(int j = 0; j < grad_x.cols; j++){
        // Retrieve a single value
        float valueX = grad_x.at<float>(i,j);
        float valueY = grad_x.at<float>(i,j);
        // Calculate the corresponding single direction, done by applying the arctangens function
        float result = fastAtan2(valueX,valueY);
        // Store in orientation matrix element
        orientation.at<float>(i,j) = result;

    }
    }


int i=0,j=0;
int x1=0,x2=0;
float results;

for(int l=0;l<96;l++) //to loop all the rows
{
    int x1=(5+(l*5)); // to get 5x5 block sizes

for(int k=0;k<64;k++)//to loop all the columns
{

    int x2=(5+(k*5)); // to get 5x5 block sizes
     results=0;


           //to get the total of 5x5 gradient values
    for(i=(x1-5); i < x1; i++){
    for(j=(x2-5); j < x2; j++){

        results=results+orientation.at<float>(i,j);
        orientation.at<float>(i,j)=0;


    }
    }

    results=results/25; //averaging the 5x5 block gradients

    orientation.at<float>((x1-3),(x2-3))=results; //to store the results in the center of the 5x5 block

}
}
results=0;


//this loop is to draw out the orientation map
for(int i=0;i<480;i++)
{
    for(int j=0;j<320;j++)
    {
        results=orientation.at<float>(i,j);
    if ((results<=22.5)&&(results>0)){
        results=0;
        img.at<int>(i,j)=255;
        img.at<int>(i,j+1)=255;
        img.at<int>(i,j+2)=255;

    }
    else if((results>22.5)&&(results<=67.5)){
        results=45;
        img.at<int>(i,j)=255;
        img.at<int>(i-1,j+1)=255;
        img.at<int>(i-2,j+2)=255;


    }
    else if((results>67.5)&&(results<=112.5)){
        results=90;
        img.at<int>(i,j)=255;
        img.at<int>(i-1,j)=255;
        img.at<int>(i-2,j)=255;

    }
    else if((results>112.5)&&(results<=157.5)){
        results=135;
        img.at<int>(i,j)=255;
        img.at<int>(i-1,j-1)=255;
        img.at<int>(i-2,j-2)=255;

    }
    else if((results>157.5)&&(results<=202.5)){
        results=180;
        img.at<int>(i,j)=255;
        img.at<int>(i,j-1)=255;
        img.at<int>(i,j-2)=255;

    }
    else if((results>202.5)&&(results<=247.5)){
        results=225;
        img.at<int>(i,j)=255;
        img.at<int>(i+1,j-1)=255;
        img.at<int>(i+2,j-2)=255;
        endx=x2-5;
        endy=x1-1;
    }
    else if((results>247.5)&&(results<=292.5)){
        results=270;
        img.at<int>(i,j)=255;
        img.at<int>(i+1,j)=255;
        img.at<int>(i+2,j)=255;

    }
    else if((results>292.5)&&(results<=337.5)){
        results=315;
        img.at<int>(i,j)=255;
        img.at<int>(i+1,j+1)=255;
        img.at<int>(i+2,j+2)=255;

    }
    else
    {
        results=0;

    }
    orientation.at<float>(i,j)=results;


    }
}

Answer 1

Here is my result: For image:

I've got result:

The code:

#include <stdio.h>
#include <stdarg.h>
#include "opencv2/opencv.hpp"
using namespace std;
using namespace cv;
int main(int argc, char* argv[])
{
    namedWindow("source");
    namedWindow("result");

    namedWindow("ang");

    Mat img=imread("D:\\ImagesForTest\\binarized_image.png",0);
    cv::threshold(img,img,128,255,cv::THRESH_BINARY);
    Mat thinned;

    thinned=img.clone(); // Just clone the input
    //Thinning(img,thinned); // Not actually needed

    cv::GaussianBlur(thinned,thinned,Size(3,3),1.0);
    Mat gx,gy,ang,mag;
    cv::Sobel(thinned,gx,CV_32FC1,1,0);
    cv::Sobel(thinned,gy,CV_32FC1,0,1);
    cv::phase(gx,gy,ang,false);
    cv::magnitude(gx,gy,mag);

    cv::normalize(mag,mag,0,1,cv::NORM_MINMAX);


    Mat angRes=Mat::zeros(img.rows*3,img.cols*3,CV_8UC1);

    for (int i=0;i< img.rows;i+=2)
    {
        for (int j=0;j< img.cols;j+=2)
        {
            int x=j*3;
            int y=i*3;

            float r=5;
            float m=r*(mag.at<float>(i,j));
            float dx=m*r*cos(ang.at<float>(i,j));
            float dy=m*r*sin(ang.at<float>(i,j));

            cv::line(angRes,cv::Point(x,y),cv::Point(x+dx,y+dy),Scalar::all(255),1,CV_AA);
        }
    }
    imshow("ang",angRes);
    imshow("source",img);
    imshow("result",thinned);
    cv::waitKey(0);
}

Enother variant (weighted block averages):

#include <stdio.h>
#include <stdarg.h>
#include "opencv2/opencv.hpp"
using namespace std;
using namespace cv;
float GetWeightedAngle(Mat& mag,Mat& ang)
{
    float res=0;
    float n=0;
    for (int i=0;i< mag.rows;++i)
    {
        for (int j=0;j< mag.cols;++j)
        {
            res+=ang.at<float>(i,j)*mag.at<float>(i,j);
            n+=mag.at<float>(i,j);
        }
    }
    res/=n;
    return res;
}


int main(int argc, char* argv[])
{
    namedWindow("source");
    namedWindow("ang");

    Mat img=imread("D:\\ImagesForTest\\binarized_image.png",0);
    cv::threshold(img,img,128,255,cv::THRESH_BINARY);
    Mat thinned;

    thinned=img.clone();
    //Thinning(img,thinned);

    //cv::GaussianBlur(thinned,thinned,Size(3,3),1.0);
    Mat gx,gy,ang,mag;
    cv::Sobel(thinned,gx,CV_32FC1,1,0,7);
    cv::Sobel(thinned,gy,CV_32FC1,0,1,7);
    cv::phase(gx,gy,ang,false);
    cv::magnitude(gx,gy,mag);

    cv::normalize(mag,mag,0,1,cv::NORM_MINMAX);


    Mat angRes=Mat::zeros(img.rows,img.cols,CV_8UC1);

    int blockSize=img.cols/15-1;
    float r=blockSize;

    for (int i=0;i< img.rows-blockSize;i+= blockSize)
    {
        for (int j=0;j< img.cols-blockSize;j+= blockSize)
        {

            float a=GetWeightedAngle(mag(Rect(j,i,blockSize,blockSize)),ang(Rect(j,i,blockSize,blockSize)));

            float dx=r*cos(a);
            float dy=r*sin(a);
            int x=j;
            int y=i;            

            cv::line(angRes,cv::Point(x,y),cv::Point(x+dx,y+dy),Scalar::all(255),1,CV_AA);
        }
    }
    imshow("ang",angRes);
    imshow("source",img);
    cv::waitKey(0);
}

It gives the result image: