How do I calculate gradient and hessian matrix by two operators?

For each point (x,y), the gradients should be calculated by two operators like:

operator x:1/4*[1,0,-1;0,0,0;1,0,-1] operator y:1/4*[-1,0,-1;0,0,0;1,0,1]

now should I first measure the gradients by its function and then convolve them with these matrices as below? (I have no image processing toolbox)

Then hessian operators are [-1,1] for x , and [-1;1] for y. Well, I know how to generate hessian matrix but don't know how to do with these operators in a way that finally I can calculate the sum of all these measurements (their dimension would be different after convolution)

Answer 1

[David Young]

On the first part of your question, I think you should not call gradient, because convolution with the operators is what calculates the gradients. This makes more sense:

opx = [1,0,-1;0,0,0;1,0,-1]/4;
opy = [-1,0,-1;0,0,0;1,0,1]/4;
Gx = conv2(im, opx);
Gy = conv2(im, opy);

However, the choice of operators puzzles me. [1 0 -1; 2 0 -1; 1 0 -1] would be a more conventional choice for an x-gradient operator.

To get the elements of the Hessian, you can convolve Gx and Gy with gradient operators, something like

Gxx = conv2(Gx, opx);

and so on. Again, the motivation for the choice of operator isn't clear - the operators [-1 1] and [-1; 1] introduce a half-pixel offset into the results, so can be an awkward choice.

An alternative approach is just to use the gradient function instead of convolution. However, convolution gives you more control, and permits the introduction of some smoothing and the use of symmetrical operators.

Depending on your data, you may well need some smoothing as well as differencing to form estimates of gradients or Hessian elements.

[Klara]

Thanks for your answer, what about Gxy and Gyx? Can I say like

Gxy=conv2(opx,opy,Gx)
Gyx=conv2(opy,opx,Gy)

[David Young]

Yes, that's correct.

[Klara]

I applied 2D gaussian filter to smooth the gradients. but when I convolve them with gradient operators, for Gxy and Gyx an error appears in my conv2 command. I tried to calculate the hessian elements by convolving the Gx and Gy with hessian operators, then I smoothed again and they are calculated very fine without any difference in dimension of the vectors. But I tested them with a reference image result, and I found it's not what I want. Have no idea!

[David Young]

Not sure what the problem is, but you could try using my smoothing+differencing functions at this location on the File Exchange.

[Edit] Just remembered you don't have the IPT. I think the only thing you need to replace to use my FEX functions is fspecial, to get the Gaussian kernel. It's not complex - let me know if you need help with that.

[Klara]

Good idea, but since each of gradients and hessian elements are 2D array, I smoothed them seperately with gsmooth1 function(one along dim1, one along dim2). I didn't get same size after smoothing again. could you help please? [consider I can use fspecial].

[David Young]

The FEX submission I linked to above does, I think, what you want.

[Klara]

It was wrong to use conv2 instead of filter2. then there is no different dimension problem.

[David Young]

filter2 just calls conv2. Both have a 'same' option which avoids a change in size; the difference is that it is the default for filter2. In both cases the array is padded with zeros to compute values near the boundaries, which may or may not be what you want.