Mapping of displacement fields

[jbanusco]

Hello, I believe there is a bug during the mapping of deformation fields.

In the .map method of the DenseFieldTransform class, I believe the locations that are mapped out-of-domain should be set to the initial coordinates instead of 0. Otherwise, if for example you want to go back to the displacement field by subtracting the reference coordinates you end-up with very large displacements instead of 0.

def map(self, x, inverse=False):
        ....
        ....
        new_map = np.vstack(tuple(
            map_coordinates(
                self._field[..., i],
                ijk.T,
                order=1,
                mode="constant",
                cval=0,
                prefilter=False,
            ) for i in range(self.reference.ndim)
        )).T

        # Now, where is 0 i set the original coordinates value -- no transformation/displacement
        new_map[np.isclose(new_map, 0)] = x[np.isclose(new_map, 0)]
        return new_map

Thanks for your work!

[effigies]

Yes, I think this is correct. I wonder if a more direct approach would be:

         return np.vstack(tuple(
-            map_coordinates(
-                self._field[..., i],
+            ijk.T + map_coordinates(
+                self._deltas[..., i],
                 ijk.T,
                 order=3,
                 mode="constant",
                 cval=0,
                 prefilter=True,
             ) for i in range(self.reference.ndim)
         )).T

By interpolating the deltas, we avoid needing a second step. We could also separate displacement fields from deformation fields, and treat displacement fields as I do and deformation fields as you do, which would prevent adding unnecessary floating point error to the mix.

[oesteban]

The rationale to adopt deformation fields (instead of displacements fields or "deltas") is float resolution. Deformation fields (SPM's preferred option) are more accurate than displacements fields (ANTs, FSL).

BTW - @jbanusco has fixed this issue (and tested himself) here - https://github.com/nipy/nitransforms/compare/master...jbanusco:nitransforms:master

I'm in the process of getting him to submit those improvements as separate PRs... :D

[effigies]

Deformation fields (SPM's preferred option) are more accurate than displacements fields (ANTs, FSL).

Is this still true when you start by converting from deltas to deformations? If you start with a deformation, then yes, you get:

interpolate(deformation, ijk)
# vs
ijk + interpolate(deformation - ijk, ijk)

If you start with deltas, then you get

interpolate(ijk + deltas, ijk)
# vs
ijk + interpolate(deltas, ijk)

I'm not sure why doing the addition in one place vs the other would accumulate more floating point error.