Jacobian test for non-Rigid registration

[mfizyczka]

Hi, In the example python code jacobian_structure_mask.py there is this line: isNonRigid = abs(math.log(abs(jacobian.data[i]))) >= 0.5

Could you please explain why 0.5 is used here as a threshold? And why do you use log here instead of simply looking on jacobian.data values different from 1 when there is some volume change?

According to Velocity 4.1 Instructions of use:

10^0.5 = 3.16 and 10^(-0.5) = 0.316 meaning 3 times decrease / increase of voxel volume which I would find a bit surprising for a rigid registration.

[mliang1987]

Hello!

Thank you for your question. First, I want to say we definitely could have made the IFU portion clearer here. It should say that the value we have is the determinant of the Jacobian matrix. We will fix this in a future release of Velocity.

To understand what that particular line does, let's look at it in context:

print('Creating mask from threshold... ')
for i in range(0,volumetricData.size[0]*volumetricData.size[1]*volumetricData.size[2]):
  isDense = volumetricData.data[i] >= thresholdValue
  isNonRigid = abs(math.log(abs(jacobian.data[i]))) >= 0.5

  # find dense material like bone that had non-rigid deformation
  mask.data[i] = isDense and isNonRigid
print('done')

From here, it looks like the intention of this mask is to identify voxels that are likely bony anatomy (isDense) and had a non-rigid registration (isNonRigid), indicating problem issues with the registration vis-a-vis bony anatomy. So, it's not that we're looking for a rigid registration, as rigid registrations have uniform Jacobian determinants (of 1!) for each and every voxel, but rather we're looking to see if a voxel that should not have changed much has indeed changed too much in the deformation.

The line:

isNonRigid = abs(math.log(abs(jacobian.data[i]))) >= 0.5

Is finding the absolute value of the natural log of the absolute value of the Jacobian determinant, which means the range of abs(jacobian.data[i] is actually from 0 < x <= 0.61 and x >= 1.65. This can be read as the voxel has either significantly shrunk or has significantly increased. However, this range is an artificial range constructed for convenience to reduce it to a single comparison. You may wish to determine your own ranges here, such as:

absLogJacobDet = abs(math.log(abs(jacobian.data[i])))
isNonRigid = (0 < absLogJacobDet  <= 0.5) or (absLogJacobDet  >= 1.5)

or something a bit more complex.

I hope this helps, and please follow-up with any questions you may have.

Michael Liang Software Engineer Varian Medical Systems

[mfizyczka]

Thank you for your response and explanation. Just to be sure - what is IFU?

[mliang1987]

Hey, sorry for the late reply. IFU here means Instructions-for-Use.