block processing to find local orientation of each block

I want to implement an algorithm in MATLAB on a Fingerprint image in order to enhance it. The algorithm starts with normalization then orientation image estimation followed by three more steps. I have completed the normalization part but stuck at the orientation step for which following are the steps: 1) Divide the image into 16x16 blocks. 2) Compute the gradient at each pixel. 3) Estimate the local orientation of each block centered at pixel(i,j) using the below

Here delX and delY are gradients at each pixel. Vx(i,j) and Vy(i,j) are orientations.

    function [imgOrient] = orientation(imgNorm)
    Grad =@(block_struct)gradient(block_struct.data)
    imgBlock=blockproc(ingNorm, [16 16], 'Grad)

Please let me know if this is correct way to use blockproc. Now I want to ask that how to save both the both the components (x as well as y ) returned by gradient function and use them in above formula. Any help would be appreciated.

[Leon]

Hello, I would like to know if you are using the original image or the normalized image for estimating the image direction field

Answer 1

[Guillaume]

Your code at present is just equivalent to

imgBlock = gradient(imgNorm);

Not really what you want.

What you want to do is pass to blockproc a function that returns a single value for each block, the orientation. You can't do it with an anonymous function (well, you could but that would be very unwieldy), so use a normal function:

function orientation = CalculateBlockOrientation(block_struct)
   %block_struct: block structure generated by blockproc
   %orientation: the orientation of the block
   [delx, dely] = gradient(block_struct.data);
   Vx = sum(sum(2*delx.*dely));
   Vy = ... %I'm sure you can work it out
   orientation = 0.5*atan(Vy/Vx);
end

And in your main function

blockorientation = blockproc(imgNorm, [16, 16], @CalculateBlockOrientation);

[Gunjan]

Thank you so much Guillaume. I have a couple of doubts in you code.

• What is 'block'? Is it same as 'block_struct.data' which I've used in my code?
• Let there are total n block then how it would keep on taking block after block.
• I also want to print these orientation (in the form of directional arrows ) in the image. Do I need to use quiver ?

[Guillaume]

Yes, sorry, completely forgot that blockproc passes a structure to the processing function. I've amended my answer to correct the mistake.

I'm not sure I understand your second question. It's blockproc's job to pass block after block to the processing function.

Yes, quiver will let you print arrows on your image. If you want the arrows in the middle of the block, you'd use 8:16:size(img, dim) for x and y, and calculate u and v using the angle returned by each block:

arrow_length = 10;
u = cos(blockorientation) * arrow_length;
v = sin(blockorientation) * arrow_length;
quiver(8:16:size(imgNorm, 2), 8:16:size(imgNorm, 1), u, v);

[Gunjan]

Thank you so much Guillaume. You helped a lot. I am new to MATLAB so facing these problems. I try to implement this. Again thanks :)

[Gunjan]

Hi Guillaume, Why did you use size function in quiver ? Is is for mentioning the coordinates over the entire image ? My code is giving different errors including one related to strfind function but I did not use it at all. What to do?

[Guillaume]

8:16:size(imgNorm, x)

will generate coordinates for the quiver from 8 to the size of the image in dimension x, in steps of 16 pixels.

There is no string manipulation involved in this so if you get an error with strfind it's nothing to do with it.

[Gunjan]

Hi Guillaume, Need your help. In order to improve the orientation I need to apply a low pass filter (say mean filter) over each of the block and based on the attached formula i need to compute the finer orientation. Here my problem is, I do not understand how to apply low pass filter for each block based on the below formula in which W is the low pass filter. How to pad a block and calculate the mean. Default size of the filter is 5x5.

[Gunjan]

function new_orientation = CalculateBlockOrientation(block_struct)

    [delx, dely] = gradient(block_struct.data);
     Vx = sum(sum(2*delx.*dely));
     Vy = sum(sum(delx^2.*dely^2));
     orientation = 0.5*atan(Vy/Vx);

     fi_x=cos(2*orientation);
     fi_y=sin(2*orientation);

     mask=(1/25)*ones(5,5);

     new_fi_x=conv2(mask,fi_x);
     new_fi_y=conv2(mask,fi_y);

     new_orientation = 0.5*atan( new_fi_y/ new_fi_x);

end

[Guillaume]

You seem to have answered your own question. Your low-pass filter looks indeed like a convolution. What's not clear from your post is what Φx and Φy are. These need to be 2D matrices for each block, not scalars. (But IIUC, Φ' is a scalar for a block)

Possibly, this is Φ:

[delx, dely] = gradient(block_struct.data);
phi_x = 2 * delx  .* dely;
phi_y = delx .^2 + dely .^2;

Therefore the convolution and calculation of Φ' would be:

w = ones(5, 5) / 25;
phi_xprime = sum(sum(conv2(w, phi_x, 'same')));
phi_yprime = sum(sum(conv2(w, phi_y, 'same')));
orientation = 0.5 * atan(phi_yprime / phi_xprime);

[Gunjan]

In order to perform the low-pass filtering, the orientation image needs to be converted into a continuous vector field, which is defined as follows: Φx(i, j) = cos(2*theta(i, j)) and Φy(i, j) = sin(2*theta(i, j)) where Φx and Φy are the x and y components of the vector field, respectively. Here, theta is 'Orientation' calculated in the above posts. Could you please tell me why to use sum while convolution function already handles this ?

[Guillaume]

I'm afraid I don't understand anymore what you are attempting to do. Within a block, θ is scalar, and hence your Φx and Φy are scalar. There is nothing to average. You could of course, filter the whole θ image but that would be done after the call to blockproc, so:

blockorientation = blockproc(imgNorm, [16, 16], @CalculateBlockOrientation); %CalculateBlockOrientation is unchanged
phi_x = cos(2*blockorientation);
phi_y = sin(2*blockorientation);
w = ones(5) / 25;
phi_x_prime = conv2(phi_x, w, 'same');
phi_y_prime = conv2(phi_y, w, 'same');

My answer was based on my understanding that the filter was to be done inside the blockproc and assumed that Φx and Φy was the thing that was summed in your original question to obtain νx and νy. The |sum(sum(| is used to convert the result of the convolution (a matrix) into a scalar, as you did originally for ν. That is what I understood:

original:

ν = Σ Σ Φ

new:

Φ' = Σ Σ (W ⊗ Φ) %convolution
ν' = Σ Σ Φ'

therefore:

&nu' = Σ Σ conv2(W, &Phi, 'same');

[Gunjan]

<<

>>

I think I'm in a big confusion. If you find time could you please look at the attached image with this post and let me know where thing are going wrong ? I do not know how to help me.

[Guillaume]

Right, so the filtering is done on the orientation and is done outside of blockproc. This corresponds to the 1st part of my answer:

blockorientation = blockproc(imgNorm, [16, 16], @CalculateBlockOrientation); %CalculateBlockOrientation is unchanged
phi_x = cos(2*blockorientation);
phi_y = sin(2*blockorientation);
w = ones(5) / 25;
phi_x_prime = conv2(phi_x, w, 'same');
phi_y_prime = conv2(phi_y, w, 'same');

[ALI Salah]

hi Guillaume the function CalculateBlockOrientation will return the oriantation of each block in normimg but the question is how can i abtain new oriantation for each block based on

thetax(i,j) = cos(2*blockorientation);

thetay(i,j) = sin(2*blockorientation);

and theta=0.5*tan((thetax(i,j))'/(thetay(i,j))');

the last 3 equation below

note that in the formula above of Gunjan on 27 Jan 2016

Answer 2

[Image Analyst]

Each block will have a variety of gradients in it. The Computer Vision System Toolbox has a HOG detector (Histogram Of Gradients), where you can find the dominant direction and other info about the texture. See http://www.mathworks.com/help/vision/examples/digit-classification-using-hog-features.html?prodcode=VP&language=en

[Gunjan]

Thank you for this information. Which type of function this gradient function do and how many neighborhoods it considers? What if I want to apply Marr-Hildreth or Sobel ? As the gradient function is slow so I want to get some knowledge about customizing this function.

[Bekhtaoui Abdelhakim]

hello mr Gunjan i have a probleme with the orientation estimation code I wonder if you could send your final code in my email [e-mail]

thank you