14 Replies  40022 Views  Last Post: 26 February 2009  11:03 PM
#1
How to calculate Hessian Matrix
Posted 29 May 2008  10:24 AM
Hessian is 2nd order differential of Gaussian. And has matrix 2x2.
Now, what number is in Hessian Matrix.
example:
Gaussian equation = F(x,y)= (1/(2*pi*T^2))exp((x^2+y^2)/2T^2)
Then,
H = [Fxx(x,y) Fxy(x,y) ;
Fyx(x,y) Fyy(x,y)]
So, what is significant number of matrix content? Then this matrix will be convoluted with image.
Please, If u know, answer that question, Thanks Before. [D'Oh!]
please, send me through:
[email protected]
Replies To: How to calculate Hessian Matrix
#2
Re: How to calculate Hessian Matrix
Posted 30 May 2008  04:36 AM
#3
Re: How to calculate Hessian Matrix
Posted 19 June 2008  08:12 AM
joske, on 30 May, 2008  04:36 AM, said:
Ya, Hessian Matrix is 2nd partial derivatives of a function. A function used is Gaussian Function.
After I calculated manually and considered x = y = 0 and sigma (T) = 1, I get the form of Hessian matrix as:
H = [Fxx Fxy ; Fyx Fyy] (mask 2x2)
H = [0.1592 0 ; 0 0.1592]
Then this matrix is convoluted with image. The result is not good, no like I hope. This matrix is convoluted with image of retinal vessel (grayscale), then the result has no uniform intensity. I hope output image has uniform intensity so that it is easily got threshold to get clear vessel segmentation.
The question is:
1. does the calculation has mistake?
2. Can U explain clearly about ur statement ".............and use the calculated derivatives to calculate the four elements of the Hessian matrix. An other solution is to calculate the second derivatives numerically for an arbitrary F(x,y)"
Thank U so much For Ur time. Glad to know You Mr. Joske
This post has been edited by jimblumberg: 30 March 2012  08:06 AM
#4
Re: How to calculate Hessian Matrix
Posted 19 June 2008  09:12 AM
Analytically: if you have to deal only with one known function, you can calculate the derivatives analytically and hardcode them in your program.
Numerically: It is very easy to calculate a derivative numerically: if you want to calculate the derivative of function f(x) at x, you just calculate f(xh) and f(x+h), where h is a very very small value. The derivative is approximately (f(x+h)  f(xh)) / (2*h). (You just calculate the slope of the function).
I entered your formula in Maple and calculated the derivatives, see attached screenshot.
F(x,y)= (1/(2*pi*T^2))*exp(((x^2+y^2)/(2*T^2)))
(Please ensure that this formula is correct and that I entered it correctly)
Edit: I evaluated the equations for x=0, y=0, T=1, as you did. I get the same results. (see screenshot)
Attached image(s)
This post has been edited by joske: 19 June 2008  09:24 AM
#5
Re: How to calculate Hessian Matrix
Posted 19 June 2008  04:39 PM
joske, on 19 Jun, 2008  09:12 AM, said:
Analytically: if you have to deal only with one known function, you can calculate the derivatives analytically and hardcode them in your program.
Numerically: It is very easy to calculate a derivative numerically: if you want to calculate the derivative of function f(x) at x, you just calculate f(xh) and f(x+h), where h is a very very small value. The derivative is approximately (f(x+h)  f(xh)) / (2*h). (You just calculate the slope of the function).
I entered your formula in Maple and calculated the derivatives, see attached screenshot.
F(x,y)= (1/(2*pi*T^2))*exp(((x^2+y^2)/(2*T^2)))
(Please ensure that this formula is correct and that I entered it correctly)
Edit: I evaluated the equations for x=0, y=0, T=1, as you did. I get the same results. (see screenshot)
Wow, You get the same matrix of mine.
Mr. Joske, would you please give me calculation again by using numerically method? "........ (f(x+h)  f(xh)) / (2*h). (You just calculate the slope of the function)........"
I dont undertand about that.
And I attach my visual c++ program and hessian matrix result
oke
Attached File(s)

program_result.doc (560K)
Number of downloads: 624 
program_c__.doc (127K)
Number of downloads: 987
This post has been edited by ardiansyahputra: 19 June 2008  04:54 PM
#9
Re: How to calculate Hessian Matrix
Posted 20 June 2008  03:07 AM
http://en.wikipedia....differentiation
Here is a small demo on how to calculate the first and second derivatives numerically with C++:
/* Numerical derivation [url="http://en.wikipedia.org/wiki/Numerical_differentiation"]http://en.wikipedia.org/wiki/Numerical_differentiation[/url] */ #include <cstdlib> #include <cstdio> #include <cmath> using namespace std; #define PAUSE {printf("Press Enter to continue...\n"); fflush(stdin); getchar(); fflush(stdin);} float pi = 3.141592654; double h = 0.00001; // define the function you want to use double F(double x, double y, double T) { return (1.0 / (2.0 * pi * T*T)) * exp(((x*x + y*y) / (2.0 * T*T))); } // first derivative to x double Fx(double x, double y, double T) { double fl = F(x  h, y, T) ; double fr = F(x + h, y, T) ; return (fr  fl) / (2.0 * h); } // first derivative to y double Fy(double x, double y, double T) { double fl = F(x, y  h, T) ; double fr = F(x, y + h, T) ; return (fr  fl) / (2.0 * h); } // second derivative of fx to x double Fxx(double x, double y, double T) { double fl = Fx(x  h, y, T) ; double fr = Fx(x + h, y, T) ; return (fr  fl) / (2.0 * h); } // second derivative of fx to y double Fxy(double x, double y, double T) { double fl = Fx(x, y  h, T) ; double fr = Fx(x, y + h, T) ; return (fr  fl) / (2.0 * h); } // second derivative of fy to x double Fyx(double x, double y, double T) { double fl = Fy(x  h, y, T) ; double fr = Fy(x + h, y, T) ; return (fr  fl) / (2.0 * h); } // second derivative of fy to x double Fyy(double x, double y, double T) { double fl = Fy(x, y  h, T) ; double fr = Fy(x, y + h, T) ; return (fr  fl) / (2.0 * h); } int main(/*int argc, char *argv[]*/) { double x = 0.0; double y = 0.0; double T = 1.0; printf("F(%g, %g, %g)=%g\n", x, y, T, F(x, y, T)); printf("\n"); printf("Fx(%g, %g, %g)=%g\n", x, y, T, Fx(x, y, T)); printf("Fy(%g, %g, %g)=%g\n", x, y, T, Fy(x, y, T)); printf("\n"); printf("Fxx(%g, %g, %g)=%g\n", x, y, T, Fxx(x, y, T)); printf("Fxy(%g, %g, %g)=%g\n", x, y, T, Fxy(x, y, T)); printf("Fyx(%g, %g, %g)=%g\n", x, y, T, Fyx(x, y, T)); printf("Fyy(%g, %g, %g)=%g\n", x, y, T, Fyy(x, y, T)); printf("\n"); PAUSE; return EXIT_SUCCESS; }
You see, using numerical derivation is very easy, and you don't have to do any difficult analytical derivations .
#10
Re: How to calculate Hessian Matrix
Posted 20 June 2008  05:25 PM
joske, on 20 Jun, 2008  03:07 AM, said:
http://en.wikipedia....differentiation
Here is a small demo on how to calculate the first and second derivatives numerically with C++:
/* Numerical derivation [url="http://en.wikipedia.org/wiki/Numerical_differentiation"]http://en.wikipedia.org/wiki/Numerical_differentiation[/url] */ #include <cstdlib> #include <cstdio> #include <cmath> using namespace std; #define PAUSE {printf("Press Enter to continue...\n"); fflush(stdin); getchar(); fflush(stdin);} float pi = 3.141592654; double h = 0.00001; // define the function you want to use double F(double x, double y, double T) { return (1.0 / (2.0 * pi * T*T)) * exp(((x*x + y*y) / (2.0 * T*T))); } // first derivative to x double Fx(double x, double y, double T) { double fl = F(x  h, y, T) ; double fr = F(x + h, y, T) ; return (fr  fl) / (2.0 * h); } // first derivative to y double Fy(double x, double y, double T) { double fl = F(x, y  h, T) ; double fr = F(x, y + h, T) ; return (fr  fl) / (2.0 * h); } // second derivative of fx to x double Fxx(double x, double y, double T) { double fl = Fx(x  h, y, T) ; double fr = Fx(x + h, y, T) ; return (fr  fl) / (2.0 * h); } // second derivative of fx to y double Fxy(double x, double y, double T) { double fl = Fx(x, y  h, T) ; double fr = Fx(x, y + h, T) ; return (fr  fl) / (2.0 * h); } // second derivative of fy to x double Fyx(double x, double y, double T) { double fl = Fy(x  h, y, T) ; double fr = Fy(x + h, y, T) ; return (fr  fl) / (2.0 * h); } // second derivative of fy to x double Fyy(double x, double y, double T) { double fl = Fy(x, y  h, T) ; double fr = Fy(x, y + h, T) ; return (fr  fl) / (2.0 * h); } int main(/*int argc, char *argv[]*/) { double x = 0.0; double y = 0.0; double T = 1.0; printf("F(%g, %g, %g)=%g\n", x, y, T, F(x, y, T)); printf("\n"); printf("Fx(%g, %g, %g)=%g\n", x, y, T, Fx(x, y, T)); printf("Fy(%g, %g, %g)=%g\n", x, y, T, Fy(x, y, T)); printf("\n"); printf("Fxx(%g, %g, %g)=%g\n", x, y, T, Fxx(x, y, T)); printf("Fxy(%g, %g, %g)=%g\n", x, y, T, Fxy(x, y, T)); printf("Fyx(%g, %g, %g)=%g\n", x, y, T, Fyx(x, y, T)); printf("Fyy(%g, %g, %g)=%g\n", x, y, T, Fyy(x, y, T)); printf("\n"); PAUSE; return EXIT_SUCCESS; }
You see, using numerical derivation is very easy, and you don't have to do any difficult analytical derivations .
It's Helpful information, Great!
Great..........
Mr. Joske, would U please add this program in convoluting between hessian matrix and original image (convolution). PLEASE, result of hessian matrix will be convoluted with original image. I want to know convolution result.
thx b4.
#11
Re: How to calculate Hessian Matrix
Posted 21 June 2008  12:56 AM
(see also http://en.wikipedia....lution_theorem)
As with numerial derivation, numerical integration is easy, have a look here: http://en.wikipedia....cal_integration.
An integral can be calculated with the rectangle rule.
/** integral to calculate int(K(x), x=x_start..x_end) */ double Kint(double x_start, double x_end) { // making the step smaller gives more accurate results (and requires more computation) double x_step = (x_end  x_start) / 100; double area = 0.0; // we start at x_start plus half the step width to calculate the function // at the middle of the elements for (double x = x_start + x_step/2; x < x_end; x += x_step) { area += K(x) * x_step; } return area; }
In your case you have to do an integral for a function with two variables
/** double integral to calculate int(K(x,y), x=x_start..x_end, y=y_start..y_end) */ double Kint(double x_start, double x_end, double y_start, double y_end) { // making the x_step smaller gives more accurate results (and requires more computation) double x_step = (x_end  x_start) / 100; double y_step = (y_end  y_start) / 100; double volume = 0.0; // we start at x_start and y_start plus half the step width to // calculate the function at the middle of the elements for (double x = x_start + x_step/2; x < x_end; x += x_step) { for (double y = y_start + y_step/2; y < y_end; y += y_step) { volume += K(x, y) * x_step * y_step; } } return volume; }
Now, to do the convolution between two functions, lets say your image Img(x,y) and the gaussian function F(x,y), you can do:
/** double convolution of Img(x,y) with F(x,y) around (x_center, y_center) */ double Kint(double x_center, double y_center) { // choose appropriate start and end values, // such that F(x_center  x, y_center  y) is about zero at the // boundaries (in this case F(100, 100) and F(100,100), etc.) double x_start = x_center  100; double x_end = x_center + 100; double y_start = y_center  100; double y_end = y_center + 100; // making the x_step smaller gives more accurate results (and requires more computation) double x_step = (x_end  x_start) / 100; double y_step = (y_end  y_start) / 100; // we start at x_start and y_start plus half the step width to // calculate the function at the middle of the elements double volume = 0.0; for (double x = x_start + x_step/2; x < x_end; x += x_step) { for (double y = y_start + y_step/2; y < y_end; y += y_step) { volume += (Img(x, y) * F(x_center  x, y_center  y)) * x_step * y_step; } } return volume; }
In fact, you take all image values around (x_center,y_center), and you multiply them with a value from F(), which is large around (x_center,y_center), and small when getting further away from the center. Finally, you sum the results up. So image values near (x_center,y_center) contribute much to the total, and image values further away contribute only little.
#12
Re: How to calculate Hessian Matrix
Posted 04 July 2008  10:14 PM
It's a long time not to make a communication with U.
I have a big question 4u One more. please answer myquestion, ya.
Myquestion have been attached in microsoft word. Please download it!
I am waiting for Ur reply. thx
MATRIKS_HESSIAN_DARI_FUNGSI_GAUSSIAN_oke.doc (51.5K)
Number of downloads: 368
#13
Re: How to calculate Hessian Matrix
Posted 05 July 2008  03:48 AM
ardiansyahputra said:
what is value of x and y?
I'm not sure if i understand what you mean. You can't "get" the values of x and y from the first derivative, you have to provide a value of x and y, and then the first derivative tells you what the slope of the function at that specific location is.
fx(x,y), the first derivative to x, gives you the slope of the function f(x,y) in the direction of the xaxis at the given position (x,y).
fy(x,y), the first derivative to y, gives you the slope of the function f(x,y) in the direction of the yaxis at the given position (x,y).
#14
Re: How to calculate Hessian Matrix
Posted 05 July 2008  08:38 PM
joske, on 5 Jul, 2008  03:48 AM, said:
ardiansyahputra said:
what is value of x and y?
I'm not sure if i understand what you mean. You can't "get" the values of x and y from the first derivative, you have to provide a value of x and y, and then the first derivative tells you what the slope of the function at that specific location is.
fx(x,y), the first derivative to x, gives you the slope of the function f(x,y) in the direction of the xaxis at the given position (x,y).
fy(x,y), the first derivative to y, gives you the slope of the function f(x,y) in the direction of the yaxis at the given position (x,y).
I hope U are not bored to receive myquestion.
Okey,
If,
fx(x,y) = 1
fy(x,y) = 0
How to get value of x and y ?
I am confused because there is exponent parameter.
Thx,
#15
Re: How to calculate Hessian Matrix
Posted 06 July 2008  03:53 AM
I have filled your equations in in Maple to solve the problem analytically, see snapshots. I think you definitely can use a program such as Maple or Mathematica for your project...
fy(x,y) = 0 has only a solution for y=0, as you can also see from the graph (where I used sigma=1).
fx(x,y) = 1 can then be calculated by filling in y=0, and filling in the value of sigma, and calculating x by the given function for x.
Attached image(s)
This post has been edited by joske: 06 July 2008  04:07 AM
#16
Re: How to calculate Hessian Matrix
Posted 06 July 2008  04:12 AM
sigma = 0.5; ' fill in your value y = 0; fx(x) = 2*x/(4*Pi*sigma^4)*Exp((x^2+y^2)/(2*sigma^2)); plot(fx(x), 1); 'we want to find where fx = 1 ' now, in the plot, click on the points where the two functions cross each other ' the tracepoint will snap to the intersection, then you can read the yvalue
This post has been edited by joske: 06 July 2008  04:12 AM
#17
Re: How to calculate Hessian Matrix
Posted 26 February 2009  10:58 PM
joske, on 20 Jun, 2008  11:56 PM, said:
(see also http://en.wikipedia....lution_theorem)
As with numerial derivation, numerical integration is easy, have a look here: http://en.wikipedia....cal_integration.
An integral can be calculated with the rectangle rule.
/** integral to calculate int(K(x), x=x_start..x_end) */ double Kint(double x_start, double x_end) { // making the step smaller gives more accurate results (and requires more computation) double x_step = (x_end  x_start) / 100; double area = 0.0; // we start at x_start plus half the step width to calculate the function // at the middle of the elements for (double x = x_start + x_step/2; x < x_end; x += x_step) { area += K(x) * x_step; } return area; }
In your case you have to do an integral for a function with two variables
/** double integral to calculate int(K(x,y), x=x_start..x_end, y=y_start..y_end) */ double Kint(double x_start, double x_end, double y_start, double y_end) { // making the x_step smaller gives more accurate results (and requires more computation) double x_step = (x_end  x_start) / 100; double y_step = (y_end  y_start) / 100; double volume = 0.0; // we start at x_start and y_start plus half the step width to // calculate the function at the middle of the elements for (double x = x_start + x_step/2; x < x_end; x += x_step) { for (double y = y_start + y_step/2; y < y_end; y += y_step) { volume += K(x, y) * x_step * y_step; } } return volume; }
Now, to do the convolution between two functions, lets say your image Img(x,y) and the gaussian function F(x,y), you can do:
/** double convolution of Img(x,y) with F(x,y) around (x_center, y_center) */ double Kint(double x_center, double y_center) { // choose appropriate start and end values, // such that F(x_center  x, y_center  y) is about zero at the // boundaries (in this case F(100, 100) and F(100,100), etc.) double x_start = x_center  100; double x_end = x_center + 100; double y_start = y_center  100; double y_end = y_center + 100; // making the x_step smaller gives more accurate results (and requires more computation) double x_step = (x_end  x_start) / 100; double y_step = (y_end  y_start) / 100; // we start at x_start and y_start plus half the step width to // calculate the function at the middle of the elements double volume = 0.0; for (double x = x_start + x_step/2; x < x_end; x += x_step) { for (double y = y_start + y_step/2; y < y_end; y += y_step) { volume += (Img(x, y) * F(x_center  x, y_center  y)) * x_step * y_step; } } return volume; }
In fact, you take all image values around (x_center,y_center), and you multiply them with a value from F(), which is large around (x_center,y_center), and small when getting further away from the center. Finally, you sum the results up. So image values near (x_center,y_center) contribute much to the total, and image values further away contribute only little.
joske, on 6 Jul, 2008  03:12 AM, said:
sigma = 0.5; ' fill in your value y = 0; fx(x) = 2*x/(4*Pi*sigma^4)*Exp((x^2+y^2)/(2*sigma^2)); plot(fx(x), 1); 'we want to find where fx = 1 ' now, in the plot, click on the points where the two functions cross each other ' the tracepoint will snap to the intersection, then you can read the yvalue
#18
Re: How to calculate Hessian Matrix
Posted 26 February 2009  11:03 PM
joske, on 6 Jul, 2008  03:12 AM, said:
sigma = 0.5; ' fill in your value y = 0; fx(x) = 2*x/(4*Pi*sigma^4)*Exp((x^2+y^2)/(2*sigma^2)); plot(fx(x), 1); 'we want to find where fx = 1 ' now, in the plot, click on the points where the two functions cross each other ' the tracepoint will snap to the intersection, then you can read the yvalue
how to find partial integral????