Smoothness of the deformation fields using the greedy method

[anorak94]

Hey I was testing your package on some simple images . After running the algorithm i checked out the warp field

And it looks something like this. I know that the greedy method does not guarantee that there should be no folds in the grid do you think this is the expected result.

I tried to implement the greedy method using SyN in jax. Although it is far from perfect right now the results look more reasonable especially the warping grid

https://colab.research.google.com/drive/1ZJYSi245WP3RMREdyxGGr1UsEYjeMWnS

But thanks for this package. This and pyrpl have really helped me with the more obscure bits such as resampling etc of displacement fields.

[anorak94]

Also just a small doubt when we are warping the images in transformer.apply transform. Why do we pass in

np.meshgrid(np.arange(w), np.arange(h)) + disp / spacing instead of np.meshgrid(np.arange(w), np.arange(h)) + disp as the identity_mappingand disp are discretised on the same size grid

[GFleishman]

Hi - I must have missed the github email about these messages, but happy to find them now. Glad that greedypy and pyrpl have been helpful for your own development! I don't use them very much anymore and they were never given the development they need for widespread use. I'm glad someone is able to look into the source code and find useful things there. Regarding your specific questions:

Based on the warp you shared, I'm assuming the circle was the fixed image and the square was the moving image. I can't say for sure without seeing the resampled image, but it seems like an accurate transform, but very non-smooth. As you say, there is nothing in the greedy algorithm which guarantees smoothness - how smooth your result is depends primarily on the regularization and gradient descent step size parameters. Circle to square is a common example I've run myself many times achieving an accurate result with a smooth transform. So either over the (admittedly long) period since I've used this code something has changed and/or your example alignment may not be parameterized well.

In greedypy there are two regularizations: smoothing of the gradient term at each iteration and smoothing of the total field at each iteration. These smoothings are parameterized via: grad_regularizer and field_regularizer which in some places are also referred to as gradient_abcd and field_abcd. Each of these parameters is a list of four numbers, [a, b, c, d] parameterizing the following differential smoothing kernel (the inverse of which is applied in Fourier space): (-a laplacian + b gradient-of-divergence + c) ** d.

If you run greedypy via greedypy_registration_method.py then the defaults are field_abcd = [0.5, 0.0, 1.0, 6.0] and gradient_abcd = [3.0, 0.0, 1.0, 2.0]. If you run via the command line interface then the defaults are field_regularizer = [1, 0, 1, 2] and grad_regularizer = [3, 0, 1, 2].

a functions more or less like the sigma of a Guassian smoothing kernel, b should always be set to 0 as I never implemented the "elastic" gradient-of-divergence term, c should basically always be set to 1 as it just prevents a divide by zero from happening, and d functions more or less like the number of times you apply a Gaussian smoothing kernel in a row.

Apart from smoothing parameters, the gradient descent step size will impact smoothness. The larger this term, the larger the magnitude of the update field - hence a more powerful smoothing would be required to straighten out any crossing displacement vectors.

[GFleishman]

Regarding your own implementation, I'm confused when you say "implement greedy using SyN" as these are two different algorithms. Greedy simply updates the forward transform each iteration with the smoothed transform gradient of the objective function. SyN is much more laborious as it simultaneously optimizes forward and inverse initial velocities ensuring they are consistent when integrated to the midpoint.

Such care in ensuring not only smoothness but numerically accurate invertibility is really only necessary if the absolute shape of objects is important (such as a brain morphometry study). But if you just want to make sure the right things are on top of each other in image A and image B, and the absolute shape of things is not so important, the extra work of SyN is not a good idea.

I'd be happy to see your results, but the link you shared is to a private google drive.

[GFleishman]

Regarding your question about applying transforms - displacement vector fields should always be encoded in physical units (e.g. millimeters), as these have real meaning that is not tied to a specific voxel grid (as voxel units would be). But the second argument to map_coordinates must be in voxel units. So, disp / spacing just converts the displacement vector field to voxel coordinates.

[anorak94]

hey thanks for your extensive comments. Regarding the warp field it was my mistake. The function i wrote for visualization expects input in the format CxHxW and i was passing in the wrong format. The field is indeed smooth.

Regarding SyN - the method implemented in

https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Registration/RegistrationMethodsv4/include/itkSyNImageRegistrationMethod.hxx

is slightly different from the one described in the original paper. The algorithm which they use to set benchmarks is the Greedy version of SyN which updates the displacement transform for fixed to middle and middle to fixed with the gradient of the objective function and then inverts the field to get the complete warp atleast thats my understanding of the algorithm.

I asked this question on the SyN website as well

https://github.com/ANTsX/ANTsPy/issues/313

I am really interested in the way smoothing computations are done in the fourier domain so far i havent found any good references for those. I also implemented lddmm using the vector momentum parameterization but in jax using autodiff for gradient computations. I was wondering what your experience with time parameterized deformable algorithms has been. Are they more accurate than greedy methods

Regarding the implementation i will make it public

[GFleishman]

I'm usually happy to chat about registration stuff :)

Another tool you might find useful for visualizing transforms is itk-snap. You can load vector fields and view as a deformation grid. Here's an example:

I see what you mean now about Greedy-SyN. I think there are 3 algorithms here: Greedy (as in this very good package), Greedy-SyN, and SyN (described in the original Avants paper).

Regarding smoothing in the Fourier domain, consider a simple 1D derivative of an image. That can be implemented by convolving with a finite difference kernel. A complicated differential operator like the simple-elastic-operator I mentioned above (and used in many registration packages) is just a bunch of simple 1D derivatives combined; so that also can be implemented by convolving with some (more complicated) finite difference kernel (e.g. here's the Laplacian finite difference formula).

Convolution in the spatial domain is equivalent to multiplication in the Fourier domain - so you can also achieve these derivatives by multiplying the Fourier transform of the image by the Fourier transform of the finite difference kernel. So to apply the inverse of a differential operator, you can multiply the Fourier transform of the image by the inverse of the operator's Fourier transform.

There is an analytical formula for the DFT of the Laplacian finite different operator. You can find it (and a more mathematical description of what I hand-waved above) in section 3.2 of the original LDDMM paper.

If what you mean by accurate is whether or not the algorithm puts the most likely corresponding objects into the same coordinates, then personally - I do not believe that time parameterized registration algorithms (LDDMM, geodesic shooting, stationary velocity fields) are any more accurate than algorithms that go strait for learning a single transform. The advantage of diffeomorphic algorithms is they provide a much stronger guarantee that the transform will be smooth (though that guarantee still breaks down due to discretization - you can get SyN or LDDMM to give you non-diffeomorphic fields if you parameterize them badly).

If you believe that more strictly limiting your optimization to the space of smoother transforms will enable you to find a more accurate alignment, then maybe it's worth it, but in for many simple image matching problems it just is not necessary.

[GFleishman]

After thinking on this a little more I want to make one more comment.

In addition to stronger smoothness guarantees, the other advantage that strongly diffeomorphic methods (i.e. time parameterized flows) have, which is probably more important than the smoothness thing, is that they are more informative models. With a simple registration (find a field matching two datasets) you are only establishing correspondence, but you're not trying to model how those two configurations transitioned from one arrangement to the other. With time parameterized flow you are modeling that change.

As an example, consider two points in 2D Euclidean space. On one hand, you can look at the difference between these two points, which is a single measurement (dx, dy). On the other hand, you can fit a line between those two points (slope, intercept) which could be used for interpolation/extrapolation and give more insight on the relationship between those two points. If you have more than 2 data points, it becomes even more useful to use a model instead of just measuring pairwise differences. Time parameterized geodesic diffeomorphic registration algorithms are the generalization of linear regression to the space of smooth transforms.

Additional models exist as well, if you happen to know that the time evolution of your shape is not geodesic. Here is a great paper on higher order regression models in diffeomorphic space.

Such methods are very theoretically advanced, so much so that the original developers spend more time thinking about the theory than the implementation. As well, such a method is mathematically beyond the typical computer science PhD - so it doesn't see much practical engineering work (good software writing) either. None the less - these are very cool and powerful ideas waiting for the right implementation and applications.