Frangi/Hessian for Anisotropic Voxel

[MOAMSA]

Hello,
I appreciate it if anyone explains how the Frangi algorithm handles anisotropic voxel?
Based on the Gradient definition (first/second-order), we have to include voxel size in different directions, but I couldn't follow this in hessian_matrix function. This function computes the Hessian matrix by convolving the image with the second derivatives of the Gaussian kernel.

Thanks

[jni]

@MOAMSA Hi! Thanks for the Q.

I must say that if we get something to work it'll be luck, based on the documentation. 😂 ie the Frangi filter does not currently support anisotropic voxels, at least not intentionally: the sigma is supposed to be a single float. However, if I follow the sequence of calls, I don't see a reason why you couldn't pass in anisotropic sigmas that ultimately get passed to ndi.gaussian_filter here:

scikit-image/skimage/feature/corner.py

Lines 187 to 188 in ed889ba

	gaussian_filtered = ndi.gaussian_filter(image, sigma=sigma,
	mode=mode, cval=cval)

In short, your sigmas should be inversely proportional to the voxel spacing. I just tried it and it actually works!

import numpy as np
from skimage import data, filters
import napari


cells = data.cells3d()
nuclei = cells[:, 1]
sigmas = np.arange(1, 10, 2)[:, np.newaxis] / [0.29, 0.26, 0.26]
frangied = filters.frangi(nuclei, sigmas)

with napari.gui_qt():
    viewer = napari.Viewer(ndisplay=3)
    viewer.add_image(
            frangied,
            colormap='magenta',
            blending='additive',
            rendering='attenuated_mip',
            contrast_limits=[0, 5e-6],
    )
    viewer.add_image(
            nuclei,
            colormap='green',
            blending='additive',
            rendering='attenuated_mip',
    )

Note that you need the skimage master version for the above to work.

I would certainly be happy with some documentation fixes, and maybe a gallery example, to make it clear that this works!

[MOAMSA]

@jni Perfect! As an alternative, do you think we can update the Hessian matrix by dividing each element by the corresponding voxel size? for example, D2F/Dx2 is divided by (voxel_x)^2...

[jni]

Ah, nice idea! Possibly we also need to do this? I need to think about it. But I don't think it's sufficient: the sigmas should definitely be scaled. Imagine a super anisotropic example, where vz=100 and vx=vy=1. You would not want to do a Gaussian blur with a sigma in np.arange(1, 10, 2) along the z dimension, right?

[MOAMSA]

Could you please take a look at this. It might be helpful.

On Napari, you haven't used scale parameter. How Napari consider the voxel size? just to be sure that you have got the proper 3D dimension...

[jni]

Could you please take a look at this. It might be helpful.

It is!

On Napari, you haven't used scale parameter. How Napari consider the voxel size? just to be sure that you have got the proper 3D dimension...

You're right. I should have added scale=[0.29, 0.26, 0.26] to the add_image calls.