From 09ef1b67657114545ccf0342d7b4a32dec0c63ec Mon Sep 17 00:00:00 2001 From: Philip Cook Date: Thu, 21 Mar 2024 14:55:43 -0400 Subject: [PATCH] ENH: Reset index after ITK filters, keep ITK object and python interface consistent --- ants/lib/LOCAL_antsImage.h | 31 +++++++++++++++++++++++++++++ ants/lib/LOCAL_cropImage.cxx | 20 ++----------------- ants/lib/LOCAL_padImage.cxx | 4 ++-- ants/utils/pad_image.py | 10 +++++----- tests/test_core_ants_image.py | 10 +++++----- tests/test_utils.py | 37 +++++++++++++++++++++++++++++++++++ 6 files changed, 82 insertions(+), 30 deletions(-) diff --git a/ants/lib/LOCAL_antsImage.h b/ants/lib/LOCAL_antsImage.h index 8cb89ee6..d2d83f86 100644 --- a/ants/lib/LOCAL_antsImage.h +++ b/ants/lib/LOCAL_antsImage.h @@ -345,5 +345,36 @@ py::tuple getSpacing( py::capsule & myPointer ) extern std::string ptrstr(py::capsule c); +/* +This function resets the region of an image to index from zero if needed. This +keeps the voxel indices in the numpy matrix consistent with the ITK image, and +also keeps the origin of physical space of the consistent with how it will be +saved as NIFTI. +*/ +template +static void FixNonZeroIndex( typename ImageType::Pointer img ) +{ + assert(img); + + typename ImageType::RegionType r = img->GetLargestPossibleRegion(); + typename ImageType::IndexType idx = r.GetIndex(); + + for (unsigned int i = 0; i < ImageType::ImageDimension; ++i) + { + // if any index is non-zero, reset the origin and region + if ( idx[i] != 0 ) + { + typename ImageType::PointType o; + img->TransformIndexToPhysicalPoint( idx, o ); + img->SetOrigin( o ); + + idx.Fill( 0 ); + r.SetIndex( idx ); + img->SetRegions( r ); + + return; + } + } +} #endif diff --git a/ants/lib/LOCAL_cropImage.cxx b/ants/lib/LOCAL_cropImage.cxx index 711cb6c6..c6120276 100644 --- a/ants/lib/LOCAL_cropImage.cxx +++ b/ants/lib/LOCAL_cropImage.cxx @@ -45,15 +45,7 @@ typename ImageType::Pointer cropImageHelper( typename ImageType::Pointer image, cropper->SetDirectionCollapseToSubmatrix(); cropper->UpdateLargestPossibleRegion(); cropper->GetOutput()->SetSpacing( image->GetSpacing() ); - typename ImageType::RegionType region = - cropper->GetOutput()->GetLargestPossibleRegion(); - typename ImageType::IndexType ind = region.GetIndex(); - typename ImageType::PointType neworig; - image->TransformIndexToPhysicalPoint( ind, neworig ); - ind.Fill(0); - region.SetIndex( ind ); - cropper->GetOutput()->SetRegions( region ); - cropper->GetOutput()->SetOrigin( neworig ); + FixNonZeroIndex( cropper->GetOutput() ); return cropper->GetOutput(); } return nullptr; @@ -97,15 +89,7 @@ typename ImageType::Pointer cropIndHelper( typename ImageType::Pointer image, cropper->SetDirectionCollapseToSubmatrix(); cropper->Update(); cropper->GetOutput()->SetSpacing( image->GetSpacing() ); - typename ImageType::RegionType region = - cropper->GetOutput()->GetLargestPossibleRegion(); - typename ImageType::IndexType ind = region.GetIndex(); - typename ImageType::PointType neworig; - image->TransformIndexToPhysicalPoint( ind, neworig ); - ind.Fill(0); - region.SetIndex( ind ); - cropper->GetOutput()->SetRegions( region ); - cropper->GetOutput()->SetOrigin( neworig ); + FixNonZeroIndex( cropper->GetOutput() ); return cropper->GetOutput(); } return nullptr; diff --git a/ants/lib/LOCAL_padImage.cxx b/ants/lib/LOCAL_padImage.cxx index b2e7d763..6cb2d4ee 100644 --- a/ants/lib/LOCAL_padImage.cxx +++ b/ants/lib/LOCAL_padImage.cxx @@ -14,7 +14,7 @@ namespace py = pybind11; template < typename ImageType > -py::capsule padImage( py::capsule & antsImage, +py::capsule padImage( py::capsule & antsImage, std::vector lowerPadDims, std::vector upperPadDims, float padValue ) @@ -47,7 +47,7 @@ py::capsule padImage( py::capsule & antsImage, padFilter->SetPadUpperBound( upperExtendRegion ); padFilter->SetConstant( padValue ); padFilter->Update(); - + FixNonZeroIndex( padFilter->GetOutput() ); return wrap< ImageType >( padFilter->GetOutput() ); } diff --git a/ants/utils/pad_image.py b/ants/utils/pad_image.py index 6116ee6b..3204c39e 100644 --- a/ants/utils/pad_image.py +++ b/ants/utils/pad_image.py @@ -19,13 +19,13 @@ def pad_image(image, shape=None, pad_width=None, value=0.0, return_padvals=False shape : tuple - if shape is given, the image will be padded in each dimension until it has this shape - - if shape is not given, the image will be padded along each - dimension to match the largest existing dimension so that it - has isotropic dimension + - if shape and pad_width are both None, the image will be padded along + each dimension to match the largest existing dimension so that it has + isotropic dimensions. pad_width : list of integers or list-of-list of integers How much to pad in each direction. If a single list is - supplied (e.g., [4,4,4]), then the image will be padded by half + supplied (e.g., [4,4,4]), then the image will be padded by half that amount on both sides. If a list-of-list is supplied (e.g., [(0,4),(0,4),(0,4)]), then the image will be padded unevenly on the different sides @@ -72,7 +72,7 @@ def pad_image(image, shape=None, pad_width=None, value=0.0, return_padvals=False libfn = utils.get_lib_fn('padImageF%i' % ndim) itkimage = libfn(image.pointer, lower_pad_vals, upper_pad_vals, value) - new_image = iio.ANTsImage(pixeltype='float', dimension=ndim, + new_image = iio.ANTsImage(pixeltype='float', dimension=ndim, components=image.components, pointer=itkimage).clone(inpixeltype) if return_padvals: return new_image, lower_pad_vals, upper_pad_vals diff --git a/tests/test_core_ants_image.py b/tests/test_core_ants_image.py index 2c7cc029..57038990 100644 --- a/tests/test_core_ants_image.py +++ b/tests/test_core_ants_image.py @@ -235,7 +235,7 @@ def test__add__(self): def test__radd__(self): if os.name == "nt": return - + #self.setUp() for img in self.imgs: # op on constant @@ -253,7 +253,7 @@ def test__radd__(self): img2 = img.clone() img2.set_spacing([2.31]*img.dimension) img3 = img + img2 - + def test__sub__(self): #self.setUp() for img in self.imgs: @@ -276,7 +276,7 @@ def test__sub__(self): def test__rsub__(self): if os.name == "nt": return - + #self.setUp() for img in self.imgs: # op on constant @@ -294,7 +294,7 @@ def test__rsub__(self): img2 = img.clone() img2.set_spacing([2.31]*img.dimension) img3 = img - img2 - + def test__mul__(self): #self.setUp() for img in self.imgs: @@ -335,7 +335,7 @@ def test__rmul__(self): img2 = img.clone() img2.set_spacing([2.31]*img.dimension) img3 = img * img2 - + def test__div__(self): #self.setUp() for img in self.imgs: diff --git a/tests/test_utils.py b/tests/test_utils.py index 8a4a350b..dfa1fefd 100644 --- a/tests/test_utils.py +++ b/tests/test_utils.py @@ -237,6 +237,43 @@ def test_decrop_image_example(self): # full image not float cropped = ants.crop_image(fi, mask.clone("unsigned int"), 1) +class TestModule_pad_image(unittest.TestCase): + def setUp(self): + self.img2d = ants.image_read(ants.get_ants_data("r16")) + self.img3d = ants.image_read(ants.get_ants_data("mni")).resample_image((2, 2, 2)) + + def tearDown(self): + pass + + def test_pad_image_example(self): + img = self.img2d.clone() + img.set_origin((0, 0)) + img.set_spacing((2, 2)) + # isotropic pad via pad_width + padded = ants.pad_image(img, pad_width=(10, 10)) + self.assertEqual(padded.shape, (img.shape[0] + 10, img.shape[1] + 10)) + for img_orig_elem, pad_orig_elem in zip(img.origin, padded.origin): + self.assertAlmostEqual(pad_orig_elem, img_orig_elem - 10, places=3) + + img = self.img2d.clone() + img.set_origin((0, 0)) + img.set_spacing((2, 2)) + img = ants.resample_image(img, (128,128), 1, 1) + # isotropic pad via shape + padded = ants.pad_image(img, shape=(160,160)) + self.assertSequenceEqual(padded.shape, (160, 160)) + + img = self.img3d.clone() + img = ants.resample_image(img, (128,160,96), 1, 1) + padded = ants.pad_image(img) + self.assertSequenceEqual(padded.shape, (160, 160, 160)) + + # pad only on superior side + img = self.img3d.clone() + padded = ants.pad_image(img, pad_width=[(0,4),(0,8),(0,12)]) + self.assertSequenceEqual(padded.shape, (img.shape[0] + 4, img.shape[1] + 8, img.shape[2] + 12)) + self.assertSequenceEqual(img.origin, padded.origin) + class TestModule_denoise_image(unittest.TestCase): def setUp(self):