2

What I am trying to do is take a numpy array representing 3D image data and calculate the hessian matrix for every voxel. My input is a matrix of shape (Z,X,Y) and I can easily take a slice along z and retrieve a single original image. 

gx, gy, gz = np.gradient(imgs)

gxx, gxy, gxz = np.gradient(gx)
gyx, gyy, gyz = np.gradient(gy)
gzx, gzy, gzz = np.gradient(gz)

And I can access the hessian for an individual voxel as follows:

x = 100
y = 100
z = 63

H = [[gxx[z][x][y], gxy[z][x][y], gxz[z][x][y]],
     [gyx[z][x][y], gyy[z][x][y], gyz[z][x][y]],
     [gzx[z][x][y], gzy[z][x][y], gzz[z][x][y]]]

But this is cumbersome and I can't easily slice the data.

I have tried using reshape as follows 

H = H.reshape(Z, X, Y, 3, 3)

But when I test this by retrieving the hessian for a specific voxel the, the value returned from the reshaped array is completely different than the original array.

I think I could use zip somehow but I have only been able to find that for making lists of tuples.

  • Bonus: If there's a faster way to accomplish this please let me know, I essentially need to calculate the three eigenvalues of the hessian matrix for every voxel in the 3D data set. Calculating the hessian values is really fast but finding the eigenvalues for a single 2D image slice takes about 20 seconds. Are there any GPUs or tensor flow accelerated libraries for image processing?
dadrake
  • 112
  • 7
  • [Shameless plug] You should try using [PyDIP](https://github.com/DIPlib/diplib). There you can do `H = dip.Hessian(img); D,V = dip.EigenDecomposition(H)`, returning a diagonal matrix image `D` with eigenvalues and a matrix image `V` with eigenvectors. In these images, `H[x,y,z]` returns the Hessian matrix for the pixel at indices `(x,y,z)`. – Cris Luengo Nov 15 '18 at 06:59

1 Answers1

2

We can use a list comprehension to get the hessians -

H_all = np.array([np.gradient(i) for i in np.gradient(imgs)]).transpose(2,3,4,0,1)

Just to give it a bit of explanation : [np.gradient(i) for i in np.gradient(imgs)] loops through the two levels of outputs from np.gradient calls, resulting in a (3 x 3) shaped tensor at the outer two axes. We need these two as the last two axes in the final output. So, we push those at the end with the transpose.

Thus, H_all holds all the hessians and hence we can extract our specific hessian given x,y,z, like so -

x = 100
y = 100
z = 63
H = H_all[z,y,x]
Divakar
  • 218,885
  • 19
  • 262
  • 358