Hi All,
I’m working on implementing 2D-3D deformable registration in ITK (v5.3.0) using a BSplineTransform to fit deformable motion, but have hit a wall optimising a BSplineTransform. The approach I’ve taken is to compose a composite transform consisting of an optimisable BSplineTransform (modelling deformable motion) and a non-optimisable Euler3D transform (for raycasting). I’ve followed along with some 2D-3D examples for rigid registration (e.g. TwoProjection2D3DRegistation), and have successfully optimised a rigid transform instead of a BSpline transform using gradient based optimisers. However, when setting the optimisable transform as a BSplineTransform, the registration/optimisation throws errors that differ with the optimiser/registration method used.
I’ve gotten the following errors:
For combination of itkImageRegistrationMethodv4 and itkLBFGSOptimizerv4: “ComputeJacobianWithRespectToPosition not yet implemented.”
For combination of itkImageRegistrationMethod and itkLBFGSOptimizer: “THE SEARCH DIRECTION IS NOT A DESCENT DIRECTION. IFLAG= -1 LINE SEARCH FAILED.”
For combination itkImageRegistrationMethod and itkRegularStepGradientDescentOptimizer: “Gradient magnitude tolerance met after 0 iterations. Gradient magnitude (0) is less than gradient magnitude tolerance (0.0001).”
N.B. Non-gradient-based optimisers (AmoebaOptimizerv4/PowellOptimizerv4) have not thrown errors, but are prohibitively slow - likely due to the combination of the relatively large number of BSpline parameters (approx 1500) that need to be optimised and ray-casting interpolator. Can anyone offer any advice/solution with regards to these errors? Thank you in advance!
Hi Dženan, thanks for your prompt reply! After some reflection and further testing, I’ve found that the gradients being very close to 0 is quite reasonable - not an error. In my case, I was registering abdominal CBCTs acquired in breath-hold to intra-fraction radiographs acquired in breath-hold (from the same treatment fraction).
In the previous tests, the BSpline Transform was quite coarse and initialised with all parameters set to 0. Given both images are in breath hold and acquired during the same fraction, it’s not surprising that the gradients are quite small - initialising the parameters using the NormalVariateGenerator, the gradients were no longer below the minimum threshold, but changes in parameter values were only very small. On another note, I was using correlation as the metric, and suspect this could also limit gradient magnitudes.
Thanks for your help!
For any that come across this issue in the future, further testing has revealed that the parameter values were not responsible (not sure what happened there - likely user error). Instead, it looks like the itkRayCastInterpolateImageFunction might be the culprit. When resampling the image, the interpolator applies the TransformPoint method to the focal point (set to be 1000 mm from isocentre in my case) to translate/rotate the point about the isocentre. No actual image resampling takes place prior to the raycasting - deformation of the moving 3D image is ignored. I suspect a custom interpolator is needed for 2D-3D deformable registration registration.
Another quick update after some more tests. I suspect that deformable 2D to 3D registration is not feasible with a standard BSpline transform (even with a custom interpolator). In the ComputeJacobianWithRespectToParameters() method, the Jacobian is evaluated in the fixed image domain (kV image), set to 0 when not within the transform’s valid region. Unfortunately, the fixed image domain can’t physically overlap with the moving image (CT / CBCT), otherwise detector panels would be invading patients.
All thoughts/ideas are welcome on potential solutions/workarounds.
Instead of using the real position of the detector panel, you could scale the distance (and accordingly size) of the panel from the source by a fixed amount, e.g. 90% or 80%. After the optimization is finished, scale it back to full distance/size. Would that work?
Hi Dženan, thanks for your support!
I’ve found solutions to the problems I was having earlier. To incorporate deformation during the raycasting process, I made a modified version of the RayCastInterpolateImageFunction that applies the transform’s TransformPoint method at each voxel sampled along the ray. To solve the zero gradients issue during registration, I made a modified version of NormalizedCorrelationImageToImageMetric that integrates the gradient image / Jacobians along the rayline instead of evaluating it at the detector pixels. Both processes can be quite slow with BSplineTransforms, but they do work. I’ve got some further testing to run to verify that the gradient calculated during the registration is correct. I’d be happy to package these methods into a module in the future. Thanks for your help!
