Use double precision by default in apply_transforms computations

ants/registration/apply_transforms.py

@@ -11,7 +11,7 @@
 def apply_transforms(fixed, moving, transformlist,
  interpolator='linear', imagetype=0,
  whichtoinvert=None, compose=None,
- defaultvalue=0, verbose=False, **kwargs):
+ defaultvalue=0, singleprecision=False, verbose=False, **kwargs):
  """
  Apply a transform list to map an image from one domain to another.
  In image registration, one computes mappings between (usually) pairs
@@ -67,6 +67,10 @@ def apply_transforms(fixed, moving, transformlist,
  defaultvalue : scalar
  Default voxel value for mappings outside the image domain.
 
+ singleprecision : boolean
+ if True, use float32 for computations. This is useful for reducing memory
+ usage for large datasets, at the cost of precision.
+
  verbose : boolean
  print command and run verbose application of transform.
 
@@ -102,16 +106,18 @@ def apply_transforms(fixed, moving, transformlist,
 
  args = [fixed, moving, transformlist, interpolator]
 
+ output_pixel_type = 'float' if singleprecision else 'double'
+
  if not isinstance(fixed, str):
  if isinstance(fixed, iio.ANTsImage) and isinstance(moving, iio.ANTsImage):
  for tl_path in transformlist:
  if not os.path.exists(tl_path):
  raise Exception('Transform %s does not exist' % tl_path)
 
  inpixeltype = fixed.pixeltype
- fixed = fixed.clone('float')
- moving = moving.clone('float')
- warpedmovout = moving.clone()
+ fixed = fixed.clone(output_pixel_type)
+ moving = moving.clone(output_pixel_type)
+ warpedmovout = moving.clone(output_pixel_type)
  f = fixed
  m = moving
  if (moving.dimension == 4) and (fixed.dimension == 3) and (imagetype == 0):
@@ -165,7 +171,7 @@ def apply_transforms(fixed, moving, transformlist,
  if verbose:
  print(myargs)
 
- processed_args = myargs + ['-z', str(1), '-v', str(myverb), '--float', str(1), '-e', str(imagetype), '-f', str(defaultvalue)]
+ processed_args = myargs + ['-z', str(1), '-v', str(myverb), '--float', str(int(singleprecision)), '-e', str(imagetype), '-f', str(defaultvalue)]
  libfn = utils.get_lib_fn('antsApplyTransforms')
  libfn(processed_args)
 
@@ -180,7 +186,7 @@ def apply_transforms(fixed, moving, transformlist,
  else:
  return 1
  else:
- args = args + ['-z', 1, '--float', 1, '-e', imagetype, '-f', defaultvalue]
+ args = args + ['-z', str(1), '--float', str(int(singleprecision)), '-e', imagetype, '-f', defaultvalue]
  processed_args = utils._int_antsProcessArguments(args)
  libfn = utils.get_lib_fn('antsApplyTransforms')
  libfn(processed_args)

tests/test_core_ants_transform.py

@@ -87,7 +87,7 @@ def test_apply(self):
  img = ants.image_read(ants.get_ants_data("r16")).clone('float')
  tx = ants.new_ants_transform(dimension=2)
  tx.set_parameters((0.9,0,0,1.1,10,11))
- img2 = tx.apply(data=img, reference=img, data_type='image') 
+ img2 = tx.apply(data=img, reference=img, data_type='image')
 
  def test_apply_to_point(self):
  tx = ants.new_ants_transform()
@@ -99,14 +99,14 @@ def test_apply_to_vector(self):
  tx = ants.new_ants_transform()
  params = tx.parameters
  tx.set_parameters(params*2)
- pt2 = tx.apply_to_vector((1,2,3)) # should be (2,4,6) 
+ pt2 = tx.apply_to_vector((1,2,3)) # should be (2,4,6)
 
  def test_apply_to_image(self):
  for ptype in self.pixeltypes:
  img = ants.image_read(ants.get_ants_data("r16")).clone(ptype)
  tx = ants.new_ants_transform(dimension=2)
  tx.set_parameters((0.9,0,0,1.1,10,11))
- img2 = tx.apply_to_image(img, img) 
+ img2 = tx.apply_to_image(img, img)
 
 
 class TestModule_ants_transform(unittest.TestCase):
@@ -154,7 +154,8 @@ def test_apply_ants_transform(self):
  img = ants.image_read(ants.get_ants_data("r16")).clone('float')
  tx = ants.new_ants_transform(dimension=2)
  tx.set_parameters((0.9,0,0,1.1,10,11))
- img2 = ants.apply_ants_transform(tx, data=img, reference=img, data_type='image') 
+ img2 = ants.apply_ants_transform(tx, data=img, reference=img, data_type='image')
+
 
  def test_apply_ants_transform_to_point(self):
  tx = ants.new_ants_transform()

tests/test_registation.py

@@ -44,11 +44,21 @@ def test_example(self):
  fixed = ants.image_read(ants.get_ants_data("r16"))
  moving = ants.image_read(ants.get_ants_data("r64"))
  fixed = ants.resample_image(fixed, (64, 64), 1, 0)
- moving = ants.resample_image(moving, (64, 64), 1, 0)
+ moving = ants.resample_image(moving, (128, 128), 1, 0)
  mytx = ants.registration(fixed=fixed, moving=moving, type_of_transform="SyN")
  mywarpedimage = ants.apply_transforms(
  fixed=fixed, moving=moving, transformlist=mytx["fwdtransforms"]
  )
+ self.assertEqual(mywarpedimage.pixeltype, moving.pixeltype)
+ self.assertTrue(ants.ants_image.image_physical_space_consistency(fixed, mywarpedimage,
+ 0.0001, datatype = False))
+
+ # Call with float precision for transforms, but should still return input type
+ mywarpedimage2 = ants.apply_transforms(
+ fixed=fixed, moving=moving, transformlist=mytx["fwdtransforms"], singleprecision=True
+ )
+ self.assertEqual(mywarpedimage2.pixeltype, moving.pixeltype)
+ self.assertAlmostEqual(mywarpedimage.sum(), mywarpedimage2.sum(), places=3)
 
  # bad interpolator
  with self.assertRaises(Exception):