From c60264663d8996611a3f055a5de54c86e2a3ea05 Mon Sep 17 00:00:00 2001 From: Thiago Franco de Moraes Date: Thu, 27 Sep 2018 15:10:36 -0300 Subject: [PATCH] Cancelable surface creation (#147) closes #142 --- invesalius/data/interpolation.pyx | 4 ++-- invesalius/data/surface.py | 537 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- invesalius/data/surface_process.py | 501 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------------------------------------------- invesalius/gui/dialogs.py | 27 +++++++++++++++++++++++++++ invesalius/utils.py | 23 +++++++++++++++++++++++ 5 files changed, 603 insertions(+), 489 deletions(-) diff --git a/invesalius/data/interpolation.pyx b/invesalius/data/interpolation.pyx index 8aaea8d..9f44648 100644 --- a/invesalius/data/interpolation.pyx +++ b/invesalius/data/interpolation.pyx @@ -4,7 +4,7 @@ import numpy as np cimport numpy as np cimport cython -from libc.math cimport floor, ceil, sqrt, fabs, round, sin, M_PI +from libc.math cimport floor, ceil, sqrt, fabs, sin, M_PI from cython.parallel import prange DEF LANCZOS_A = 4 @@ -81,7 +81,7 @@ cdef double[64][64] temp = [ @cython.cdivision(True) @cython.wraparound(False) cdef double nearest_neighbour_interp(image_t[:, :, :] V, double x, double y, double z) nogil: - return V[round(z), round(y), round(x)] + return V[(z), (y), (x)] @cython.boundscheck(False) # turn of bounds-checking for entire function @cython.cdivision(True) diff --git a/invesalius/data/surface.py b/invesalius/data/surface.py index fbc41ae..56ee6d6 100644 --- a/invesalius/data/surface.py +++ b/invesalius/data/surface.py @@ -17,6 +17,7 @@ # detalhes. #-------------------------------------------------------------------------- +import functools import multiprocessing import os import plistlib @@ -24,10 +25,19 @@ import random import shutil import sys import tempfile +import time +import traceback import weakref +try: + import queue +except ImportError: + import Queue as queue + import vtk import wx +import wx.lib.agw.genericmessagedialog as GMD + from wx.lib.pubsub import pub as Publisher if sys.platform == 'win32': @@ -48,9 +58,12 @@ import invesalius.data.surface_process as surface_process import invesalius.utils as utl import invesalius.data.vtk_utils as vu +from invesalius.gui import dialogs + from invesalius.data import cy_mesh # TODO: Verificar ReleaseDataFlagOn and SetSource + class Surface(): """ Represent both vtkPolyData and associated properties. @@ -438,10 +451,81 @@ class SurfaceManager(): #### #(mask_index, surface_name, quality, fill_holes, keep_largest) + def _on_complete_surface_creation(self, args, overwrite, surface_name, colour, dialog): + surface_filename, surface_measures = args + print(surface_filename, surface_measures) + reader = vtk.vtkXMLPolyDataReader() + reader.SetFileName(surface_filename) + reader.Update() + polydata = reader.GetOutput() + + # Map polygonal data (vtkPolyData) to graphics primitives. + mapper = vtk.vtkPolyDataMapper() + mapper.SetInputData(polydata) + mapper.ScalarVisibilityOff() + # mapper.ReleaseDataFlagOn() + mapper.ImmediateModeRenderingOn() # improve performance + + # Represent an object (geometry & properties) in the rendered scene + actor = vtk.vtkActor() + actor.SetMapper(mapper) + del mapper + #Create Surface instance + if overwrite: + surface = Surface(index = self.last_surface_index) + else: + surface = Surface(name=surface_name) + surface.colour = colour + surface.polydata = polydata + surface.volume = surface_measures['volume'] + surface.area = surface_measures['area'] + del polydata + + # Set actor colour and transparency + actor.GetProperty().SetColor(colour) + actor.GetProperty().SetOpacity(1-surface.transparency) + + prop = actor.GetProperty() + + interpolation = int(ses.Session().surface_interpolation) + + prop.SetInterpolation(interpolation) + + proj = prj.Project() + if overwrite: + proj.ChangeSurface(surface) + else: + index = proj.AddSurface(surface) + surface.index = index + self.last_surface_index = index + + session = ses.Session() + session.ChangeProject() + + Publisher.sendMessage('Load surface actor into viewer', actor=actor) + + # Send actor by pubsub to viewer's render + if overwrite and self.actors_dict.keys(): + old_actor = self.actors_dict[self.last_surface_index] + Publisher.sendMessage('Remove surface actor from viewer', actor=old_actor) + + # Save actor for future management tasks + self.actors_dict[surface.index] = actor + Publisher.sendMessage('Update surface info in GUI', surface=surface) + Publisher.sendMessage('End busy cursor') + + dialog.running = False + + def _on_callback_error(self, e, dialog=None): + dialog.running = False + msg = utl.log_traceback(e) + dialog.error = msg + def AddNewActor(self, slice_, mask, surface_parameters): """ Create surface actor, save into project and send it to viewer. """ + t_init = time.time() matrix = slice_.matrix filename_img = slice_.matrix_filename spacing = slice_.spacing @@ -470,9 +554,6 @@ class SurfaceManager(): smooth_relaxation_factor = const.SURFACE_QUALITY[quality][2] decimate_reduction = const.SURFACE_QUALITY[quality][3] - #if imagedata_resolution: - #imagedata = iu.ResampleImage3D(imagedata, imagedata_resolution) - pipeline_size = 4 if decimate_reduction: pipeline_size += 1 @@ -483,10 +564,6 @@ class SurfaceManager(): if keep_largest: pipeline_size += 1 - ## Update progress value in GUI - UpdateProgress = vu.ShowProgress(pipeline_size) - UpdateProgress(0, _("Creating 3D surface...")) - language = ses.Session().language if (prj.Project().original_orientation == const.CORONAL): @@ -496,220 +573,59 @@ class SurfaceManager(): n_processors = multiprocessing.cpu_count() - pipe_in, pipe_out = multiprocessing.Pipe() o_piece = 1 - piece_size = 2000 + piece_size = 20 n_pieces = int(round(matrix.shape[0] / piece_size + 0.5, 0)) - q_in = multiprocessing.Queue() - q_out = multiprocessing.Queue() - - p = [] - for i in range(n_processors): - sp = surface_process.SurfaceProcess(pipe_in, filename_img, - matrix.shape, matrix.dtype, - mask.temp_file, - mask.matrix.shape, - mask.matrix.dtype, - spacing, - mode, min_value, max_value, - decimate_reduction, - smooth_relaxation_factor, - smooth_iterations, language, - flip_image, q_in, q_out, - algorithm != 'Default', - algorithm, - imagedata_resolution) - p.append(sp) - sp.start() - - for i in range(n_pieces): - init = i * piece_size - end = init + piece_size + o_piece - roi = slice(init, end) - q_in.put(roi) - print("new_piece", roi) - - for i in p: - q_in.put(None) - - none_count = 1 - while 1: - msg = pipe_out.recv() - if(msg is None): - none_count += 1 - else: - UpdateProgress(msg[0]/(n_pieces * pipeline_size), msg[1]) - - if none_count > n_pieces: - break + filenames = [] + pool = multiprocessing.Pool(processes=min(n_pieces, n_processors)) + manager = multiprocessing.Manager() + msg_queue = manager.Queue(1) - polydata_append = vtk.vtkAppendPolyData() - # polydata_append.ReleaseDataFlagOn() - t = n_pieces - while t: - filename_polydata = q_out.get() + # If InVesalius is running without GUI + if wx.GetApp() is None: + for i in range(n_pieces): + init = i * piece_size + end = init + piece_size + o_piece + roi = slice(init, end) + print("new_piece", roi) + f = pool.apply_async(surface_process.create_surface_piece, + args = (filename_img, matrix.shape, matrix.dtype, + mask.temp_file, mask.matrix.shape, + mask.matrix.dtype, roi, spacing, mode, + min_value, max_value, decimate_reduction, + smooth_relaxation_factor, + smooth_iterations, language, flip_image, + algorithm != 'Default', algorithm, + imagedata_resolution), + callback=lambda x: filenames.append(x)) + + while len(filenames) != n_pieces: + time.sleep(0.25) + + f = pool.apply_async(surface_process.join_process_surface, + args=(filenames, algorithm, smooth_iterations, + smooth_relaxation_factor, + decimate_reduction, keep_largest, + fill_holes, options, msg_queue)) + + while not f.ready(): + time.sleep(0.25) + + try: + surface_filename, surface_measures = f.get() + except Exception as e: + print(_("InVesalius was not able to create the surface")) + print(traceback.print_exc()) + return reader = vtk.vtkXMLPolyDataReader() - reader.SetFileName(filename_polydata) - # reader.ReleaseDataFlagOn() + reader.SetFileName(surface_filename) reader.Update() - # reader.GetOutput().ReleaseDataFlagOn() polydata = reader.GetOutput() - # polydata.SetSource(None) - - polydata_append.AddInputData(polydata) - del reader - del polydata - t -= 1 - - polydata_append.Update() - # polydata_append.GetOutput().ReleaseDataFlagOn() - polydata = polydata_append.GetOutput() - #polydata.Register(None) - # polydata.SetSource(None) - del polydata_append - - if algorithm == 'ca_smoothing': - normals = vtk.vtkPolyDataNormals() - normals_ref = weakref.ref(normals) - normals_ref().AddObserver("ProgressEvent", lambda obj,evt: - UpdateProgress(normals_ref(), _("Creating 3D surface..."))) - normals.SetInputData(polydata) - # normals.ReleaseDataFlagOn() - #normals.SetFeatureAngle(80) - #normals.AutoOrientNormalsOn() - normals.ComputeCellNormalsOn() - # normals.GetOutput().ReleaseDataFlagOn() - normals.Update() - del polydata - polydata = normals.GetOutput() - # polydata.SetSource(None) - del normals - - clean = vtk.vtkCleanPolyData() - # clean.ReleaseDataFlagOn() - # clean.GetOutput().ReleaseDataFlagOn() - clean_ref = weakref.ref(clean) - clean_ref().AddObserver("ProgressEvent", lambda obj,evt: - UpdateProgress(clean_ref(), _("Creating 3D surface..."))) - clean.SetInputData(polydata) - clean.PointMergingOn() - clean.Update() - - del polydata - polydata = clean.GetOutput() - # polydata.SetSource(None) - del clean - - # try: - # polydata.BuildLinks() - # except TypeError: - # polydata.BuildLinks(0) - # polydata = ca_smoothing.ca_smoothing(polydata, options['angle'], - # options['max distance'], - # options['min weight'], - # options['steps']) - - mesh = cy_mesh.Mesh(polydata) - cy_mesh.ca_smoothing(mesh, options['angle'], - options['max distance'], - options['min weight'], - options['steps']) - # polydata = mesh.to_vtk() - - # polydata.SetSource(None) - # polydata.DebugOn() - else: - #smoother = vtk.vtkWindowedSincPolyDataFilter() - smoother = vtk.vtkSmoothPolyDataFilter() - smoother_ref = weakref.ref(smoother) - smoother_ref().AddObserver("ProgressEvent", lambda obj,evt: - UpdateProgress(smoother_ref(), _("Creating 3D surface..."))) - smoother.SetInputData(polydata) - smoother.SetNumberOfIterations(smooth_iterations) - smoother.SetRelaxationFactor(smooth_relaxation_factor) - smoother.SetFeatureAngle(80) - #smoother.SetEdgeAngle(90.0) - #smoother.SetPassBand(0.1) - smoother.BoundarySmoothingOn() - smoother.FeatureEdgeSmoothingOn() - #smoother.NormalizeCoordinatesOn() - #smoother.NonManifoldSmoothingOn() - # smoother.ReleaseDataFlagOn() - # smoother.GetOutput().ReleaseDataFlagOn() - smoother.Update() - del polydata - polydata = smoother.GetOutput() - #polydata.Register(None) - # polydata.SetSource(None) - del smoother - - - if decimate_reduction: - print("Decimating", decimate_reduction) - decimation = vtk.vtkQuadricDecimation() - # decimation.ReleaseDataFlagOn() - decimation.SetInputData(polydata) - decimation.SetTargetReduction(decimate_reduction) - decimation_ref = weakref.ref(decimation) - decimation_ref().AddObserver("ProgressEvent", lambda obj,evt: - UpdateProgress(decimation_ref(), _("Creating 3D surface..."))) - #decimation.PreserveTopologyOn() - #decimation.SplittingOff() - #decimation.BoundaryVertexDeletionOff() - # decimation.GetOutput().ReleaseDataFlagOn() - decimation.Update() - del polydata - polydata = decimation.GetOutput() - #polydata.Register(None) - # polydata.SetSource(None) - del decimation - - #to_measure.Register(None) - # to_measure.SetSource(None) - if keep_largest: - conn = vtk.vtkPolyDataConnectivityFilter() - conn.SetInputData(polydata) - conn.SetExtractionModeToLargestRegion() - conn_ref = weakref.ref(conn) - conn_ref().AddObserver("ProgressEvent", lambda obj,evt: - UpdateProgress(conn_ref(), _("Creating 3D surface..."))) - conn.Update() - # conn.GetOutput().ReleaseDataFlagOn() - del polydata - polydata = conn.GetOutput() - #polydata.Register(None) - # polydata.SetSource(None) - del conn - - #Filter used to detect and fill holes. Only fill boundary edges holes. - #TODO: Hey! This piece of code is the same from - #polydata_utils.FillSurfaceHole, we need to review this. - if fill_holes: - filled_polydata = vtk.vtkFillHolesFilter() - # filled_polydata.ReleaseDataFlagOn() - filled_polydata.SetInputData(polydata) - filled_polydata.SetHoleSize(300) - filled_polydata_ref = weakref.ref(filled_polydata) - filled_polydata_ref().AddObserver("ProgressEvent", lambda obj,evt: - UpdateProgress(filled_polydata_ref(), _("Creating 3D surface..."))) - filled_polydata.Update() - # filled_polydata.GetOutput().ReleaseDataFlagOn() - del polydata - polydata = filled_polydata.GetOutput() - #polydata.Register(None) - # polydata.SetSource(None) - # polydata.DebugOn() - del filled_polydata - - to_measure = polydata - - # If InVesalius is running without GUI - if wx.GetApp() is None: proj = prj.Project() #Create Surface instance if overwrite: @@ -720,115 +636,108 @@ class SurfaceManager(): index = proj.AddSurface(surface) surface.index = index self.last_surface_index = index + surface.colour = colour surface.polydata = polydata + surface.volume = surface_measures['volume'] + surface.area = surface_measures['area'] # With GUI else: - normals = vtk.vtkPolyDataNormals() - # normals.ReleaseDataFlagOn() - normals_ref = weakref.ref(normals) - normals_ref().AddObserver("ProgressEvent", lambda obj,evt: - UpdateProgress(normals_ref(), _("Creating 3D surface..."))) - normals.SetInputData(polydata) - normals.SetFeatureAngle(80) - normals.AutoOrientNormalsOn() - # normals.GetOutput().ReleaseDataFlagOn() - normals.Update() - del polydata - polydata = normals.GetOutput() - #polydata.Register(None) - # polydata.SetSource(None) - del normals - - # Improve performance - stripper = vtk.vtkStripper() - # stripper.ReleaseDataFlagOn() - stripper_ref = weakref.ref(stripper) - stripper_ref().AddObserver("ProgressEvent", lambda obj,evt: - UpdateProgress(stripper_ref(), _("Creating 3D surface..."))) - stripper.SetInputData(polydata) - stripper.PassThroughCellIdsOn() - stripper.PassThroughPointIdsOn() - # stripper.GetOutput().ReleaseDataFlagOn() - stripper.Update() - del polydata - polydata = stripper.GetOutput() - #polydata.Register(None) - # polydata.SetSource(None) - del stripper - - # Map polygonal data (vtkPolyData) to graphics primitives. - mapper = vtk.vtkPolyDataMapper() - mapper.SetInputData(polydata) - mapper.ScalarVisibilityOff() - # mapper.ReleaseDataFlagOn() - mapper.ImmediateModeRenderingOn() # improve performance - - # Represent an object (geometry & properties) in the rendered scene - actor = vtk.vtkActor() - actor.SetMapper(mapper) - del mapper - #Create Surface instance - if overwrite: - surface = Surface(index = self.last_surface_index) - else: - surface = Surface(name=surface_name) - surface.colour = colour - surface.polydata = polydata - del polydata - - # Set actor colour and transparency - actor.GetProperty().SetColor(colour) - actor.GetProperty().SetOpacity(1-surface.transparency) - - prop = actor.GetProperty() - - interpolation = int(ses.Session().surface_interpolation) - - prop.SetInterpolation(interpolation) - - proj = prj.Project() - if overwrite: - proj.ChangeSurface(surface) - else: - index = proj.AddSurface(surface) - surface.index = index - self.last_surface_index = index - - session = ses.Session() - session.ChangeProject() - - measured_polydata = vtk.vtkMassProperties() - # measured_polydata.ReleaseDataFlagOn() - measured_polydata.SetInputData(to_measure) - volume = float(measured_polydata.GetVolume()) - area = float(measured_polydata.GetSurfaceArea()) - surface.volume = volume - surface.area = area - self.last_surface_index = surface.index - del measured_polydata - del to_measure - - Publisher.sendMessage('Load surface actor into viewer', actor=actor) - - # Send actor by pubsub to viewer's render - if overwrite and self.actors_dict.keys(): - old_actor = self.actors_dict[self.last_surface_index] - Publisher.sendMessage('Remove surface actor from viewer', actor=old_actor) - - # Save actor for future management tasks - self.actors_dict[surface.index] = actor - - Publisher.sendMessage('Update surface info in GUI', surface=surface) - - #When you finalize the progress. The bar is cleaned. - UpdateProgress = vu.ShowProgress(1) - UpdateProgress(0, _("Ready")) - Publisher.sendMessage('Update status text in GUI', label=_("Ready")) - - Publisher.sendMessage('End busy cursor') - del actor + sp = dialogs.SurfaceProgressWindow() + for i in range(n_pieces): + init = i * piece_size + end = init + piece_size + o_piece + roi = slice(init, end) + print("new_piece", roi) + try: + f = pool.apply_async(surface_process.create_surface_piece, + args = (filename_img, matrix.shape, matrix.dtype, + mask.temp_file, mask.matrix.shape, + mask.matrix.dtype, roi, spacing, mode, + min_value, max_value, decimate_reduction, + smooth_relaxation_factor, + smooth_iterations, language, flip_image, + algorithm != 'Default', algorithm, + imagedata_resolution), + callback=lambda x: filenames.append(x), + error_callback=functools.partial(self._on_callback_error, + dialog=sp)) + # python2 + except TypeError: + f = pool.apply_async(surface_process.create_surface_piece, + args = (filename_img, matrix.shape, matrix.dtype, + mask.temp_file, mask.matrix.shape, + mask.matrix.dtype, roi, spacing, mode, + min_value, max_value, decimate_reduction, + smooth_relaxation_factor, + smooth_iterations, language, flip_image, + algorithm != 'Default', algorithm, + imagedata_resolution), + callback=lambda x: filenames.append(x)) + + while len(filenames) != n_pieces: + if sp.WasCancelled() or not sp.running: + break + time.sleep(0.25) + sp.Update(_("Creating 3D surface...")) + wx.Yield() + + if not sp.WasCancelled() or sp.running: + try: + f = pool.apply_async(surface_process.join_process_surface, + args=(filenames, algorithm, smooth_iterations, + smooth_relaxation_factor, + decimate_reduction, keep_largest, + fill_holes, options, msg_queue), + callback=functools.partial(self._on_complete_surface_creation, + overwrite=overwrite, + surface_name=surface_name, + colour=colour, + dialog=sp), + error_callback=functools.partial(self._on_callback_error, + dialog=sp)) + # python2 + except TypeError: + f = pool.apply_async(surface_process.join_process_surface, + args=(filenames, algorithm, smooth_iterations, + smooth_relaxation_factor, + decimate_reduction, keep_largest, + fill_holes, options, msg_queue), + callback=functools.partial(self._on_complete_surface_creation, + overwrite=overwrite, + surface_name=surface_name, + colour=colour, + dialog=sp)) + + while sp.running: + if sp.WasCancelled(): + break + time.sleep(0.25) + try: + msg = msg_queue.get_nowait() + sp.Update(msg) + except: + sp.Update(None) + wx.Yield() + + t_end = time.time() + print("Elapsed time - {}".format(t_end-t_init)) + sp.Close() + if sp.error: + dlg = GMD.GenericMessageDialog(None, sp.error, + "Exception!", + wx.OK|wx.ICON_ERROR) + dlg.ShowModal() + del sp + + pool.close() + pool.terminate() + del pool + del manager + del msg_queue + import gc + gc.collect() def UpdateSurfaceInterpolation(self): interpolation = int(ses.Session().surface_interpolation) diff --git a/invesalius/data/surface_process.py b/invesalius/data/surface_process.py index 1f325ea..171e109 100644 --- a/invesalius/data/surface_process.py +++ b/invesalius/data/surface_process.py @@ -1,14 +1,22 @@ import multiprocessing +import os import tempfile import time +try: + import queue +except ImportError: + import Queue as queue + import numpy import vtk import invesalius.i18n as i18n import invesalius.data.converters as converters +from invesalius.data import cy_mesh # import invesalius.data.imagedata_utils as iu +import weakref from scipy import ndimage # TODO: Code duplicated from file {imagedata_utils.py}. @@ -32,182 +40,329 @@ def ResampleImage3D(imagedata, value): return resample.GetOutput() -class SurfaceProcess(multiprocessing.Process): - - def __init__(self, pipe, filename, shape, dtype, mask_filename, - mask_shape, mask_dtype, spacing, mode, min_value, max_value, - decimate_reduction, smooth_relaxation_factor, - smooth_iterations, language, flip_image, q_in, q_out, - from_binary, algorithm, imagedata_resolution): - - multiprocessing.Process.__init__(self) - self.pipe = pipe - self.spacing = spacing - self.filename = filename - self.mode = mode - self.min_value = min_value - self.max_value = max_value - self.decimate_reduction = decimate_reduction - self.smooth_relaxation_factor = smooth_relaxation_factor - self.smooth_iterations = smooth_iterations - self.language = language - self.flip_image = flip_image - self.q_in = q_in - self.q_out = q_out - self.dtype = dtype - self.shape = shape - self.from_binary = from_binary - self.algorithm = algorithm - self.imagedata_resolution = imagedata_resolution - - self.mask_filename = mask_filename - self.mask_shape = mask_shape - self.mask_dtype = mask_dtype - - def run(self): - if self.from_binary: - self.mask = numpy.memmap(self.mask_filename, mode='r', - dtype=self.mask_dtype, - shape=self.mask_shape) - else: - self.image = numpy.memmap(self.filename, mode='r', dtype=self.dtype, - shape=self.shape) - self.mask = numpy.memmap(self.mask_filename, mode='r', - dtype=self.mask_dtype, - shape=self.mask_shape) - - while 1: - roi = self.q_in.get() - if roi is None: - break - self.CreateSurface(roi) - - def SendProgress(self, obj, msg): - prog = obj.GetProgress() - self.pipe.send([prog, msg]) - - def CreateSurface(self, roi): - if self.from_binary: - a_mask = numpy.array(self.mask[roi.start + 1: roi.stop + 1, + +def create_surface_piece(filename, shape, dtype, mask_filename, mask_shape, + mask_dtype, roi, spacing, mode, min_value, max_value, + decimate_reduction, smooth_relaxation_factor, + smooth_iterations, language, flip_image, + from_binary, algorithm, imagedata_resolution): + if from_binary: + mask = numpy.memmap(mask_filename, mode='r', + dtype=mask_dtype, + shape=mask_shape) + a_mask = numpy.array(mask[roi.start + 1: roi.stop + 1, + 1:, 1:]) + image = converters.to_vtk(a_mask, spacing, roi.start, + "AXIAL") + del a_mask + else: + image = numpy.memmap(filename, mode='r', dtype=dtype, + shape=shape) + mask = numpy.memmap(mask_filename, mode='r', + dtype=mask_dtype, + shape=mask_shape) + a_image = numpy.array(image[roi]) + + if algorithm == u'InVesalius 3.b2': + a_mask = numpy.array(mask[roi.start + 1: roi.stop + 1, 1:, 1:]) - image = converters.to_vtk(a_mask, self.spacing, roi.start, + a_image[a_mask == 1] = a_image.min() - 1 + a_image[a_mask == 254] = (min_value + max_value) / 2.0 + + image = converters.to_vtk(a_image, spacing, roi.start, "AXIAL") + + gauss = vtk.vtkImageGaussianSmooth() + gauss.SetInputData(image) + gauss.SetRadiusFactor(0.3) + gauss.ReleaseDataFlagOn() + gauss.Update() + + del image + image = gauss.GetOutput() + del gauss del a_mask else: - a_image = numpy.array(self.image[roi]) - - if self.algorithm == u'InVesalius 3.b2': - a_mask = numpy.array(self.mask[roi.start + 1: roi.stop + 1, - 1:, 1:]) - a_image[a_mask == 1] = a_image.min() - 1 - a_image[a_mask == 254] = (self.min_value + self.max_value) / 2.0 - - image = converters.to_vtk(a_image, self.spacing, roi.start, - "AXIAL") - - gauss = vtk.vtkImageGaussianSmooth() - gauss.SetInputData(image) - gauss.SetRadiusFactor(0.3) - gauss.ReleaseDataFlagOn() - gauss.Update() - - del image - image = gauss.GetOutput() - del gauss - del a_mask - else: - image = converters.to_vtk(a_image, self.spacing, roi.start, - "AXIAL") - del a_image - - if self.imagedata_resolution: - # image = iu.ResampleImage3D(image, self.imagedata_resolution) - image = ResampleImage3D(image, self.imagedata_resolution) - - flip = vtk.vtkImageFlip() - flip.SetInputData(image) - flip.SetFilteredAxis(1) - flip.FlipAboutOriginOn() - flip.ReleaseDataFlagOn() - flip.Update() - - del image - image = flip.GetOutput() - del flip - - #filename = tempfile.mktemp(suffix='_%s.vti' % (self.pid)) - #writer = vtk.vtkXMLImageDataWriter() - #writer.SetInput(mask_vtk) - #writer.SetFileName(filename) - #writer.Write() - - #print "Writing piece", roi, "to", filename - - # Create vtkPolyData from vtkImageData - #print "Generating Polydata" - #if self.mode == "CONTOUR": - #print "Contour" - contour = vtk.vtkContourFilter() - contour.SetInputData(image) - #contour.SetInput(flip.GetOutput()) - if self.from_binary: - contour.SetValue(0, 127) # initial threshold - else: - contour.SetValue(0, self.min_value) # initial threshold - contour.SetValue(1, self.max_value) # final threshold - contour.ComputeScalarsOn() - contour.ComputeGradientsOn() - contour.ComputeNormalsOn() - contour.ReleaseDataFlagOn() - contour.Update() - #contour.AddObserver("ProgressEvent", lambda obj,evt: - # self.SendProgress(obj, _("Generating 3D surface..."))) - polydata = contour.GetOutput() - del image - del contour - - #else: #mode == "GRAYSCALE": - #mcubes = vtk.vtkMarchingCubes() - #mcubes.SetInput(flip.GetOutput()) - #mcubes.SetValue(0, self.min_value) - #mcubes.SetValue(1, self.max_value) - #mcubes.ComputeScalarsOff() - #mcubes.ComputeGradientsOff() - #mcubes.ComputeNormalsOff() - #mcubes.AddObserver("ProgressEvent", lambda obj,evt: - #self.SendProgress(obj, _("Generating 3D surface..."))) - #polydata = mcubes.GetOutput() - - #triangle = vtk.vtkTriangleFilter() - #triangle.SetInput(polydata) - #triangle.AddObserver("ProgressEvent", lambda obj,evt: - #self.SendProgress(obj, _("Generating 3D surface..."))) - #triangle.Update() - #polydata = triangle.GetOutput() - - #if self.decimate_reduction: - - ##print "Decimating" - #decimation = vtk.vtkDecimatePro() - #decimation.SetInput(polydata) - #decimation.SetTargetReduction(0.3) - #decimation.AddObserver("ProgressEvent", lambda obj,evt: - #self.SendProgress(obj, _("Generating 3D surface..."))) - ##decimation.PreserveTopologyOn() - #decimation.SplittingOff() - #decimation.BoundaryVertexDeletionOff() - #polydata = decimation.GetOutput() - - self.pipe.send(None) - - filename = tempfile.mktemp(suffix='_%s.vtp' % (self.pid)) - writer = vtk.vtkXMLPolyDataWriter() - writer.SetInputData(polydata) - writer.SetFileName(filename) - writer.Write() - - print("Writing piece", roi, "to", filename) + image = converters.to_vtk(a_image, spacing, roi.start, + "AXIAL") + del a_image + + if imagedata_resolution: + image = ResampleImage3D(image, imagedata_resolution) + + flip = vtk.vtkImageFlip() + flip.SetInputData(image) + flip.SetFilteredAxis(1) + flip.FlipAboutOriginOn() + flip.ReleaseDataFlagOn() + flip.Update() + + del image + image = flip.GetOutput() + del flip + + contour = vtk.vtkContourFilter() + contour.SetInputData(image) + if from_binary: + contour.SetValue(0, 127) # initial threshold + else: + contour.SetValue(0, min_value) # initial threshold + contour.SetValue(1, max_value) # final threshold + # contour.ComputeScalarsOn() + # contour.ComputeGradientsOn() + # contour.ComputeNormalsOn() + contour.ReleaseDataFlagOn() + contour.Update() + + polydata = contour.GetOutput() + del image + del contour + + filename = tempfile.mktemp(suffix='_%d_%d.vtp' % (roi.start, roi.stop)) + writer = vtk.vtkXMLPolyDataWriter() + writer.SetInputData(polydata) + writer.SetFileName(filename) + writer.Write() + + print("Writing piece", roi, "to", filename) + print("MY PID MC", os.getpid()) + return filename + + +def join_process_surface(filenames, algorithm, smooth_iterations, smooth_relaxation_factor, decimate_reduction, keep_largest, fill_holes, options, msg_queue): + def send_message(msg): + try: + msg_queue.put_nowait(msg) + except queue.Full as e: + print(e) + + send_message('Joining surfaces ...') + polydata_append = vtk.vtkAppendPolyData() + for f in filenames: + reader = vtk.vtkXMLPolyDataReader() + reader.SetFileName(f) + reader.Update() + + polydata = reader.GetOutput() + + polydata_append.AddInputData(polydata) + del reader del polydata - del writer - self.q_out.put(filename) + polydata_append.Update() + # polydata_append.GetOutput().ReleaseDataFlagOn() + polydata = polydata_append.GetOutput() + #polydata.Register(None) + # polydata.SetSource(None) + del polydata_append + + send_message('Cleaning surface ...') + clean = vtk.vtkCleanPolyData() + # clean.ReleaseDataFlagOn() + # clean.GetOutput().ReleaseDataFlagOn() + clean_ref = weakref.ref(clean) + # clean_ref().AddObserver("ProgressEvent", lambda obj,evt: + # UpdateProgress(clean_ref(), _("Creating 3D surface..."))) + clean.SetInputData(polydata) + clean.PointMergingOn() + clean.Update() + + del polydata + polydata = clean.GetOutput() + # polydata.SetSource(None) + del clean + + if algorithm == 'ca_smoothing': + send_message('Calculating normals ...') + normals = vtk.vtkPolyDataNormals() + normals_ref = weakref.ref(normals) + # normals_ref().AddObserver("ProgressEvent", lambda obj,evt: + # UpdateProgress(normals_ref(), _("Creating 3D surface..."))) + normals.SetInputData(polydata) + # normals.ReleaseDataFlagOn() + #normals.SetFeatureAngle(80) + #normals.AutoOrientNormalsOn() + normals.ComputeCellNormalsOn() + # normals.GetOutput().ReleaseDataFlagOn() + normals.Update() + del polydata + polydata = normals.GetOutput() + # polydata.SetSource(None) + del normals + + clean = vtk.vtkCleanPolyData() + # clean.ReleaseDataFlagOn() + # clean.GetOutput().ReleaseDataFlagOn() + clean_ref = weakref.ref(clean) + # clean_ref().AddObserver("ProgressEvent", lambda obj,evt: + # UpdateProgress(clean_ref(), _("Creating 3D surface..."))) + clean.SetInputData(polydata) + clean.PointMergingOn() + clean.Update() + + del polydata + polydata = clean.GetOutput() + # polydata.SetSource(None) + del clean + + # try: + # polydata.BuildLinks() + # except TypeError: + # polydata.BuildLinks(0) + # polydata = ca_smoothing.ca_smoothing(polydata, options['angle'], + # options['max distance'], + # options['min weight'], + # options['steps']) + + send_message('Context Aware smoothing ...') + mesh = cy_mesh.Mesh(polydata) + cy_mesh.ca_smoothing(mesh, options['angle'], + options['max distance'], + options['min weight'], + options['steps']) + # polydata = mesh.to_vtk() + + # polydata.SetSource(None) + # polydata.DebugOn() + else: + #smoother = vtk.vtkWindowedSincPolyDataFilter() + send_message('Smoothing ...') + smoother = vtk.vtkSmoothPolyDataFilter() + smoother_ref = weakref.ref(smoother) + # smoother_ref().AddObserver("ProgressEvent", lambda obj,evt: + # UpdateProgress(smoother_ref(), _("Creating 3D surface..."))) + smoother.SetInputData(polydata) + smoother.SetNumberOfIterations(smooth_iterations) + smoother.SetRelaxationFactor(smooth_relaxation_factor) + smoother.SetFeatureAngle(80) + #smoother.SetEdgeAngle(90.0) + #smoother.SetPassBand(0.1) + smoother.BoundarySmoothingOn() + smoother.FeatureEdgeSmoothingOn() + #smoother.NormalizeCoordinatesOn() + #smoother.NonManifoldSmoothingOn() + # smoother.ReleaseDataFlagOn() + # smoother.GetOutput().ReleaseDataFlagOn() + smoother.Update() + del polydata + polydata = smoother.GetOutput() + #polydata.Register(None) + # polydata.SetSource(None) + del smoother + + + if decimate_reduction: + print("Decimating", decimate_reduction) + send_message('Decimating ...') + decimation = vtk.vtkQuadricDecimation() + # decimation.ReleaseDataFlagOn() + decimation.SetInputData(polydata) + decimation.SetTargetReduction(decimate_reduction) + decimation_ref = weakref.ref(decimation) + # decimation_ref().AddObserver("ProgressEvent", lambda obj,evt: + # UpdateProgress(decimation_ref(), _("Creating 3D surface..."))) + #decimation.PreserveTopologyOn() + #decimation.SplittingOff() + #decimation.BoundaryVertexDeletionOff() + # decimation.GetOutput().ReleaseDataFlagOn() + decimation.Update() + del polydata + polydata = decimation.GetOutput() + #polydata.Register(None) + # polydata.SetSource(None) + del decimation + + #to_measure.Register(None) + # to_measure.SetSource(None) + + if keep_largest: + send_message('Finding the largest ...') + conn = vtk.vtkPolyDataConnectivityFilter() + conn.SetInputData(polydata) + conn.SetExtractionModeToLargestRegion() + conn_ref = weakref.ref(conn) + # conn_ref().AddObserver("ProgressEvent", lambda obj,evt: + # UpdateProgress(conn_ref(), _("Creating 3D surface..."))) + conn.Update() + # conn.GetOutput().ReleaseDataFlagOn() + del polydata + polydata = conn.GetOutput() + #polydata.Register(None) + # polydata.SetSource(None) + del conn + + #Filter used to detect and fill holes. Only fill boundary edges holes. + #TODO: Hey! This piece of code is the same from + #polydata_utils.FillSurfaceHole, we need to review this. + if fill_holes: + send_message('Filling holes ...') + filled_polydata = vtk.vtkFillHolesFilter() + # filled_polydata.ReleaseDataFlagOn() + filled_polydata.SetInputData(polydata) + filled_polydata.SetHoleSize(300) + filled_polydata_ref = weakref.ref(filled_polydata) + # filled_polydata_ref().AddObserver("ProgressEvent", lambda obj,evt: + # UpdateProgress(filled_polydata_ref(), _("Creating 3D surface..."))) + filled_polydata.Update() + # filled_polydata.GetOutput().ReleaseDataFlagOn() + del polydata + polydata = filled_polydata.GetOutput() + #polydata.Register(None) + # polydata.SetSource(None) + # polydata.DebugOn() + del filled_polydata + + to_measure = polydata + + normals = vtk.vtkPolyDataNormals() + # normals.ReleaseDataFlagOn() + # normals_ref = weakref.ref(normals) + # normals_ref().AddObserver("ProgressEvent", lambda obj,evt: + # UpdateProgress(normals_ref(), _("Creating 3D surface..."))) + normals.SetInputData(polydata) + # normals.SetFeatureAngle(80) + # normals.SplittingOff() + # normals.AutoOrientNormalsOn() + # normals.GetOutput().ReleaseDataFlagOn() + normals.Update() + del polydata + polydata = normals.GetOutput() + #polydata.Register(None) + # polydata.SetSource(None) + del normals + + + # Improve performance + stripper = vtk.vtkStripper() + # stripper.ReleaseDataFlagOn() + # stripper_ref = weakref.ref(stripper) + # stripper_ref().AddObserver("ProgressEvent", lambda obj,evt: + # UpdateProgress(stripper_ref(), _("Creating 3D surface..."))) + stripper.SetInputData(polydata) + stripper.PassThroughCellIdsOn() + stripper.PassThroughPointIdsOn() + # stripper.GetOutput().ReleaseDataFlagOn() + stripper.Update() + del polydata + polydata = stripper.GetOutput() + #polydata.Register(None) + # polydata.SetSource(None) + del stripper + + send_message('Calculating area and volume ...') + measured_polydata = vtk.vtkMassProperties() + measured_polydata.SetInputData(to_measure) + measured_polydata.Update() + volume = float(measured_polydata.GetVolume()) + area = float(measured_polydata.GetSurfaceArea()) + del measured_polydata + + filename = tempfile.mktemp(suffix='_full.vtp') + writer = vtk.vtkXMLPolyDataWriter() + writer.SetInputData(polydata) + writer.SetFileName(filename) + writer.Write() + del writer + + print("MY PID", os.getpid()) + return filename, {'volume': volume, 'area': area} diff --git a/invesalius/gui/dialogs.py b/invesalius/gui/dialogs.py index aa2dde1..891fa76 100644 --- a/invesalius/gui/dialogs.py +++ b/invesalius/gui/dialogs.py @@ -3523,3 +3523,30 @@ class ObjectCalibrationDialog(wx.Dialog): def GetValue(self): return self.obj_fiducials, self.obj_orients, self.obj_ref_id, self.obj_name + + +class SurfaceProgressWindow(object): + def __init__(self): + self.title = "InVesalius 3" + self.msg = _("Creating 3D surface ...") + self.style = wx.PD_APP_MODAL | wx.PD_APP_MODAL | wx.PD_CAN_ABORT | wx.PD_ELAPSED_TIME + self.dlg = wx.ProgressDialog(self.title, + self.msg, + parent=None, + style=self.style) + self.running = True + self.error = None + self.dlg.Show() + + def WasCancelled(self): + # print("Cancelled?", self.dlg.WasCancelled()) + return self.dlg.WasCancelled() + + def Update(self, msg=None, value=None): + if msg is None: + self.dlg.Pulse() + else: + self.dlg.Pulse(msg) + + def Close(self): + self.dlg.Destroy() diff --git a/invesalius/utils.py b/invesalius/utils.py index 9787fbe..0bdc615 100644 --- a/invesalius/utils.py +++ b/invesalius/utils.py @@ -22,8 +22,10 @@ import sys import re import locale import math +import traceback from distutils.version import LooseVersion +from functools import wraps import numpy as np @@ -462,3 +464,24 @@ def encode(text, encoding, *args): return text.encode(encoding, *args) except AttributeError: return text + + +def timing(f): + @wraps(f) + def wrapper(*args, **kwargs): + start = time.time() + result = f(*args, **kwargs) + end = time.time() + print('{} elapsed time: {}'.format(f.__name__, end-start)) + return result + return wrapper + + +def log_traceback(ex): + if hasattr(ex, '__traceback__'): + ex_traceback = ex.__traceback__ + else: + _, _, ex_traceback = sys.exc_info() + tb_lines = [line.rstrip('\n') for line in + traceback.format_exception(ex.__class__, ex, ex_traceback)] + return ''.join(tb_lines) -- libgit2 0.21.2