Merge pull request #588 from ANTsX/reset_index_after_filter

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

Author: cookpa

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