diff --git a/invesalius/data/imagedata_utils.py b/invesalius/data/imagedata_utils.py index 55e3a90..37e76cc 100644 --- a/invesalius/data/imagedata_utils.py +++ b/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 -- libgit2 0.21.2