[ITK] Obtaining Hessian matrix values

[Francois Budin]

Hello Nicolas,

Those are very good questions.
itkImageSSRTD33 is not currently wrapped in PyBuffer. I just submitted a
patch to take care of this.
However, you can use your other approach to get the pixel values, waiting
for the patch to be merged into ITK.

The order of the different components can be found by comparing the two
ways to access the values in a Symmetric Second Rank Tensor:
Following your code, you can add the following lines:
image=hessian_filter.GetOutput()
t=image.GetPixel([0,0,0])
print("%f,%f,%f,%f,%f,%f"%(t[0],t[1],t[2],t[3],t[4],t[5]))
print("%f,%f,%f,%f,%f,%f"%(t(0,0),t(0,1),t(0,2),t(1,1),t(1,2),t(2,2)))

You can see that the components are ordered the following way:
0 1 2
1 3 4
2 4 5

To extract the boundaries of a region:
size=hessian_filter.GetOutput().GetBufferedRegion().GetSize()
print("%d,%d,%d"%(size[0],size[1],size[2]))

To extract the eigenvectors and eigenvalues, you can do the following:
v=itk.FixedArray[itk.D,3]()
e=itk.Matrix[itk.D,3,3]()
t.ComputeEigenAnalysis(v,e)
print("%f,%f,%f"%(v[0],v[1],v[2]))
arr=itk.GetArrayFromVnlMatrix(e.GetVnlMatrix().as_matrix())
print arr[0,0]

Hope this helps,
Francois

On Fri, Mar 24, 2017 at 4:52 AM, Nicolas Cedilnik <nicolas.cedilnik at inria.fr
> wrote:

> Hello,
>
> I am very new to ITK and looking for help.
> I am trying to obtain hessian matrix values, in order to extract
> information about the detected vessels orientation. I'd be fine directly
> obtaining eigenvalues and eigenvectors.
> I have compiled ITK with python wrappers and the BridgeNumpy module.
>
> What I've managed to do so far:
>
> # read my image
> pixelType = itk.F
> imageType = itk.Image[pixelType, 3]
> readerType = itk.ImageFileReader[imageType]
> reader = readerType.New()
> reader.SetFileName("/path/to/my/image.mha")
> # compute Hessian matrix
> hessian_filter = itk.HessianRecursiveGaussianImageFilter[imageType].New()
> hessian_filter.SetInput(reader.GetOutput())
> # compute vesselness image
> vesselness_filter = itk.Hessian3DToVesselnessMeasureImageFilter[itk.F].New()
> vesselness_filter.SetInput(hessian_filter.GetOutput())
> vesselness_filter.Update()
> # convert vesselness image to numpy array
> itk_py_converter = itk.PyBuffer[imageType]
> vesselness_ndarray = itk_py_converter.GetArrayFromImage(vesselness_filter.GetOutput())
> # output: x,y,z numpy array corresponding to my input image, which is a good start
>
>
> Unfortunately, if I try to convert the Hessian matrix the same way:
>
> buf = itk.PyBuffer[itk.Image.SSRTD33]
> hess_matrix_ndarray = buf.GetArrayFromImage(hessian_filter.GetOutput())
>
> I run into this error:
>
> (...)
> KeyError: "itkTemplate : No template [<class 'itkImagePython.itkImageSSRTD33'>] for the itk::PyBuffer class"
>
> *Is this because I did not compile/install properly or because this type
> of pixel (itkImageSSRTD33) is not implemented in the PyBuffer module?*
>
> Anyway I tried a different approach to obtain hessian matrix values:
>
> hessian_filter.GetOutput().GetPixel([x,y,z]).GetNthComponent(n)
>
> which works but I don't know *in what order are** the 6 components* and
> didn't even find out *how to extract the image boundaries as python
> standard types?*
>
> In: hessian_filter.GetOutput().GetBufferedRegion()
> Out: itkImageRegion3([0, 0, 0], [480, 384, 255])  # this looks OK but...
> In: hessian_filter.GetOutput().GetBufferedRegion()[1]
> Out: TypeError: 'itkImageRegion3' object does not support indexing  # DAMN!
>
> As the vesselness filter performs eigenanalysis of the hessian matrix,*
> is there a way to extract the eigenvectors from there*?
>
> If you've read this far, thanks a lot. Please point me the right
> direction... I used SimpleITK before and was less confused but
> unfortunately it is my understanding that it is not possible to use hessian
> filters in SimpleITK at the moment.
>
> Regards,
>
> --
> Nicolas
>
> _______________________________________________
> Community mailing list
> Community at itk.org
> http://public.kitware.com/mailman/listinfo/community
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20170324/a6282c61/attachment.html>