# Compute a hessian matrix on a image?

Is it possible for Mathmatica to load a gray image and apply a hessian matrix on it?... It used image processing to find good features.

• Is "Apply a Hessian" the same as "Calculate the second derivatives matrices"? May 15, 2015 at 20:17
• yes @belisarius
– SURF
May 15, 2015 at 22:46
• It looks like you might need ImageFilter, but there are also many built in filters that may do what you want. May 16, 2015 at 6:10

I normally use something like this:

img = ExampleData[{"TestImage", "Airport"}];

{h[xx], h[xy], h[yy]} =
GaussianFilter[ImageData[img], 5, #] & /@ {{0, 2}, {1, 1}, {2, 0}};


(Note: This line can be written much more tersely. But this way it should be clear what's going on: We're applying a Gaussian derivative filter, to get the 2nd order derivatives of the image, and assign the 3 results to the 3 variables h[xx], h[xy] and h[yy].

Now, h[xx], h[xy], h[yy] contain the 3 independent components of the Hessian at each pixel.

Then you can do symbolic calculations using simliar symbols, e.g. to calculate the eigenvalues of a generic symmetric 2x2 matrix:

eigenvalues = FullSimplify@Eigenvalues[{{m[xx], m[xy]}, {m[xy], m[yy]}}]


$\left\{\frac{1}{2} \left(-\sqrt{(m(\text{xx})-m(\text{yy}))^2+4 m(\text{xy})^2}+m(\text{xx})+m(\text{yy})\right),\frac{1}{2} \left(\sqrt{(m(\text{xx})-m(\text{yy}))^2+4 m(\text{xy})^2}+m(\text{xx})+m(\text{yy})\right)\right\}$

And then simply replace m with h in the symbolic result to apply it to the image Hessian:

imgEigenvalues = eigenvalues /. m -> h;


Now imgEigenvalues contains a 2d array for each of the two eigenvalues of the Hessian at every pixel:

GraphicsRow[Image /@ Rescale[imgEigenvalues]]


You can use that to search for image features, e.g. "where is the second eigenvalue larger than some value":

HighlightImage[img, Binarize[Image[imgEigenvalues[[2]]], 0.01]]


• Can you not just use DerivativeFilter or GradientFilter ? I believe there's also an example in the documentation: ridgeFilter[img_, s_: 1] := Module[{data = ImageData@img, Lxx,Lxy, Lyy},{Lxx, Lxy, Lyy} = DerivativeFilter[data, {{0, 2}, {1, 1}, {2,0}},s];Image[Chop[s^(3/2)/2 (Sqrt[(Lxx - Lyy)^2 + 4 Lxy^2] - Lxx - Lyy)]]] May 16, 2015 at 14:55
• @Histograms: As far as I understand the documentation, GradientFilter basically uses GaussianFilter, and it can only calculate the 1st derivative. You could use DerivativeFilter, of course. I generally start with GaussianFilter, as it's in some sense the most "natural" kernel and has (at least theoretically) certain unique advantages, especially if you calculate multiple scales. (The reasoning is too long for a comment, and honestly, I'd have to look most of it up myself, but I can give references if you're interested.) May 16, 2015 at 17:47
• Would it be possible mark those which is above a certain threshold for both eigenvalues?..
– SURF
May 16, 2015 at 18:16
• @SURF: Sure, you could use e.g. minEigenvalue = MapThread[Min, imgEigenvalues, 2]; to calculate the per-pixel min of the two eigenvalues, then binarize the result that. May 16, 2015 at 18:53
• Could these values also be calculated for different scales?... such as the input image still is the same, but the hessian is calculated on different scales..
– SURF
May 16, 2015 at 20:35