Here’s what I found when looking at the displacement fields.
The displacement field mf_g1.nii
had one slice with several displacement values with a magnitude of 1e6 or higher. I’m guessing that these values were errors of some kind? They were far from the “cardiac box”, so I set them to zero and then the inversion worked fine.
The displacement field mf_g2.nii
was more interesting. At one edge of the “cardiac box,” the Jacobian determinant was negative. Essentially, the displacement field folds over itself at this edge and becomes non-invertible. More on that later, but first here is the Python code that I used for the check:
# use vtkImagePointDataIterator to iterate through the displacement field
image_iter = vtk.vtkImagePointIterator(image_data_sptr)
# double jac[3][3]
jac = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]
# double p2[3]
p2 = [0.0, 0.0, 0.0]
# iterate through all sample points in the image data
while not image_iter.IsAtEnd():
# get current position in image (in mm)
p1 = image_iter.GetPosition()
# evaluate Jacobian of bspline_transform at this point
bspline_transform.InternalTransformDerivative(p1, p2, jac)
# determinant
det = (jac[0][0]*jac[1][1]*jac[2][2] +
jac[1][0]*jac[2][1]*jac[0][2] +
jac[2][0]*jac[0][1]*jac[1][2] -
jac[0][0]*jac[2][1]*jac[1][2] -
jac[1][0]*jac[0][1]*jac[2][2] -
jac[2][0]*jac[1][1]*jac[0][2])
if det < 0.0:
idx = image_iter.GetIndex()
print(idx, p1, det)
image_iter.Next()
So more about the fold. Everywhere outside of the “cardiac box”, the displacements are set to zero. The problem is this: at the edge of the box there are outward displacements on the order of 3 mm, but the sample spacing is only 2 mm. Since the displacement drops from 3 mm to 0 mm over a distance of 2 mm, there is a fold and therefore the one-to-one mapping between the original coordinate space and the deformed coordinate space is lost.
It looks to me that the “cardiac box” was too small when the measurements were done, resulting in displacements too large at the edge of the box. If you modify the data so that the displacements decay towards zero more gradually, the folding can be eliminated.
If the lack of invertibility at the edge doesn’t bother you, then you can just turn off the warnings when you invert the transform:
vtkObject::GlobalWarningDisplayOff();
// compute inverse grid here
vtkObject::GlobalWarningDisplayOn();