Resampling: Ghost slices due to rounding errors

dicom2nifti/convert_generic.py

@@ -435,8 +435,10 @@ def _convert_slice_incement_inconsistencies(dicom_input):
         if data.ndim > 3:  # do not squeeze single slice data
             data = data.squeeze()
         current_volume = nibabel.Nifti1Image(data, affine)
-        slice_increment = numpy.linalg.norm(current_volume.header.get_zooms())
-        voxel_sizes['%.5f' % slice_increment] = current_volume.header.get_zooms()
+        # Compute voxel size in float64, because header.get_zooms() is rounded to float32
+        slice_voxel_size = numpy.linalg.norm(affine[:3, :3], axis=0)
+        slice_increment = numpy.linalg.norm(slice_voxel_size)
+        voxel_sizes['%.5f' % slice_increment] = slice_voxel_size
         slice_increments.extend([slice_increment] * (len(dicom_slices) - 1))
         slice_incement_niftis.append(current_volume)
 

dicom2nifti/resample.py

@@ -50,9 +50,10 @@ def resample_nifti_images(nifti_images, voxel_size=None):
 
     # get the smallest voxelsize and use that
     if voxel_size is None:
-        voxel_size = nifti_images[0].header.get_zooms()
-        for nifti_image in nifti_images[1:]:
-            voxel_size = numpy.minimum(voxel_size, nifti_image.header.get_zooms())
+        # Compute voxel size from float64 affine, because header.get_zooms() is rounded to float32
+        voxel_sizes = [numpy.linalg.norm(nifti_image.affine[:3, :3], axis=0)                       
+                       for nifti_image in nifti_images]
+        voxel_size = numpy.array(voxel_sizes).min(axis=0)
 
     x_axis_world = numpy.transpose(numpy.dot(nifti_images[0].affine, [[1], [0], [0], [0]]))[0, :3]
     y_axis_world = numpy.transpose(numpy.dot(nifti_images[0].affine, [[0], [1], [0], [0]]))[0, :3]
@@ -99,7 +100,10 @@ def resample_nifti_images(nifti_images, voxel_size=None):
              min_projected[2] * z_axis_world
 
     new_voxelsize = voxel_size
-    new_shape = numpy.ceil(new_size_mm / new_voxelsize).astype(numpy.int16) + 1
+    new_shape_float = new_size_mm / new_voxelsize + 1
+    new_shape = numpy.ceil(
+        numpy.round(new_shape_float, decimals=2)  # Floor values below 0.005
+    ).astype(numpy.int16)
 
     new_affine = _create_affine(x_axis_world, y_axis_world, z_axis_world, origin, voxel_size)
 

tests/test_resample.py

@@ -0,0 +1,23 @@
+import numpy
+from nibabel import Nifti1Image
+
+from dicom2nifti.resample import resample_nifti_images
+
+
+def test_resample():
+    z_spacing = 0.45  # In float64: 0.45000000000000001
+    nifti = Nifti1Image(
+        dataobj=numpy.ones((3, 3, 3)),
+        affine=numpy.array(
+            [
+                [1, 0, 0, 0],
+                [0, 1, 0, 0],
+                [0, 0, z_spacing, 0],
+                [0, 0, 0, 1],
+            ]
+        ),
+    )
+    voxel_size = nifti.header.get_zooms()  # Rounds to float32: 0.44999998807907104492
+    nifti_out = resample_nifti_images([nifti], voxel_size)
+    # shape was (3, 3, 4) before rounding fix in resample_nifti_images()
+    assert nifti_out.shape == (3, 3, 3)