diff --git a/invesalius/control.py b/invesalius/control.py index f437fa2..b4f3387 100644 --- a/invesalius/control.py +++ b/invesalius/control.py @@ -824,7 +824,6 @@ class Controller(): self.matrix, scalar_range, self.filename = image_utils.dcm2memmap(filelist, size, orientation, resolution_percentage) - print(xyspacing, zspacing) if orientation == 'AXIAL': spacing = xyspacing[0], xyspacing[1], zspacing elif orientation == 'CORONAL': @@ -832,7 +831,10 @@ class Controller(): elif orientation == 'SAGITTAL': spacing = zspacing, xyspacing[1], xyspacing[0] else: - self.matrix, spacing, scalar_range, self.filename = image_utils.dcmmf2memmap(filelist[0], orientation) + print(">>>>>> filelist", filelist) + self.matrix, scalar_range, spacing, self.filename = image_utils.dcmmf2memmap(filelist[0], orientation) + + print(">>>>>> spacing", spacing) self.Slice = sl.Slice() self.Slice.matrix = self.matrix diff --git a/invesalius/data/converters.py b/invesalius/data/converters.py index 3206ff9..418b4f2 100644 --- a/invesalius/data/converters.py +++ b/invesalius/data/converters.py @@ -17,10 +17,13 @@ # detalhes. #-------------------------------------------------------------------------- -import numpy +import gdcm +import numpy as np import vtk + from vtk.util import numpy_support + def to_vtk(n_array, spacing, slice_number, orientation, origin=(0, 0, 0), padding=(0, 0, 0)): if orientation == "SAGITTAL": @@ -84,3 +87,40 @@ def np_rgba_to_vtk(n_array, spacing=(1.0, 1.0, 1.0)): image.GetPointData().SetScalars(v_image) return image + + +# Based on http://gdcm.sourceforge.net/html/ConvertNumpy_8py-example.html +def gdcm_to_numpy(image, apply_intercep_scale=True): + map_gdcm_np = { + gdcm.PixelFormat.UINT8 :np.uint8, + gdcm.PixelFormat.INT8 :np.int8, + #gdcm.PixelFormat.UINT12 :np.uint12, + #gdcm.PixelFormat.INT12 :np.int12, + gdcm.PixelFormat.UINT16 :np.uint16, + gdcm.PixelFormat.INT16 :np.int16, + gdcm.PixelFormat.UINT32 :np.uint32, + gdcm.PixelFormat.INT32 :np.int32, + #gdcm.PixelFormat.FLOAT16:np.float16, + gdcm.PixelFormat.FLOAT32:np.float32, + gdcm.PixelFormat.FLOAT64:np.float64, + } + + pf = image.GetPixelFormat() + if image.GetNumberOfDimensions() == 3: + shape = image.GetDimension(2), image.GetDimension(1), image.GetDimension(0), pf.GetSamplesPerPixel() + else: + shape = image.GetDimension(1), image.GetDimension(0), pf.GetSamplesPerPixel() + dtype = map_gdcm_np[pf.GetScalarType()] + gdcm_array = image.GetBuffer() + np_array = np.frombuffer(gdcm_array.encode('utf-8', errors="surrogateescape"), dtype=dtype) + np_array.shape = shape + np_array = np_array.squeeze() + + if apply_intercep_scale: + shift = image.GetIntercept() + scale = image.GetSlope() + output = np.empty_like(np_array, np.int16) + output[:] = scale * np_array + shift + return output + else: + return np_array diff --git a/invesalius/data/imagedata_utils.py b/invesalius/data/imagedata_utils.py index 8d84b50..167d5e5 100644 --- a/invesalius/data/imagedata_utils.py +++ b/invesalius/data/imagedata_utils.py @@ -23,9 +23,10 @@ import sys import tempfile import gdcm +import imageio import numpy +import numpy as np import vtk -import vtkgdcm from wx.lib.pubsub import pub as Publisher from scipy.ndimage import shift, zoom @@ -124,20 +125,15 @@ def resize_slice(im_array, resolution_percentage): 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() +def read_dcm_slice_as_np2(filename, resolution_percentage=1.0): + reader = gdcm.ImageReader() + reader.SetFileName(filename) + reader.Read() + image = reader.GetImage() + output = converters.gdcm_to_numpy(image) 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 + output = zoom(output, resolution_percentage) + return output def FixGantryTilt(matrix, spacing, tilt): @@ -240,17 +236,6 @@ def View(imagedata): import time time.sleep(10) -def ViewGDCM(imagedata): - viewer = vtkgdcm.vtkImageColorViewer() - viewer.SetInput(reader.GetOutput()) - viewer.SetColorWindow(500.) - viewer.SetColorLevel(50.) - viewer.Render() - - import time - time.sleep(5) - - def ExtractVOI(imagedata,xi,xf,yi,yf,zi,zf): """ @@ -265,267 +250,33 @@ def ExtractVOI(imagedata,xi,xf,yi,yf,zi,zf): return voi.GetOutput() -def create_dicom_thumbnails(filename, window=None, level=None): - rvtk = vtkgdcm.vtkGDCMImageReader() - rvtk.SetFileName(utils.encode(filename, const.FS_ENCODE)) - rvtk.Update() - - img = rvtk.GetOutput() +def create_dicom_thumbnails(image, window=None, level=None): + pf = image.GetPixelFormat() + np_image = converters.gdcm_to_numpy(image, pf.GetSamplesPerPixel() == 1) if window is None or level is None: - _min, _max = img.GetScalarRange() + _min, _max = np_image.min(), np_image.max() window = _max - _min level = _min + window / 2 - dx, dy, dz = img.GetDimensions() - - if dz > 1: + if image.GetNumberOfDimensions() >= 3: thumbnail_paths = [] - for i in range(dz): - img_slice = ExtractVOI(img, 0, dx-1, 0, dy-1, i, i+1) - - colorer = vtk.vtkImageMapToWindowLevelColors() - colorer.SetInputData(img_slice) - colorer.SetWindow(window) - colorer.SetLevel(level) - colorer.SetOutputFormatToRGB() - colorer.Update() - - resample = vtk.vtkImageResample() - resample.SetInputData(colorer.GetOutput()) - resample.SetAxisMagnificationFactor ( 0, 0.25 ) - resample.SetAxisMagnificationFactor ( 1, 0.25 ) - resample.SetAxisMagnificationFactor ( 2, 1 ) - resample.Update() - - thumbnail_path = tempfile.mktemp() - - write_png = vtk.vtkPNGWriter() - write_png.SetInputData(resample.GetOutput()) - write_png.SetFileName(thumbnail_path) - write_png.Write() - + for i in range(np_image.shape[0]): + thumb_image = zoom(np_image[i], 0.25) + thumb_image = np.array(get_LUT_value_255(thumb_image, window, level), dtype=np.uint8) + thumbnail_path = tempfile.mktemp(prefix='thumb_', suffix='.png') + imageio.imsave(thumbnail_path, thumb_image) thumbnail_paths.append(thumbnail_path) - return thumbnail_paths else: - colorer = vtk.vtkImageMapToWindowLevelColors() - colorer.SetInputData(img) - colorer.SetWindow(window) - colorer.SetLevel(level) - colorer.SetOutputFormatToRGB() - colorer.Update() - - resample = vtk.vtkImageResample() - resample.SetInputData(colorer.GetOutput()) - resample.SetAxisMagnificationFactor ( 0, 0.25 ) - resample.SetAxisMagnificationFactor ( 1, 0.25 ) - resample.SetAxisMagnificationFactor ( 2, 1 ) - resample.Update() - - thumbnail_path = tempfile.mktemp() - - write_png = vtk.vtkPNGWriter() - write_png.SetInputData(resample.GetOutput()) - write_png.SetFileName(thumbnail_path) - write_png.Write() - - return thumbnail_path - - -def CreateImageData(filelist, zspacing, xyspacing,size, - bits, use_dcmspacing): - message = _("Generating multiplanar visualization...") - - if not const.VTK_WARNING: - log_path = os.path.join(inv_paths.USER_LOG_DIR, 'vtkoutput.txt') - fow = vtk.vtkFileOutputWindow() - fow.SetFileName(log_path) - ow = vtk.vtkOutputWindow() - ow.SetInstance(fow) - - x,y = size - px, py = utils.predict_memory(len(filelist), x, y, bits) - - utils.debug("Image Resized to >>> %f x %f" % (px, py)) - - if (x == px) and (y == py): - const.REDUCE_IMAGEDATA_QUALITY = 0 - else: - const.REDUCE_IMAGEDATA_QUALITY = 1 - - if not(const.REDUCE_IMAGEDATA_QUALITY): - update_progress= vtk_utils.ShowProgress(1, dialog_type = "ProgressDialog") - - array = vtk.vtkStringArray() - for x in range(len(filelist)): - array.InsertValue(x,filelist[x]) - - reader = vtkgdcm.vtkGDCMImageReader() - reader.SetFileNames(array) - reader.AddObserver("ProgressEvent", lambda obj,evt: - update_progress(reader,message)) - reader.Update() - - # The zpacing is a DicomGroup property, so we need to set it - imagedata = vtk.vtkImageData() - imagedata.DeepCopy(reader.GetOutput()) - if (use_dcmspacing): - spacing = xyspacing - spacing[2] = zspacing - else: - spacing = imagedata.GetSpacing() - - imagedata.SetSpacing(spacing[0], spacing[1], zspacing) - else: - - update_progress= vtk_utils.ShowProgress(2*len(filelist), - dialog_type = "ProgressDialog") - - # Reformat each slice and future append them - appender = vtk.vtkImageAppend() - appender.SetAppendAxis(2) #Define Stack in Z - - - # Reformat each slice - for x in range(len(filelist)): - # TODO: We need to check this automatically according - # to each computer's architecture - # If the resolution of the matrix is too large - reader = vtkgdcm.vtkGDCMImageReader() - reader.SetFileName(filelist[x]) - reader.AddObserver("ProgressEvent", lambda obj,evt: - update_progress(reader,message)) - reader.Update() - - if (use_dcmspacing): - spacing = xyspacing - spacing[2] = zspacing - else: - spacing = reader.GetOutput().GetSpacing() - - tmp_image = vtk.vtkImageData() - tmp_image.DeepCopy(reader.GetOutput()) - tmp_image.SetSpacing(spacing[0], spacing[1], zspacing) - tmp_image.Update() - - #Resample image in x,y dimension - slice_imagedata = ResampleImage2D(tmp_image, px, py, update_progress) - #Stack images in Z axes - appender.AddInput(slice_imagedata) - #appender.AddObserver("ProgressEvent", lambda obj,evt:update_progress(appender)) - appender.Update() - - spacing = appender.GetOutput().GetSpacing() - - # The zpacing is a DicomGroup property, so we need to set it - imagedata = vtk.vtkImageData() - imagedata.DeepCopy(appender.GetOutput()) - imagedata.SetSpacing(spacing[0], spacing[1], zspacing) - - imagedata.AddObserver("ProgressEvent", lambda obj,evt: - update_progress(imagedata,message)) - imagedata.Update() - - return imagedata - - -class ImageCreator: - def __init__(self): - self.running = True - Publisher.subscribe(self.CancelImageDataLoad, "Cancel DICOM load") - - def CancelImageDataLoad(self, evt_pusub): - utils.debug("Canceling") - self.running = False - - def CreateImageData(self, filelist, zspacing, size, bits): - message = _("Generating multiplanar visualization...") - - if not const.VTK_WARNING: - log_path = os.path.join(inv_paths.USER_LOG_DIR, 'vtkoutput.txt') - fow = vtk.vtkFileOutputWindow() - fow.SetFileName(log_path) - ow = vtk.vtkOutputWindow() - ow.SetInstance(fow) - - x,y = size - px, py = utils.predict_memory(len(filelist), x, y, bits) - utils.debug("Image Resized to >>> %f x %f" % (px, py)) - - if (x == px) and (y == py): - const.REDUCE_IMAGEDATA_QUALITY = 0 + thumbnail_path = tempfile.mktemp(prefix='thumb_', suffix='.png') + if pf.GetSamplesPerPixel() == 1: + thumb_image = zoom(np_image, 0.25) + thumb_image = np.array(get_LUT_value_255(thumb_image, window, level), dtype=np.uint8) else: - const.REDUCE_IMAGEDATA_QUALITY = 1 - - if not(const.REDUCE_IMAGEDATA_QUALITY): - update_progress= vtk_utils.ShowProgress(1, dialog_type = "ProgressDialog") - - array = vtk.vtkStringArray() - for x in range(len(filelist)): - if not self.running: - return False - array.InsertValue(x,filelist[x]) - - if not self.running: - return False - reader = vtkgdcm.vtkGDCMImageReader() - reader.SetFileNames(array) - reader.AddObserver("ProgressEvent", lambda obj,evt: - update_progress(reader,message)) - reader.Update() - - if not self.running: - reader.AbortExecuteOn() - return False - # The zpacing is a DicomGroup property, so we need to set it - imagedata = vtk.vtkImageData() - imagedata.DeepCopy(reader.GetOutput()) - spacing = imagedata.GetSpacing() - imagedata.SetSpacing(spacing[0], spacing[1], zspacing) - else: - - update_progress= vtk_utils.ShowProgress(2*len(filelist), - dialog_type = "ProgressDialog") - - # Reformat each slice and future append them - appender = vtk.vtkImageAppend() - appender.SetAppendAxis(2) #Define Stack in Z - - - # Reformat each slice - for x in range(len(filelist)): - # TODO: We need to check this automatically according - # to each computer's architecture - # If the resolution of the matrix is too large - if not self.running: - return False - reader = vtkgdcm.vtkGDCMImageReader() - reader.SetFileName(filelist[x]) - reader.AddObserver("ProgressEvent", lambda obj,evt: - update_progress(reader,message)) - reader.Update() - - #Resample image in x,y dimension - slice_imagedata = ResampleImage2D(reader.GetOutput(), px, py, update_progress) - #Stack images in Z axes - appender.AddInput(slice_imagedata) - #appender.AddObserver("ProgressEvent", lambda obj,evt:update_progress(appender)) - appender.Update() - - # The zpacing is a DicomGroup property, so we need to set it - if not self.running: - return False - imagedata = vtk.vtkImageData() - imagedata.DeepCopy(appender.GetOutput()) - spacing = imagedata.GetSpacing() - - imagedata.SetSpacing(spacing[0], spacing[1], zspacing) - - imagedata.AddObserver("ProgressEvent", lambda obj,evt: - update_progress(imagedata,message)) - imagedata.Update() + thumb_image = zoom(np_image, (0.25, 0.25, 1)) + imageio.imsave(thumbnail_path, thumb_image) + return thumbnail_path - return imagedata def bitmap2memmap(files, slice_size, orientation, spacing, resolution_percentage): """ @@ -626,7 +377,6 @@ def bitmap2memmap(files, slice_size, orientation, spacing, resolution_percentage return matrix, scalar_range, temp_file - def dcm2memmap(files, slice_size, orientation, resolution_percentage): """ From a list of dicom files it creates memmap file in the temp folder and @@ -635,7 +385,7 @@ 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) + first_slice = read_dcm_slice_as_np2(files[0], resolution_percentage) slice_size = first_slice.shape[::-1] temp_file = tempfile.mktemp() @@ -649,7 +399,7 @@ def dcm2memmap(files, slice_size, orientation, resolution_percentage): matrix = numpy.memmap(temp_file, mode='w+', dtype='int16', shape=shape) for n, f in enumerate(files): - im_array = read_dcm_slice_as_np(f, resolution_percentage) + im_array = read_dcm_slice_as_np2(f, resolution_percentage)[::-1] if orientation == 'CORONAL': matrix[:, shape[1] - n - 1, :] = im_array @@ -669,37 +419,35 @@ def dcm2memmap(files, slice_size, orientation, resolution_percentage): def dcmmf2memmap(dcm_file, orientation): - r = vtkgdcm.vtkGDCMImageReader() - r.SetFileName(dcm_file) - r.Update() - + reader = gdcm.ImageReader() + reader.SetFileName(dcm_file) + reader.Read() + image = reader.GetImage() + xs, ys, zs = image.GetSpacing() + pf = image.GetPixelFormat() + np_image = converters.gdcm_to_numpy(image, pf.GetSamplesPerPixel() == 1) temp_file = tempfile.mktemp() - - o = r.GetOutput() - x, y, z = o.GetDimensions() - spacing = o.GetSpacing() - - matrix = numpy.memmap(temp_file, mode='w+', dtype='int16', shape=(z, y, x)) - - d = numpy_support.vtk_to_numpy(o.GetPointData().GetScalars()) - d.shape = z, y, x + matrix = numpy.memmap(temp_file, mode='w+', dtype='int16', shape=np_image.shape) + print("Number of dimensions", np_image.shape) + z, y, x = np_image.shape if orientation == 'CORONAL': + spacing = xs, zs, ys matrix.shape = y, z, x for n in range(z): - matrix[:, n, :] = d[n] + matrix[:, n, :] = np_image[n][::-1] elif orientation == 'SAGITTAL': - matrix.shape = x, z, y + spacing = zs, ys, xs + matrix.shape = y, x, z for n in range(z): - matrix[:, :, n] = d[n] + matrix[:, :, n] = np_image[n][::-1] else: - matrix[:] = d + spacing = xs, ys, zs + matrix[:] = np_image[:, ::-1, :] matrix.flush() scalar_range = matrix.min(), matrix.max() - print("ORIENTATION", orientation) - - return matrix, spacing, scalar_range, temp_file + return matrix, scalar_range, spacing, temp_file def img2memmap(group): @@ -753,3 +501,14 @@ def imgnormalize(data, srange=(0, 255)): datan = datan.astype(numpy.int16) return datan + + +def get_LUT_value_255(data, window, level): + shape = data.shape + data_ = data.ravel() + data = np.piecewise(data_, + [data_ <= (level - 0.5 - (window-1)/2), + data_ > (level - 0.5 + (window-1)/2)], + [0, 255, lambda data_: ((data_ - (level - 0.5))/(window-1) + 0.5)*(255)]) + data.shape = shape + return data diff --git a/invesalius/gui/bitmap_preview_panel.py b/invesalius/gui/bitmap_preview_panel.py index 12b8968..6564bb0 100644 --- a/invesalius/gui/bitmap_preview_panel.py +++ b/invesalius/gui/bitmap_preview_panel.py @@ -1,6 +1,5 @@ import wx import vtk -import vtkgdcm import time import numpy diff --git a/invesalius/gui/dicom_preview_panel.py b/invesalius/gui/dicom_preview_panel.py index 8d66d67..ffb8806 100644 --- a/invesalius/gui/dicom_preview_panel.py +++ b/invesalius/gui/dicom_preview_panel.py @@ -35,7 +35,8 @@ import invesalius.constants as const import invesalius.reader.dicom_reader as dicom_reader import invesalius.data.vtk_utils as vtku import invesalius.utils as utils -import vtkgdcm +from invesalius.data import imagedata_utils +from invesalius.data import converters if sys.platform == 'win32': try: @@ -890,20 +891,20 @@ class SingleImagePreview(wx.Panel): reader.Update() image = reader.GetOutput() - else: - rdicom = vtkgdcm.vtkGDCMImageReader() + filename = dicom.image.file if _has_win32api: - rdicom.SetFileName(win32api.GetShortPathName(dicom.image.file).encode(const.FS_ENCODE)) - else: - rdicom.SetFileName(dicom.image.file) - rdicom.Update() + filename = win32api.GetShortPathName(filename).encode(const.FS_ENCODE) + + np_image = imagedata_utils.read_dcm_slice_as_np2(filename) + print(">>> spacing", dicom.image.spacing) + vtk_image = converters.to_vtk(np_image, dicom.image.spacing, 0, 'AXIAL') # ADJUST CONTRAST window_level = dicom.image.level window_width = dicom.image.window colorer = vtk.vtkImageMapToWindowLevelColors() - colorer.SetInputConnection(rdicom.GetOutputPort()) + colorer.SetInputData(vtk_image) colorer.SetWindow(float(window_width)) colorer.SetLevel(float(window_level)) colorer.Update() diff --git a/invesalius/reader/dicom.py b/invesalius/reader/dicom.py index 41d570f..72586e7 100644 --- a/invesalius/reader/dicom.py +++ b/invesalius/reader/dicom.py @@ -19,7 +19,6 @@ #--------------------------------------------------------------------- import time #import gdcm -#import vtkgdcm import sys import invesalius.utils as utils import invesalius.constants as const diff --git a/invesalius/reader/dicom_reader.py b/invesalius/reader/dicom_reader.py index 9208c69..d4ada61 100644 --- a/invesalius/reader/dicom_reader.py +++ b/invesalius/reader/dicom_reader.py @@ -24,7 +24,6 @@ import sys from multiprocessing import cpu_count import vtk -import vtkgdcm import gdcm from wx.lib.pubsub import pub as Publisher @@ -202,14 +201,11 @@ class LoadDicom: level = None window = None - if _has_win32api: - thumbnail_path = imagedata_utils.create_dicom_thumbnails(win32api.GetShortPathName(self.filepath), window, level) - else: - thumbnail_path = imagedata_utils.create_dicom_thumbnails(self.filepath, window, level) + img = reader.GetImage() + thumbnail_path = imagedata_utils.create_dicom_thumbnails(img, window, level) #------ Verify the orientation -------------------------------- - img = reader.GetImage() direc_cosines = img.GetDirectionCosines() orientation = gdcm.Orientation() try: -- libgit2 0.21.2