Create border closed (#149)

* Created a function to pad image * taking padding in consideration when converting to vtk * Better normals computation * Added option to close border holes

Create border closed (#149)
* Created a function to pad image * taking padding in consideration when converting to vtk * Better normals computation * Added option to close border holes
Thiago Franco de Moraes · GitHub
1 parent ccf02f83
Showing 6 changed files with 146 additions and 72 deletions Show diff stats
invesalius/data/converters.py
invesalius/data/imagedata_utils.py
invesalius/data/surface.py
invesalius/data/surface_process.py
invesalius/gui/dialogs.py
invesalius/gui/preferences.py
@@ -21,7 +21,7 @@ import numpy
 import vtk
 from vtk.util import numpy_support
-def to_vtk(n_array, spacing, slice_number, orientation):
+def to_vtk(n_array, spacing, slice_number, orientation, origin=(0, 0, 0), padding=(0, 0, 0)):
     if orientation == "SAGITTAL":
         orientation = "SAGITAL"
@@ -32,20 +32,22 @@ def to_vtk(n_array, spacing, slice_number, orientation):
         dy, dx = n_array.shape
         dz = 1
+    px, py, pz = padding
+
     v_image = numpy_support.numpy_to_vtk(n_array.flat)
     if orientation == 'AXIAL':
-        extent = (0, dx -1, 0, dy -1, slice_number, slice_number + dz - 1)
+        extent = (0 - px , dx -1 - px, 0 - py, dy - 1 - py, slice_number - pz, slice_number + dz - 1 - pz)
     elif orientation == 'SAGITAL':
         dx, dy, dz = dz, dx, dy
-        extent = (slice_number, slice_number + dx - 1, 0, dy - 1, 0, dz - 1)
+        extent = (slice_number - px, slice_number + dx - 1 - px, 0 - py, dy - 1 - py, 0 - pz, dz - 1 - pz)
     elif orientation == 'CORONAL':
         dx, dy, dz = dx, dz, dy
-        extent = (0, dx - 1, slice_number, slice_number + dy - 1, 0, dz - 1)
+        extent = (0 - px, dx - 1 - px, slice_number - py, slice_number + dy - 1 - py, 0 - pz, dz - 1 - pz)
     # Generating the vtkImageData
     image = vtk.vtkImageData()
-    image.SetOrigin(0, 0, 0)
+    image.SetOrigin(origin)
     image.SetSpacing(spacing)
     image.SetDimensions(dx, dy, dz)
     # SetNumberOfScalarComponents and SetScalrType were replaced by
@@ -751,3 +751,22 @@ def imgnormalize(data, srange=(0, 255)):
     datan = datan.astype(numpy.int16)
     return datan
+
+
+def pad_image(image, pad_value, pad_bottom, pad_top):
+    dz, dy, dx = image.shape
+    z_iadd = 0
+    z_eadd = 0
+    if pad_bottom:
+        z_iadd = 1
+        dz += 1
+    if pad_top:
+        z_eadd = 1
+        dz += 1
+    new_shape = dz, dy + 2, dx + 2
+
+    paded_image = numpy.empty(shape=new_shape, dtype=image.dtype)
+    paded_image[:] = pad_value
+    paded_image[z_iadd: z_iadd + image.shape[0], 1:-1, 1:-1] = image
+
+    return paded_image
@@ -538,6 +538,8 @@ class SurfaceManager():
         fill_holes = surface_parameters['options']['fill']
         keep_largest = surface_parameters['options']['keep_largest']
+        fill_border_holes = surface_parameters['options'].get('fill_border_holes', True)
+
         mode = 'CONTOUR' # 'GRAYSCALE'
         min_value, max_value = mask.threshold_range
         colour = mask.colour[:3]
@@ -598,7 +600,7 @@ class SurfaceManager():
                                              smooth_relaxation_factor,
                                              smooth_iterations, language, flip_image,
                                              algorithm != 'Default', algorithm,
-                                             imagedata_resolution),
+                                             imagedata_resolution, fill_border_holes),
                                      callback=lambda x: filenames.append(x))
             while len(filenames) != n_pieces:
@@ -659,7 +661,7 @@ class SurfaceManager():
                                                  smooth_relaxation_factor,
                                                  smooth_iterations, language, flip_image,
                                                  algorithm != 'Default', algorithm,
-                                                 imagedata_resolution),
+                                                 imagedata_resolution, fill_border_holes),
                                          callback=lambda x: filenames.append(x),
                                          error_callback=functools.partial(self._on_callback_error,
                                                                           dialog=sp))
@@ -673,7 +675,7 @@ class SurfaceManager():
                                                  smooth_relaxation_factor,
                                                  smooth_iterations, language, flip_image,
                                                  algorithm != 'Default', algorithm,
-                                                 imagedata_resolution),
+                                                 imagedata_resolution, fill_border_holes),
                                          callback=lambda x: filenames.append(x))
             while len(filenames) != n_pieces:
@@ -14,7 +14,7 @@ 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 invesalius.data.imagedata_utils as iu
 import weakref
 from scipy import ndimage
@@ -45,15 +45,25 @@ 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):
+                         from_binary, algorithm, imagedata_resolution, fill_border_holes):
+
+    pad_bottom = (roi.start == 0)
+    pad_top = (roi.stop >= shape[0])
+
+    if fill_border_holes:
+        padding = (1, 1, pad_bottom)
+    else:
+        padding = (0, 0, 0)
+
     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")
+        if fill_border_holes:
+            a_mask = iu.pad_image(mask[roi.start + 1: roi.stop + 1, 1:, 1:], 0, pad_bottom, pad_top)
+        else:
+            a_mask = numpy.array(mask[roi.start + 1: roi.stop + 1, 1:, 1:])
+        image =  converters.to_vtk(a_mask, spacing, roi.start, "AXIAL", padding=padding)
         del a_mask
     else:
         image = numpy.memmap(filename, mode='r', dtype=dtype,
@@ -61,16 +71,21 @@ def create_surface_piece(filename, shape, dtype, mask_filename, mask_shape,
         mask = numpy.memmap(mask_filename, mode='r',
                                  dtype=mask_dtype,
                                  shape=mask_shape)
-        a_image = numpy.array(image[roi])
+        if fill_border_holes:
+            a_image = iu.pad_image(image[roi], numpy.iinfo(image.dtype).min, pad_bottom, pad_top)
+        else:
+            a_image = numpy.array(image[roi])
+        #  if z_iadd:
+            #  a_image[0, 1:-1, 1:-1] = image[0]
+        #  if z_eadd:
+            #  a_image[-1, 1:-1, 1:-1] = image[-1]
         if algorithm == u'InVesalius 3.b2':
-            a_mask = numpy.array(mask[roi.start + 1: roi.stop + 1,
-                                           1:, 1:])
+            a_mask = numpy.array(mask[roi.start + 1: roi.stop + 1, 1:, 1:])
             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")
+            image =  converters.to_vtk(a_image, spacing, roi.start, "AXIAL", padding=padding)
             gauss = vtk.vtkImageGaussianSmooth()
             gauss.SetInputData(image)
@@ -83,8 +98,11 @@ def create_surface_piece(filename, shape, dtype, mask_filename, mask_shape,
             del gauss
             del a_mask
         else:
-            image = converters.to_vtk(a_image, spacing, roi.start,
-                                       "AXIAL")
+            #  if z_iadd:
+                #  origin = -spacing[0], -spacing[1], -spacing[2]
+            #  else:
+                #  origin = 0, -spacing[1], -spacing[2]
+            image = converters.to_vtk(a_image, spacing, roi.start, "AXIAL", padding=padding)
         del a_image
     if imagedata_resolution:
@@ -97,6 +115,11 @@ def create_surface_piece(filename, shape, dtype, mask_filename, mask_shape,
     flip.ReleaseDataFlagOn()
     flip.Update()
+    #  writer = vtk.vtkXMLImageDataWriter()
+    #  writer.SetFileName('/tmp/camboja.vti')
+    #  writer.SetInputData(flip.GetOutput())
+    #  writer.Write()
+
     del image
     image = flip.GetOutput()
     del flip
@@ -224,34 +247,34 @@ def join_process_surface(filenames, algorithm, smooth_iterations, smooth_relaxat
         #  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:
+    #  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 not decimate_reduction:
         print("Decimating", decimate_reduction)
         send_message('Decimating ...')
         decimation = vtk.vtkQuadricDecimation()
@@ -320,9 +343,11 @@ def join_process_surface(filenames, algorithm, smooth_iterations, smooth_relaxat
     #  normals_ref().AddObserver("ProgressEvent", lambda obj,evt:
                     #  UpdateProgress(normals_ref(), _("Creating 3D surface...")))
     normals.SetInputData(polydata)
-    #  normals.SetFeatureAngle(80)
-    #  normals.SplittingOff()
-    #  normals.AutoOrientNormalsOn()
+    normals.SetFeatureAngle(80)
+    normals.SplittingOn()
+    normals.AutoOrientNormalsOn()
+    normals.NonManifoldTraversalOn()
+    normals.ComputeCellNormalsOn()
     #  normals.GetOutput().ReleaseDataFlagOn()
     normals.Update()
     del polydata
@@ -332,22 +357,22 @@ def join_process_surface(filenames, algorithm, smooth_iterations, smooth_relaxat
     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
+    #  # 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()
@@ -1662,7 +1662,10 @@ class SurfaceCreationOptionsPanel(wx.Panel):
                               (combo_quality, 0, flag_button, 0)])
-        # LINES 4 and 5: Checkboxes
+        # LINES 4, 5 and 6: Checkboxes
+        check_box_border_holes = wx.CheckBox(self, -1, _("Fill border holes"))
+        check_box_border_holes.SetValue(False)
+        self.check_box_border_holes = check_box_border_holes
         check_box_holes = wx.CheckBox(self, -1, _("Fill holes"))
         check_box_holes.SetValue(False)
         self.check_box_holes = check_box_holes
@@ -1673,6 +1676,7 @@ class SurfaceCreationOptionsPanel(wx.Panel):
         # Merge all sizers and checkboxes
         sizer = wx.BoxSizer(wx.VERTICAL)
         sizer.Add(fixed_sizer, 0, wx.TOP|wx.RIGHT|wx.LEFT|wx.GROW|wx.EXPAND, 5)
+        sizer.Add(check_box_border_holes, 0, wx.RIGHT|wx.LEFT, 5)
         sizer.Add(check_box_holes, 0, wx.RIGHT|wx.LEFT, 5)
         sizer.Add(check_box_largest, 0, wx.RIGHT|wx.LEFT, 5)
@@ -1687,11 +1691,13 @@ class SurfaceCreationOptionsPanel(wx.Panel):
         mask_index = self.combo_mask.GetSelection()
         surface_name = self.text.GetValue()
         quality = const.SURFACE_QUALITY_LIST[self.combo_quality.GetSelection()]
+        fill_border_holes = self.check_box_border_holes.GetValue()
         fill_holes = self.check_box_holes.GetValue()
         keep_largest = self.check_box_largest.GetValue()
         return {"index": mask_index,
                 "name": surface_name,
                 "quality": quality,
+                "fill_border_holes": fill_border_holes,
                 "fill": fill_holes,
                 "keep_largest": keep_largest,
                 "overwrite": False}
@@ -15,7 +15,7 @@ except ImportError: # if it&#39;s not there locally, try the wxPython lib.
 class Preferences(wx.Dialog):
     def __init__( self, parent, id = ID, title = _("Preferences"), size=wx.DefaultSize,\
-                                pos=wx.DefaultPosition, style=wx.DEFAULT_DIALOG_STYLE):
+                                pos=wx.DefaultPosition, style=wx.DEFAULT_DIALOG_STYLE | wx.RESIZE_BORDER):
         try:
             pre = wx.PreDialog()
@@ -36,14 +36,16 @@ class Preferences(wx.Dialog):
         else:
             self.book = fnb.FlatNotebook(self, wx.ID_ANY, agwStyle=bookStyle)
-        sizer.Add(self.book, 80, wx.EXPAND|wx.ALL)
+        sizer.Add(self.book, 1, wx.EXPAND|wx.ALL)
         self.pnl_viewer2d = Viewer2D(self)
         self.pnl_viewer3d = Viewer3D(self)
+        #  self.pnl_surface = SurfaceCreation(self)
         self.pnl_language = Language(self)
         self.book.AddPage(self.pnl_viewer2d, _("2D Visualization"))
         self.book.AddPage(self.pnl_viewer3d, _("3D Visualization"))
+        #  self.book.AddPage(self.pnl_surface, _("Surface creation"))
         self.book.AddPage(self.pnl_language, _("Language"))
         line = wx.StaticLine(self, -1, size=(20,-1), style=wx.LI_HORIZONTAL)
@@ -59,7 +61,7 @@ class Preferences(wx.Dialog):
         btnsizer.Realize()
-        sizer.Add(btnsizer, 10, wx.GROW|wx.ALIGN_CENTER_VERTICAL|wx.RIGHT|wx.TOP|wx.BOTTOM, 5)
+        sizer.Add(btnsizer, 0, wx.GROW|wx.ALIGN_CENTER_VERTICAL|wx.RIGHT|wx.TOP|wx.BOTTOM, 5)
         self.SetSizer(sizer)
         sizer.Fit(self)
@@ -214,3 +216,21 @@ class Language(wx.Panel):
         locales = self.lg.GetLocalesKey()
         selection = locales.index(language)
         self.cmb_lang.SetSelection(int(selection))
+
+
+
+class SurfaceCreation(wx.Panel):
+    def __init__(self, parent):
+        wx.Panel.__init__(self, parent)
+        self.rb_fill_border = wx.RadioBox(self, -1, "Fill border holes", choices=[_('Yes'), _('No')], style=wx.RA_SPECIFY_COLS | wx.NO_BORDER)
+
+        sizer = wx.BoxSizer(wx.VERTICAL)
+        sizer.Add(self.rb_fill_border)
+
+        self.SetSizerAndFit(sizer)
+
+    def GetSelection(self):
+        return {}
+
+    def LoadSelection(self, values):
+        pass
	@@ -751,3 +751,22 @@ def imgnormalize(data, srange=(0, 255)):		@@ -751,3 +751,22 @@ def imgnormalize(data, srange=(0, 255)):
751	datan = datan.astype(numpy.int16)	751	datan = datan.astype(numpy.int16)
752		752
753	return datan	753	return datan
		754	+
		755	+
		756	+def pad_image(image, pad_value, pad_bottom, pad_top):
		757	+ dz, dy, dx = image.shape
		758	+ z_iadd = 0
		759	+ z_eadd = 0
		760	+ if pad_bottom:
		761	+ z_iadd = 1
		762	+ dz += 1
		763	+ if pad_top:
		764	+ z_eadd = 1
		765	+ dz += 1
		766	+ new_shape = dz, dy + 2, dx + 2
		767	+
		768	+ paded_image = numpy.empty(shape=new_shape, dtype=image.dtype)
		769	+ paded_image[:] = pad_value
		770	+ paded_image[z_iadd: z_iadd + image.shape[0], 1:-1, 1:-1] = image
		771	+
		772	+ return paded_image