Merge pull request #588 from ANTsX/reset_index_after_filter

ENH: Reset index after ITK filters, keep ITK object and python interface consistent

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 <typename ImageType>
+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

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<ImageType>( 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<ImageType>( cropper->GetOutput() );
  return cropper->GetOutput();
  }
  return nullptr;

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<int> lowerPadDims,
  std::vector<int> upperPadDims,
  float padValue )
@@ -47,7 +47,7 @@ py::capsule padImage( py::capsule & antsImage,
  padFilter->SetPadUpperBound( upperExtendRegion );
  padFilter->SetConstant( padValue );
  padFilter->Update();
-
+ FixNonZeroIndex<ImageType>( padFilter->GetOutput() );
  return wrap< ImageType >( padFilter->GetOutput() );
 }
 

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

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:

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):