Single-precision displacement field transforms: a float32 field pays a float64 conversion today

[vboussot]

SimpleITK wraps only TransformBaseTemplate<double>, so a float32 displacement field handed to DisplacementFieldTransform is converted to float64 and interleaved at that width before anything is resampled.

At the scale of registration fields this is the dominant cost of the operation. Measured on a 122-row region of a large field (SimpleITK 2.5.5): building the double-precision transform allocates 6.9 GiB and takes 5.0 s, which is 39% of everything that region costs end to end in our pipeline. The field itself is float32 on disk and float32 is plenty for the geometry.

sitk.Warp escapes because it is templated on the field's own pixel type, but it only covers the narrowest case: field grid == output grid, whole map, linear interpolation. Any composite or resampled-field use pays the float64 build.

ITK already instantiates DisplacementFieldTransform<float, N>. Exposing it (or a Resample overload accepting a sitkVectorFloat32 field directly) would halve the memory and remove the conversion for every consumer that resamples through stored fields. Happy to help test a PR if the wrapping approach is agreed on.

[blowekamp]

Thank you for reporting this issue. Can you please provide a gist, of these cases to have concrete code to work with and compare?

[vboussot]

Here is a gist reproducing the cases: https://gist.github.com/vboussot/6d74bae6f5d27bd0422c38cb9c6898cb

It starts from one float32 3-component field and does three things (measured here on SimpleITK 2.5.5, a 512^3 x 3 field):

1. DisplacementFieldTransform(field32) throws immediately: Expected input displacement field image must be of pixel type: vector of 64-bit float. So double precision is a hard API requirement, not an internal preference.
2. The only way through is Cast(field, sitkVectorFloat64) first, then build the transform and resample: +3.50 GiB peak / 0.95 s, and the cast doubles the field buffer from 1.50 to 3.00 GiB.
3. sitk.Warp on the same float32 field keeps it float32: +0.50 GiB / 0.52 s, but Warp only covers the case where the field is the output grid (order 1, whole map, nothing to compose), so any composite or resampled-field use is stuck on path 2.

ITK already instantiates DisplacementFieldTransform<float, N>; exposing it through SimpleITK (or a Resample overload accepting a sitkVectorFloat32 field) would remove the mandatory cast and halve the field buffer for that path. This is the same double-precision coercion you traced on the file-read side in InsightSoftwareConsortium/ITK#6830. Happy to help test a wrapping.