# How to calculate Hessian Matrix

Page 1 of 1

## 14 Replies - 32166 Views - Last Post: 26 February 2009 - 11:03 PMRate Topic: //<![CDATA[ rating = new ipb.rating( 'topic_rate_', { url: 'http://www.dreamincode.net/forums/index.php?app=forums&module=ajax&section=topics&do=rateTopic&t=53310&amp;s=ee6fb810a740da7e0dc199f97fbde1d3&md5check=' + ipb.vars['secure_hash'], cur_rating: 0, rated: 0, allow_rate: 0, multi_rate: 1, show_rate_text: true } ); //]]>

### #1 ardiansyahputra

Reputation: 2
• Posts: 8
• Joined: 29-May 08

# How to calculate Hessian Matrix

Posted 29 May 2008 - 10:24 AM

Hai, I am confused about Hessian.
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.
poetra_ardiansyah@yahoo.co.id

Is This A Good Question/Topic? 0

## Replies To: How to calculate Hessian Matrix

### #2 joske

• D.I.C Regular

Reputation: 43
• Posts: 297
• Joined: 04-September 07

## Re: How to calculate Hessian Matrix

Posted 30 May 2008 - 04:36 AM

The hessian matrix contains the second partial derivatives of the function F(x,y). If you want to calculate this for the particular function that you give, you can calculate the the second derivatives analytically, 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). What do you need?

### #3 ardiansyahputra

Reputation: 2
• Posts: 8
• Joined: 29-May 08

## Re: How to calculate Hessian Matrix

Posted 19 June 2008 - 08:12 AM

joske, on 30 May, 2008 - 04:36 AM, said:

The hessian matrix contains the second partial derivatives of the function F(x,y). If you want to calculate this for the particular function that you give, you can calculate the the second derivatives analytically, 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). What do you need?

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 joske

• D.I.C Regular

Reputation: 43
• Posts: 297
• Joined: 04-September 07

## Re: How to calculate Hessian Matrix

Posted 19 June 2008 - 09:12 AM

Well, you can solve your problem in two ways:
Analytically: if you have to deal only with one known function, you can calculate the derivatives analytically and hard-code 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(x-h) and f(x+h), where h is a very very small value. The derivative is approximately (f(x+h) - f(x-h)) / (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 ardiansyahputra

Reputation: 2
• Posts: 8
• Joined: 29-May 08

## Re: How to calculate Hessian Matrix

Posted 19 June 2008 - 04:39 PM

joske, on 19 Jun, 2008 - 09:12 AM, said:

Well, you can solve your problem in two ways:
Analytically: if you have to deal only with one known function, you can calculate the derivatives analytically and hard-code 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(x-h) and f(x+h), where h is a very very small value. The derivative is approximately (f(x+h) - f(x-h)) / (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(x-h)) / (2*h). (You just calculate the slope of the function)........"

And I attach my visual c++ program and hessian matrix result

oke

#### Attached File(s)

This post has been edited by ardiansyahputra: 19 June 2008 - 04:54 PM

### #9 joske

• D.I.C Regular

Reputation: 43
• Posts: 297
• Joined: 04-September 07

## 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 ardiansyahputra

Reputation: 2
• Posts: 8
• Joined: 29-May 08

## 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 .

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 joske

• D.I.C Regular

Reputation: 43
• Posts: 297
• Joined: 04-September 07

## Re: How to calculate Hessian Matrix

Posted 21 June 2008 - 12:56 AM

The convolution you want to do is in fact an integral of your gaussian fuction, and the image values.

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 ardiansyahputra

Reputation: 2
• Posts: 8
• Joined: 29-May 08

## Re: How to calculate Hessian Matrix

Posted 04 July 2008 - 10:14 PM

hI, Mr. Joske.
It's a long time not to make a communication with U.
I have a big question 4u One more. please answer myquestion, ya.
I am waiting for Ur reply. thx
MATRIKS_HESSIAN_DARI_FUNGSI_GAUSSIAN_oke.doc (51.5K)

### #13 joske

• D.I.C Regular

Reputation: 43
• Posts: 297
• Joined: 04-September 07

## Re: How to calculate Hessian Matrix

Posted 05 July 2008 - 03:48 AM

ardiansyahputra said:

how to get value of x and y from first order derivative of x and first order derivative of y ?

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 x-axis 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 y-axis at the given position (x,y).

### #14 ardiansyahputra

Reputation: 2
• Posts: 8
• Joined: 29-May 08

## Re: How to calculate Hessian Matrix

Posted 05 July 2008 - 08:38 PM

joske, on 5 Jul, 2008 - 03:48 AM, said:

ardiansyahputra said:

how to get value of x and y from first order derivative of x and first order derivative of y ?

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 x-axis 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 y-axis 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 joske

• D.I.C Regular

Reputation: 43
• Posts: 297
• Joined: 04-September 07

## 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 joske

• D.I.C Regular

Reputation: 43
• Posts: 297
• Joined: 04-September 07

## Re: How to calculate Hessian Matrix

Posted 06 July 2008 - 04:12 AM

So, to calculate the actual solution (without having Maple), you can for example use the free program SpeQ Mathematics and enter:
```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 y-value

```

This post has been edited by joske: 06 July 2008 - 04:12 AM

### #17 VRVNTH

Reputation: 0
• Posts: 2
• Joined: 26-February 09

## Re: How to calculate Hessian Matrix

Posted 26 February 2009 - 10:58 PM

joske, on 20 Jun, 2008 - 11:56 PM, said:

The convolution you want to do is in fact an integral of your gaussian fuction, and the image values.

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:

So, to calculate the actual solution (without having Maple), you can for example use the free program SpeQ Mathematics and enter:
```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 y-value

```

### #18 VRVNTH

Reputation: 0
• Posts: 2
• Joined: 26-February 09

## Re: How to calculate Hessian Matrix

Posted 26 February 2009 - 11:03 PM

joske, on 6 Jul, 2008 - 03:12 AM, said:

So, to calculate the actual solution (without having Maple), you can for example use the free program SpeQ Mathematics and enter:
```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 y-value

```

how to find partial integral????