Using the resample method to calculate the final image size when rescaling (closes #126)

Thiago Franco de Moraes
1 parent 6d9e7e2a
Showing 1 changed file with 54 additions and 60 deletions Show diff stats
invesalius/data/imagedata_utils.py
@@ -27,7 +27,7 @@ import vtk
 import vtkgdcm
 from wx.lib.pubsub import pub as Publisher
-from scipy.ndimage import shift
+from scipy.ndimage import shift, zoom
 from vtk.util import numpy_support
 import invesalius.constants as const
@@ -69,20 +69,20 @@ def ResampleImage2D(imagedata, px=None, py=None, resolution_percentage = None,
     dimensions = imagedata.GetDimensions()
     if resolution_percentage:
-        px = math.ceil(dimensions[0] * resolution_percentage)
-        py = math.ceil(dimensions[1] * resolution_percentage)
-
-    if abs(extent[1]-extent[3]) < abs(extent[3]-extent[5]):
-        f = extent[1]
-    elif abs(extent[1]-extent[5]) < abs(extent[1] - extent[3]):
-        f = extent[1]
-    elif abs(extent[3]-extent[5]) < abs(extent[1] - extent[3]):
-        f = extent[3]
+        factor_x = resolution_percentage
+        factor_y = resolution_percentage
     else:
-        f = extent[1]
+        if abs(extent[1]-extent[3]) < abs(extent[3]-extent[5]):
+            f = extent[1]
+        elif abs(extent[1]-extent[5]) < abs(extent[1] - extent[3]):
+            f = extent[1]
+        elif abs(extent[3]-extent[5]) < abs(extent[1] - extent[3]):
+            f = extent[3]
+        else:
+            f = extent[1]
-    factor_x = px/float(f+1)
-    factor_y = py/float(f+1)
+        factor_x = px/float(f+1)
+        factor_y = py/float(f+1)
     resample = vtk.vtkImageResample()
     resample.SetInputData(imagedata)
@@ -98,6 +98,35 @@ def ResampleImage2D(imagedata, px=None, py=None, resolution_percentage = None,
     return resample.GetOutput()
+
+def resize_slice(im_array, resolution_percentage):
+    """
+    Uses ndimage.zoom to resize a slice.
+
+    input:
+        im_array: slice as a numpy array.
+        resolution_percentage: percentage of resize.
+    """
+    out = zoom(im_array, resolution_percentage, im_array.dtype, order=2)
+    return out
+
+
+def read_dcm_slice_as_np(filename, resolution_percentage=1.0):
+    """
+    read a dicom slice file and return the slice as numpy ndarray
+    """
+    dcm_reader = vtkgdcm.vtkGDCMImageReader()
+    dcm_reader.SetFileName(filename)
+    dcm_reader.Update()
+    image = dcm_reader.GetOutput()
+    if resolution_percentage < 1.0:
+        image = ResampleImage2D(image, resolution_percentage=resolution_percentage)
+    dx, dy, dz = image.GetDimensions()
+    im_array = numpy_support.vtk_to_numpy(image.GetPointData().GetScalars())
+    im_array.shape = dy, dx
+    return im_array
+
+
 def FixGantryTilt(matrix, spacing, tilt):
     """
     Fix gantry tilt given a vtkImageData and the tilt value. Return new
@@ -525,70 +554,35 @@ def dcm2memmap(files, slice_size, orientation, resolution_percentage):
     message = _("Generating multiplanar visualization...")
     update_progress= vtk_utils.ShowProgress(len(files) - 1, dialog_type = "ProgressDialog")
+    first_slice = read_dcm_slice_as_np(files[0], resolution_percentage)
+    slice_size = first_slice.shape[::-1]
+
     temp_file = tempfile.mktemp()
     if orientation == 'SAGITTAL':
-        if resolution_percentage == 1.0:
-            shape = slice_size[0], slice_size[1], len(files)
-        else:
-            shape = int(math.ceil(slice_size[0]*resolution_percentage)),\
-                    int(math.ceil(slice_size[1]*resolution_percentage)), len(files)
-
+        shape = slice_size[0], slice_size[1], len(files)
     elif orientation == 'CORONAL':
-        if resolution_percentage == 1.0:
-            shape = slice_size[1], len(files), slice_size[0]
-        else:
-            shape = int(math.ceil(slice_size[1]*resolution_percentage)), len(files),\
-                                        int(math.ceil(slice_size[0]*resolution_percentage))
+        shape = slice_size[1], len(files), slice_size[0]
     else:
-        if resolution_percentage == 1.0:
-            shape = len(files), slice_size[1], slice_size[0]
-        else:
-            shape = len(files), int(math.ceil(slice_size[1]*resolution_percentage)),\
-                                int(math.ceil(slice_size[0]*resolution_percentage))
+        shape = len(files), slice_size[1], slice_size[0]
     matrix = numpy.memmap(temp_file, mode='w+', dtype='int16', shape=shape)
-    dcm_reader = vtkgdcm.vtkGDCMImageReader()
-    cont = 0
-    max_scalar = None
-    min_scalar = None
-
     for n, f in enumerate(files):
-        dcm_reader.SetFileName(f)
-        dcm_reader.Update()
-        image = dcm_reader.GetOutput()
-
-        if resolution_percentage != 1.0:
-            image_resized = ResampleImage2D(image, px=None, py=None,\
-                                resolution_percentage = resolution_percentage, update_progress = None)
-
-            image = image_resized
-
-        min_aux, max_aux = image.GetScalarRange()
-        if min_scalar is None or min_aux < min_scalar:
-            min_scalar = min_aux
-
-        if max_scalar is None or max_aux > max_scalar:
-            max_scalar = max_aux
+        im_array = read_dcm_slice_as_np(f, resolution_percentage)
-        array = numpy_support.vtk_to_numpy(image.GetPointData().GetScalars())
         if orientation == 'CORONAL':
-            array.shape = matrix.shape[0], matrix.shape[2]
-            matrix[:, shape[1] - n - 1, :] = array
+            matrix[:, shape[1] - n - 1, :] = im_array
         elif orientation == 'SAGITTAL':
-            array.shape = matrix.shape[0], matrix.shape[1]
             # TODO: Verify if it's necessary to add the slices swapped only in
             # sagittal rmi or only in # Rasiane's case or is necessary in all
             # sagittal cases.
-            matrix[:, :, n] = array
+            matrix[:, :, n] = im_array
         else:
-            array.shape = matrix.shape[1], matrix.shape[2]
-            matrix[n] = array
-        update_progress(cont,message)
-        cont += 1
+            matrix[n] = im_array
+        update_progress(n, message)
     matrix.flush()
-    scalar_range = min_scalar, max_scalar
+    scalar_range = matrix.min(), matrix.max()
     return matrix, scalar_range, temp_file