resample_bounding_box sizes a region by walking the boundary of the transformed grid. For a linear map that is exact. For a non-linear one it reports what the map does to the boundary and nothing about what it does strictly inside, so a displacement that is zero on the grid's boundary voxels and large in the middle leaves the region unchanged. resample then crops the moving image to that region, and the pixels the bump displaces onto come back as the default value.
Measured
2D, fixed grid 64x64, moving 256x256, unit scale, zero translation, padding 1. A raised-cosine bump of radius 20 centred on the grid, peak displacement (y, x) = (+60, -60), exactly zero on the boundary voxels. Handed in as an itk.DisplacementFieldTransform[itk.F, 2]:
itk transform maps fixed point (x=32, y=32) -> (x=-28.000, y=92.000)
itk transform maps fixed point (x=0, y=0) -> (x=0.000, y=0.000)
ITK-transform branch : start_index {'y': -1, 'x': -1} size {'y': 66, 'x': 66}
identity transform : start_index {'y': -1, 'x': -1} size {'y': 66, 'x': 66}
The region for the field is the region for the identity. The true reach is [0, 92.6] on y and [-28.6, 63] on x, so it misses 28.6 px above on y and 27.6 px below on x.
End to end, resample against a single whole-image ITK call: 63 of 4096 output pixels differ, by up to 213.5. The worst is (y=40, x=39), where resample gives 0 and the whole-image call gives 213.5, reading moving index (67.0, 12.0), outside the reported region.
The same without an interpolation oracle: overwrite every moving pixel outside the reported region with a sentinel and resample again. 64 of 4096 output pixels change, so the resample provably reads outside the region it reported.
Moving the bump onto a boundary row makes the region a superset again, which is what identifies the boundary walk as the cause. It reproduces in 3D (32^3 fixed, 128^3 moving, amplitude 30: misses 14.3 px on z, 13.3 on y, 14.3 on x).
The RFC-5 field path is fixed; this is the ITK one
#695 fixes this for an RFC-5 displacements or coordinates transform, by reading the range of the field's values and moving and widening the region by it. A field is read through a kernel that is non-negative and sums to one, so a displacement anywhere lies between the smallest and largest of the values around it: a bound rather than a sample of the boundary.
An ITK transform can carry the same bound, and from data already in hand:
- a
DisplacementFieldTransform holds its field in parameters, so the range is one pass over that array;
- a
BSplineTransform holds coefficients on a control grid, and the cubic B-spline basis is non-negative and sums to one, so the same argument applies to them;
- a
CompositeTransform folds the bounds of its stages, an affine stage carrying an inner residual through |A| @ residual.
Everything else stays as it is: a linear ITK transform is already exact under the boundary walk.
Affects
itk_transform_resample and itk_transform_resample_bounding_box on main today, for any non-linear transform, which includes the CompositeTransform an Elastix registration returns. The TypeScript resampleBoundingBox has the same walk.
resample_bounding_boxsizes a region by walking the boundary of the transformed grid. For a linear map that is exact. For a non-linear one it reports what the map does to the boundary and nothing about what it does strictly inside, so a displacement that is zero on the grid's boundary voxels and large in the middle leaves the region unchanged.resamplethen crops the moving image to that region, and the pixels the bump displaces onto come back as the default value.Measured
2D, fixed grid 64x64, moving 256x256, unit scale, zero translation, padding 1. A raised-cosine bump of radius 20 centred on the grid, peak displacement (y, x) = (+60, -60), exactly zero on the boundary voxels. Handed in as an
itk.DisplacementFieldTransform[itk.F, 2]:The region for the field is the region for the identity. The true reach is
[0, 92.6]on y and[-28.6, 63]on x, so it misses 28.6 px above on y and 27.6 px below on x.End to end,
resampleagainst a single whole-image ITK call: 63 of 4096 output pixels differ, by up to 213.5. The worst is(y=40, x=39), whereresamplegives 0 and the whole-image call gives 213.5, reading moving index(67.0, 12.0), outside the reported region.The same without an interpolation oracle: overwrite every moving pixel outside the reported region with a sentinel and resample again. 64 of 4096 output pixels change, so the resample provably reads outside the region it reported.
Moving the bump onto a boundary row makes the region a superset again, which is what identifies the boundary walk as the cause. It reproduces in 3D (32^3 fixed, 128^3 moving, amplitude 30: misses 14.3 px on z, 13.3 on y, 14.3 on x).
The RFC-5 field path is fixed; this is the ITK one
#695 fixes this for an RFC-5
displacementsorcoordinatestransform, by reading the range of the field's values and moving and widening the region by it. A field is read through a kernel that is non-negative and sums to one, so a displacement anywhere lies between the smallest and largest of the values around it: a bound rather than a sample of the boundary.An ITK transform can carry the same bound, and from data already in hand:
DisplacementFieldTransformholds its field inparameters, so the range is one pass over that array;BSplineTransformholds coefficients on a control grid, and the cubic B-spline basis is non-negative and sums to one, so the same argument applies to them;CompositeTransformfolds the bounds of its stages, an affine stage carrying an inner residual through|A| @ residual.Everything else stays as it is: a linear ITK transform is already exact under the boundary walk.
Affects
itk_transform_resampleanditk_transform_resample_bounding_boxonmaintoday, for any non-linear transform, which includes theCompositeTransforman Elastix registration returns. The TypeScriptresampleBoundingBoxhas the same walk.