diff --git a/invesalius/control.py b/invesalius/control.py index b7a3c85..ff270b9 100755 --- a/invesalius/control.py +++ b/invesalius/control.py @@ -464,7 +464,7 @@ class Controller(): filename = filename.replace("/", "") #Fix problem case other/Skull_DICOM dirpath = session.CreateProject(filename) - proj.SavePlistProject(dirpath, filename) + #proj.SavePlistProject(dirpath, filename) def OnOpenDicomGroup(self, pubsub_evt): group, interval, file_range = pubsub_evt.data diff --git a/invesalius/data/mask.py b/invesalius/data/mask.py index da2353f..62ef4ac 100644 --- a/invesalius/data/mask.py +++ b/invesalius/data/mask.py @@ -42,6 +42,7 @@ class Mask(): self.edition_threshold_range = [const.THRESHOLD_OUTVALUE, const.THRESHOLD_INVALUE] self.is_shown = 1 self.edited_points = {} + self.was_edited = False def SavePlist(self, filename): mask = {} diff --git a/invesalius/data/slice_.py b/invesalius/data/slice_.py index ef8841c..f3e3688 100644 --- a/invesalius/data/slice_.py +++ b/invesalius/data/slice_.py @@ -285,6 +285,7 @@ class Slice(object): mask = self.buffer_slices[orientation].mask image = self.buffer_slices[orientation].image thresh_min, thresh_max = self.current_mask.edition_threshold_range + self.current_mask.was_edited = True if hasattr(position, '__iter__'): py, px = position @@ -515,7 +516,9 @@ class Slice(object): If slice_number is None then all the threshold is calculated for all slices, otherwise only to indicated slice. """ + self.current_mask.was_edited = False thresh_min, thresh_max = threshold_range + print "Threshold" if self.current_mask.index == index: # TODO: find out a better way to do threshold @@ -588,23 +591,26 @@ class Slice(object): #--------------------------------------------------------------------------- def CreateSurfaceFromIndex(self, pubsub_evt): - mask_index, overwrite_surface = pubsub_evt.data - + mask_index, overwrite_surface, algorithm, options = pubsub_evt.data proj = Project() mask = proj.mask_dict[mask_index] # This is very important. Do not use masks' imagedata. It would mess up # surface quality event when using contour - colour = mask.colour - threshold = mask.threshold_range - edited_points = mask.edited_points - - self.SetMaskThreshold(mask.index, threshold) + #self.SetMaskThreshold(mask.index, threshold) + for n in xrange(1, mask.matrix.shape[0]): + if mask.matrix[n, 0, 0] == 0: + m = mask.matrix[n, 1:, 1:] + mask.matrix[n, 1:, 1:] = self.do_threshold_to_a_slice(self.matrix[n-1], m) mask.matrix.flush() - ps.Publisher().sendMessage('Create surface', (self.matrix, self.matrix_filename, mask, self.spacing)) + ps.Publisher().sendMessage('Create surface', (algorithm, options, + self.matrix, + self.matrix_filename, + mask, self.spacing, + overwrite_surface)) def GetOutput(self): return self.blend_filter.GetOutput() diff --git a/invesalius/data/surface.py b/invesalius/data/surface.py index f05352c..65cc758 100644 --- a/invesalius/data/surface.py +++ b/invesalius/data/surface.py @@ -21,8 +21,6 @@ import multiprocessing import os import plistlib import random -import sys -import tempfile import vtk import wx.lib.pubsub as ps @@ -35,6 +33,7 @@ import session as ses import surface_process import utils as utl import vtk_utils as vu +import data.ca_smoothing as ca_smoothing class Surface(): """ @@ -363,9 +362,11 @@ class SurfaceManager(): """ Create surface actor, save into project and send it to viewer. """ - matrix, filename_img, mask, spacing = pubsub_evt.data + algorithm, options, matrix, filename_img, mask, spacing, overwrite = pubsub_evt.data min_value, max_value = mask.threshold_range - fill_holes = True + fill_holes = False + + mask.matrix.flush() #if len(surface_data) == 5: #imagedata, colour, [min_value, max_value], \ @@ -382,7 +383,7 @@ class SurfaceManager(): mode = 'CONTOUR' # 'GRAYSCALE' quality=_('Optimal *') - keep_largest = True + keep_largest = False surface_name = "" colour = mask.colour @@ -411,8 +412,6 @@ class SurfaceManager(): language = ses.Session().language - overwrite = 0 - if (prj.Project().original_orientation == const.CORONAL): flip_image = False else: @@ -422,7 +421,7 @@ class SurfaceManager(): pipe_in, pipe_out = multiprocessing.Pipe() o_piece = 1 - piece_size = 20 + piece_size = 2000 n_pieces = int(round(matrix.shape[0] / piece_size + 0.5, 0)) @@ -432,10 +431,18 @@ class SurfaceManager(): p = [] for i in xrange(n_processors): sp = surface_process.SurfaceProcess(pipe_in, filename_img, - matrix.shape, matrix.dtype, spacing, - mode, min_value, max_value, - decimate_reduction, smooth_relaxation_factor, - smooth_iterations, language, flip_image, q_in, q_out) + 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, + mask.was_edited and algorithm != u'InVesalius 3.b2', + algorithm) p.append(sp) sp.start() @@ -471,28 +478,65 @@ class SurfaceManager(): polydata_append.AddInput(reader.GetOutput()) t -= 1 + polydata_append.Update() polydata = polydata_append.GetOutput() - clean = vtk.vtkCleanPolyData() - clean.AddObserver("ProgressEvent", lambda obj,evt: - UpdateProgress(obj, _("Generating 3D surface..."))) - clean.SetInput(polydata) - clean.PointMergingOn() - clean.Update() - polydata = clean.GetOutput() - smoother = vtk.vtkWindowedSincPolyDataFilter() - smoother.AddObserver("ProgressEvent", lambda obj,evt: - UpdateProgress(obj, _("Generating 3D surface..."))) - smoother.SetInput(polydata) - smoother.SetNumberOfIterations(smooth_iterations) - smoother.SetFeatureAngle(120) - smoother.BoundarySmoothingOn() - smoother.SetPassBand(0.1) - #smoother.FeatureEdgeSmoothingOn() - #smoother.NonManifoldSmoothingOn() - #smoother.NormalizeCoordinatesOn() - smoother.Update() - polydata = smoother.GetOutput() + if algorithm == u'Context aware smoothing': + normals = vtk.vtkPolyDataNormals() + normals.AddObserver("ProgressEvent", lambda obj,evt: + UpdateProgress(obj, _("Generating 3D surface..."))) + normals.SetInput(polydata) + #normals.SetFeatureAngle(80) + #normals.AutoOrientNormalsOn() + normals.ComputeCellNormalsOn() + normals.Update() + polydata = normals.GetOutput() + + clean = vtk.vtkCleanPolyData() + clean.AddObserver("ProgressEvent", lambda obj,evt: + UpdateProgress(obj, _("Generating 3D surface..."))) + clean.SetInput(polydata) + clean.PointMergingOn() + clean.Update() + polydata = clean.GetOutput() + + polydata.BuildLinks() + polydata = ca_smoothing.ca_smoothing(polydata, options['angle'], + options['max distance'], + options['min weight'], + options['steps']) + + else: + smoother = vtk.vtkWindowedSincPolyDataFilter() + smoother.AddObserver("ProgressEvent", lambda obj,evt: + UpdateProgress(obj, _("Generating 3D surface..."))) + smoother.SetInput(polydata) + smoother.SetNumberOfIterations(smooth_iterations) + smoother.SetFeatureAngle(120) + smoother.SetEdgeAngle(90.0) + smoother.BoundarySmoothingOn() + smoother.SetPassBand(0.1) + #smoother.FeatureEdgeSmoothingOn() + #smoother.NonManifoldSmoothingOn() + #smoother.NormalizeCoordinatesOn() + smoother.Update() + polydata = smoother.GetOutput() + + + if decimate_reduction: + print "Decimating", decimate_reduction + decimation = vtk.vtkQuadricDecimation() + decimation.SetInput(polydata) + decimation.SetTargetReduction(decimate_reduction) + decimation.AddObserver("ProgressEvent", lambda obj,evt: + UpdateProgress(obj, _("Generating 3D surface..."))) + #decimation.PreserveTopologyOn() + #decimation.SplittingOff() + #decimation.BoundaryVertexDeletionOff() + decimation.Update() + polydata = decimation.GetOutput() + + to_measure = polydata if keep_largest: conn = vtk.vtkPolyDataConnectivityFilter() @@ -500,17 +544,19 @@ class SurfaceManager(): conn.SetExtractionModeToLargestRegion() conn.AddObserver("ProgressEvent", lambda obj,evt: UpdateProgress(obj, _("Generating 3D surface..."))) + conn.Update() polydata = conn.GetOutput() - # Filter used to detect and fill holes. Only fill boundary edges holes. + #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. + #polydata_utils.FillSurfaceHole, we need to review this. if fill_holes: filled_polydata = vtk.vtkFillHolesFilter() filled_polydata.SetInput(polydata) filled_polydata.SetHoleSize(300) filled_polydata.AddObserver("ProgressEvent", lambda obj,evt: UpdateProgress(obj, _("Generating 3D surface..."))) + filled_polydata.Update() polydata = filled_polydata.GetOutput() normals = vtk.vtkPolyDataNormals() @@ -572,7 +618,7 @@ class SurfaceManager(): # The following lines have to be here, otherwise all volumes disappear measured_polydata = vtk.vtkMassProperties() - measured_polydata.SetInput(smoother.GetOutput()) + measured_polydata.SetInput(to_measure) volume = measured_polydata.GetVolume() surface.volume = volume self.last_surface_index = surface.index diff --git a/invesalius/data/surface_process.py b/invesalius/data/surface_process.py index 8d28aab..3844430 100644 --- a/invesalius/data/surface_process.py +++ b/invesalius/data/surface_process.py @@ -11,9 +11,11 @@ from scipy import ndimage class SurfaceProcess(multiprocessing.Process): - def __init__(self, pipe, filename, shape, dtype, spacing, mode, min_value, max_value, + 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): + smooth_iterations, language, flip_image, q_in, q_out, + from_binary, algorithm): multiprocessing.Process.__init__(self) self.pipe = pipe @@ -31,10 +33,26 @@ class SurfaceProcess(multiprocessing.Process): self.q_out = q_out self.dtype = dtype self.shape = shape + self.from_binary = from_binary + self.algorithm = algorithm + + self.mask_filename = mask_filename + self.mask_shape = mask_shape + self.mask_dtype = mask_dtype def run(self): - self.mask = numpy.memmap(self.filename, mode='r', dtype=self.dtype, - shape=self.shape) + 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: @@ -46,59 +64,95 @@ class SurfaceProcess(multiprocessing.Process): self.pipe.send([prog, msg]) def CreateSurface(self, roi): - smoothed = numpy.array(self.mask[roi]) - image = converters.to_vtk(smoothed, self.spacing, roi.start, + if self.from_binary: + a_mask = numpy.array(self.mask[roi.start + 1: roi.stop + 1, + 1:, 1:]) + image = converters.to_vtk(a_mask, self.spacing, roi.start, "AXIAL") + 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.SetInput(image) + gauss.SetRadiusFactor(0.3) + gauss.Update() + + image = gauss.GetOutput() + else: + image = converters.to_vtk(a_image, self.spacing, roi.start, + "AXIAL") + flip = vtk.vtkImageFlip() flip.SetInput(image) flip.SetFilteredAxis(1) flip.FlipAboutOriginOn() flip.Update() + + #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": + #if self.mode == "CONTOUR": #print "Contour" - contour = vtk.vtkContourFilter() - contour.SetInput(flip.GetOutput()) + contour = vtk.vtkContourFilter() + 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.AddObserver("ProgressEvent", lambda obj,evt: - self.SendProgress(obj, _("Generating 3D surface..."))) - polydata = contour.GetOutput() - 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: + contour.ComputeScalarsOn() + contour.ComputeGradientsOn() + contour.ComputeNormalsOn() + contour.AddObserver("ProgressEvent", lambda obj,evt: self.SendProgress(obj, _("Generating 3D surface..."))) - #decimation.PreserveTopologyOn() - decimation.SplittingOff() - decimation.BoundaryVertexDeletionOff() - polydata = decimation.GetOutput() + polydata = contour.GetOutput() + #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) @@ -108,4 +162,6 @@ class SurfaceProcess(multiprocessing.Process): writer.SetFileName(filename) writer.Write() + print "Writing piece", roi, "to", filename + self.q_out.put(filename) diff --git a/invesalius/gui/dialogs.py b/invesalius/gui/dialogs.py index 201973c..ac51199 100644 --- a/invesalius/gui/dialogs.py +++ b/invesalius/gui/dialogs.py @@ -24,6 +24,7 @@ import sys import wx from wx.lib import masked +from wx.lib.agw import floatspin from wx.lib.wordwrap import wordwrap import wx.lib.pubsub as ps @@ -852,3 +853,114 @@ def ExportPicture(type_=""): return () +class CAOptions(wx.Panel): + ''' + Options related to Context aware algorithm: + Angle: The min angle to a vertex to be considered a staircase vertex; + Max distance: The max distance a normal vertex must be to calculate its + weighting; + Min Weighting: The min weight a vertex must have; + Steps: The number of iterations the smoothing algorithm have to do. + ''' + def __init__(self, parent): + wx.Panel.__init__(self, parent, -1) + self._build_widgets() + + def _build_widgets(self): + self.angle = floatspin.FloatSpin(self, -1, value=0.7, min_val=0.0, + max_val=1.0, increment=0.1, + digits=1) + + self.max_distance = floatspin.FloatSpin(self, -1, value=3.0, min_val=0.0, + max_val=100.0, increment=0.1, + digits=2) + + self.min_weight = floatspin.FloatSpin(self, -1, value=0.2, min_val=0.0, + max_val=1.0, increment=0.1, + digits=1) + + self.steps = wx.SpinCtrl(self, -1, value='10', min=1, max=100) + + layout_sizer = wx.FlexGridSizer(rows=4, cols=2, hgap=5, vgap=5) + layout_sizer.Add(wx.StaticText(self, -1, u'Angle:'), 0, wx.EXPAND) + layout_sizer.Add(self.angle, 0, wx.EXPAND) + layout_sizer.Add(wx.StaticText(self, -1, u'Max. distance:'), 0, wx.EXPAND) + layout_sizer.Add(self.max_distance, 0, wx.EXPAND) + layout_sizer.Add(wx.StaticText(self, -1, u'Min. weight:'), 0, wx.EXPAND) + layout_sizer.Add(self.min_weight, 0, wx.EXPAND) + layout_sizer.Add(wx.StaticText(self, -1, u'N. steps:'), 0, wx.EXPAND) + layout_sizer.Add(self.steps, 0, wx.EXPAND) + + self.main_sizer = wx.StaticBoxSizer(wx.StaticBox(self, -1, 'Context aware options'), wx.VERTICAL) + self.main_sizer.Add(layout_sizer, 0, wx.EXPAND | wx.ALL, 5) + self.SetSizer(self.main_sizer) + +class SurfaceDialog(wx.Dialog): + ''' + This dialog is only shown when the mask whose surface will be generate was + edited. So far, the only options available are the choice of method to + generate the surface, Binary or `Context aware smoothing', and options from + `Context aware smoothing' + ''' + def __init__(self): + wx.Dialog.__init__(self, None, -1, u'Surface generation options') + self._build_widgets() + self._bind_wx() + + def _build_widgets(self): + btn_ok = wx.Button(self, wx.ID_OK) + btn_cancel = wx.Button(self, wx.ID_CANCEL) + btn_sizer = wx.StdDialogButtonSizer() + btn_sizer.AddButton(btn_ok) + btn_sizer.AddButton(btn_cancel) + btn_sizer.Realize() + + self.alg_types = {0: u'Context aware smoothing', 1: u'Binary', + 2: u'InVesalius 3.b2'} + self.cb_types = wx.ComboBox(self, -1, self.alg_types[0], + choices=[self.alg_types[i] for i in sorted(self.alg_types)], + style=wx.CB_READONLY) + + self.ca_options = CAOptions(self) + + method_sizer = wx.BoxSizer(wx.HORIZONTAL) + method_sizer.Add(wx.StaticText(self, -1, u'Method:'), 0, + wx.EXPAND | wx.ALL, 5) + method_sizer.Add(self.cb_types, 0, wx.EXPAND) + + self.main_sizer = wx.BoxSizer(wx.VERTICAL) + self.main_sizer.Add(method_sizer, 0, wx.EXPAND | wx.ALL, 5) + self.main_sizer.Add(self.ca_options, 0, wx.EXPAND | wx.ALL, 5) + self.main_sizer.Add(btn_sizer, 0, wx.EXPAND | wx.ALL, 5) + + self.SetSizer(self.main_sizer) + self.Fit() + + def _bind_wx(self): + self.cb_types.Bind(wx.EVT_COMBOBOX, self._set_cb_types) + + def _set_cb_types(self, evt): + print evt.GetString() + if self.alg_types[evt.GetSelection()] == self.alg_types[0]: + self.ca_options.Enable() + else: + self.ca_options.Disable() + evt.Skip() + + def GetAlgorithmSelected(self): + try: + return self.alg_types[self.cb_types.GetSelection()] + except KeyError: + return self.alg_types[0] + + def GetOptions(self): + if self.GetAlgorithmSelected() == self.alg_types[0]: + options = {'angle': self.ca_options.angle.GetValue(), + 'max distance': self.ca_options.max_distance.GetValue(), + 'min weight': self.ca_options.min_weight.GetValue(), + 'steps': self.ca_options.steps.GetValue()} + else: + options = {} + return options + + diff --git a/invesalius/gui/task_slice.py b/invesalius/gui/task_slice.py index 9bb4d01..e0d2147 100644 --- a/invesalius/gui/task_slice.py +++ b/invesalius/gui/task_slice.py @@ -25,6 +25,7 @@ import wx.lib.platebtn as pbtn import wx.lib.pubsub as ps import data.mask as mask +import data.slice_ as slice_ import constants as const import gui.dialogs as dlg import gui.widgets.gradient as grad @@ -143,10 +144,20 @@ class InnerTaskPanel(wx.Panel): def OnButtonNextTask(self, evt): overwrite = self.check_box.IsChecked() + algorithm = '' + options = {} if self.GetMaskSelected() != -1: + sl = slice_.Slice() + if sl.current_mask.was_edited: + dlgs = dlg.SurfaceDialog() + if dlgs.ShowModal() == wx.ID_OK: + algorithm = dlgs.GetAlgorithmSelected() + options = dlgs.GetOptions() + else: + return ps.Publisher().sendMessage('Create surface from index', (self.GetMaskSelected(), - overwrite)) + overwrite, algorithm, options)) ps.Publisher().sendMessage('Fold surface task') else: dlg.InexistentMask() -- libgit2 0.21.2