diff --git a/environment.yml b/environment.yml index 0b008ab48..07c70e36a 100644 --- a/environment.yml +++ b/environment.yml @@ -1,6 +1,7 @@ channels: - conda-forge - bioconda + - pytorch dependencies: - python==3.11 - cython==3.0.10 diff --git a/invesalius/constants.py b/invesalius/constants.py index c6d2f448a..fec239c0c 100644 --- a/invesalius/constants.py +++ b/invesalius/constants.py @@ -621,6 +621,7 @@ ID_MASK_3D_PREVIEW = wx.NewIdRef() ID_MASK_3D_RELOAD = wx.NewIdRef() ID_MASK_3D_AUTO_RELOAD = wx.NewIdRef() +ID_MASK_3D_EDIT = wx.NewIdRef() ID_GOTO_SLICE = wx.NewIdRef() ID_GOTO_COORD = wx.NewIdRef() @@ -648,6 +649,7 @@ STATE_MEASURE_DENSITY_POLYGON = 1011 STATE_NAVIGATION = 1012 STATE_REGISTRATION = 1013 +STATE_MASK_3D_EDIT = 1014 SLICE_STATE_CROSS = 3006 SLICE_STATE_SCROLL = 3007 @@ -727,6 +729,8 @@ STATE_REGISTRATION: 3, # Override all other states when in navigation mode. STATE_NAVIGATION: 4, + # TODO(Henrique): check whether we need to override the previous states + STATE_MASK_3D_EDIT: 3, } # ------------ Prefereces options key ------------ @@ -1050,3 +1054,9 @@ # Pedal KEYSTROKE_PEDAL_ENABLED = True KEYSTROKE_PEDAL_KEY = wx.WXK_F21 + +# Mask 3D Edit modes + +MASK_3D_EDIT_INCLUDE = 0 +MASK_3D_EDIT_EXCLUDE = 1 +MASK_3D_EDIT_OP_NAME = [_("Include Inside"), _("Exclude Inside")] diff --git a/invesalius/data/polygon_select.py b/invesalius/data/polygon_select.py new file mode 100644 index 000000000..f79704cd1 --- /dev/null +++ b/invesalius/data/polygon_select.py @@ -0,0 +1,71 @@ +# -------------------------------------------------------------------------- +# Software: InVesalius - Software de Reconstrucao 3D de Imagens Medicas +# Copyright: (C) 2001 Centro de Pesquisas Renato Archer +# Homepage: http://www.softwarepublico.gov.br +# Contact: invesalius@cti.gov.br +# License: GNU - GPL 2 (LICENSE.txt/LICENCA.txt) +# -------------------------------------------------------------------------- +# Este programa e software livre; voce pode redistribui-lo e/ou +# modifica-lo sob os termos da Licenca Publica Geral GNU, conforme +# publicada pela Free Software Foundation; de acordo com a versao 2 +# da Licenca. +# +# Este programa eh distribuido na expectativa de ser util, mas SEM +# QUALQUER GARANTIA; sem mesmo a garantia implicita de +# COMERCIALIZACAO ou de ADEQUACAO A QUALQUER PROPOSITO EM +# PARTICULAR. Consulte a Licenca Publica Geral GNU para obter mais +# detalhes. +# -------------------------------------------------------------------------- + +import invesalius.constants as const +from invesalius.gui.widgets.canvas_renderer import CanvasHandlerBase, Polygon + +# import invesalius.session as ses + + +class PolygonSelectCanvas(CanvasHandlerBase): + """ + Inspired on PolygonDensityMeasure, stores and renders a polygon for a canvas + """ + + def __init__(self, colour=(255, 0, 0, 255), interactive=True): + super(PolygonSelectCanvas, self).__init__(None) + self.parent = None + self.children = [] + self.layer = 0 + + self.colour = colour + self.points = [] + + self.complete = False + + self.location = const.SURFACE + self.index = 0 + + self.polygon = Polygon(self, fill=False, closed=False, line_colour=self.colour, is_3d=False) + self.polygon.layer = 1 + self.add_child(self.polygon) + + self.text_box = None + self.interactive = interactive + + def draw_to_canvas(self, gc, canvas): + pass + + def insert_point(self, point): + self.polygon.append_point(point) + + def complete_polygon(self): + if len(self.polygon.points) >= 3: + self.polygon.closed = True + self.complete = True + + def IsComplete(self): + return self.complete + + def SetVisibility(self, value): + self.polygon.visible = value + + def set_interactive(self, value): + self.interactive = bool(value) + self.polygon.interactive = self.interactive diff --git a/invesalius/data/styles_3d.py b/invesalius/data/styles_3d.py index 3205bc884..5c7c86453 100644 --- a/invesalius/data/styles_3d.py +++ b/invesalius/data/styles_3d.py @@ -28,10 +28,24 @@ import invesalius.constants as const import invesalius.project as prj +import invesalius.segmentation.mask_3d_edit as m3e +from invesalius.data.polygon_select import PolygonSelectCanvas from invesalius.pubsub import pub as Publisher PROP_MEASURE = 0.8 +import numpy as np + + +# TODO: find a better place +def vtkarray_to_numpy(m): + nm = np.zeros((4, 4)) + for i in range(4): + for j in range(4): + nm[i, j] = m.GetElement(i, j) + + return nm + class Base3DInteractorStyle(vtkInteractorStyleTrackballCamera): def __init__(self, viewer): @@ -620,6 +634,94 @@ def OnScrollBackward(self, evt, obj): super().OnScrollBackward(evt, obj) +class Mask3DEditorInteractorStyle(DefaultInteractorStyle): + """ + Interactor style for selecting a polygon of interest and performing a mesh edit based on that. + """ + + def __init__(self, viewer): + super().__init__(viewer) + self.viewer = viewer + + # Manages the mask operations + self.mask3deditor = m3e.Mask3DEditor() + + self.picker = vtkCellPicker() + self.picker.PickFromListOn() + + self.has_open_poly = False + self.poly = None + + self.RemoveObservers("LeftButtonPressEvent") + self.AddObserver("LeftButtonPressEvent", self.OnInsertPolygonPoint) + self.AddObserver("RightButtonPressEvent", self.OnInsertPolygon) + + Publisher.subscribe(self.ClearPolygons, "M3E clear polygons") + + def CleanUp(self): + self.RemoveObservers("LeftButtonPressEvent") + self.RemoveObservers("RightButtonPressEvent") + + def ClearPolygons(self): + self.viewer.canvas.draw_list.clear() + self.viewer.UpdateCanvas() + + def OnInsertPolygonPoint(self, obj, evt): + interactor = self.viewer.interactor + interactor.Render() + + mouse_x, mouse_y = self.viewer.interactor.GetEventPosition() + + if not self.has_open_poly: + self.poly = PolygonSelectCanvas() + self.has_open_poly = True + self.viewer.canvas.draw_list.append(self.poly) + + self.poly.insert_point((mouse_x, mouse_y)) + self.viewer.UpdateCanvas() + + def __get_cam_parameters(self): + w, h = self.viewer.GetSize() + self.viewer.ren.Render() + cam = self.viewer.ren.GetActiveCamera() + near, far = cam.GetClippingRange() + + # This flip around the Y axis was done to countereffect the flip that vtk performs + # in volume.py:780. If we do not flip back, what is being displayed is flipped, + # although the actual coordinates are the initial ones, so the cutting gets wrong + # after rotations around y or x. + inv_Y_matrix = np.eye(4) + inv_Y_matrix[1, 1] = -1 + + # Composite transform world to viewport (projection * view) + M = cam.GetCompositeProjectionTransformMatrix(w / float(h), near, far) + M = vtkarray_to_numpy(M) + M = M @ inv_Y_matrix + + MV = cam.GetViewTransformMatrix() + MV = vtkarray_to_numpy(MV) + MV = MV @ inv_Y_matrix + + params = { + "model_to_screen": M, + "model_view": MV, + "resolution": (w, h), + "clipping_range": (near, far), + } + + return params + + def OnInsertPolygon(self, obj, evt): + self.poly.complete_polygon() + self.has_open_poly = False + + self.viewer.UpdateCanvas() + self.viewer.ren.Render() + + Publisher.sendMessage("M3E add polygon", points=self.poly.polygon.points) + Publisher.sendMessage("M3E set camera", params=self.__get_cam_parameters()) + + class Styles: styles = { const.STATE_DEFAULT: DefaultInteractorStyle, @@ -634,6 +736,7 @@ class Styles: const.SLICE_STATE_CROSS: CrossInteractorStyle, const.STATE_NAVIGATION: NavigationInteractorStyle, const.STATE_REGISTRATION: RegistrationInteractorStyle, + const.STATE_MASK_3D_EDIT: Mask3DEditorInteractorStyle, } @classmethod diff --git a/invesalius/gui/task_slice.py b/invesalius/gui/task_slice.py index e9ef6ccb2..3c56edbed 100644 --- a/invesalius/gui/task_slice.py +++ b/invesalius/gui/task_slice.py @@ -758,6 +758,42 @@ def __init__(self, parent): self.gradient_thresh = gradient_thresh self.bind_evt_gradient = True + ## LINE 5 + m3ediv = wx.StaticLine( + self, wx.ID_ANY, wx.DefaultPosition, wx.DefaultSize, wx.LI_HORIZONTAL + ) + + cbox_mask_edit_3d = wx.CheckBox(self, -1, "Edit in 3D") + btn_do_3d_edit = wx.Button(self, -1, _("Cut")) + combo_mask_edit_3d_op = wx.ComboBox( + self, -1, "", choices=const.MASK_3D_EDIT_OP_NAME, style=wx.CB_DROPDOWN | wx.CB_READONLY + ) + combo_mask_edit_3d_op.SetSelection(const.MASK_3D_EDIT_INCLUDE) + + line5 = wx.BoxSizer(wx.HORIZONTAL) + # line5.Add(txt_mesh_edit_3d) + line5.Add(cbox_mask_edit_3d) + line5.Add(btn_do_3d_edit) + line5.Add(combo_mask_edit_3d_op) + + self.cbox_mask_edit_3d = cbox_mask_edit_3d + self.btn_do_3d_edit = btn_do_3d_edit + self.combo_mask_edit_3d_op = combo_mask_edit_3d_op + + ## LINE 6 + cbox_mask_edit_3d_depth = wx.CheckBox(self, -1, "Cut with depth") + spin_mask_edit_3d_depth = wx.SpinCtrlDouble(self, value="1.0", min=0.0, max=1.0, inc=0.05) + btn_clear_3d_poly = wx.Button(self, -1, _("Clear Polygons")) + + line6 = wx.BoxSizer(wx.HORIZONTAL) + line6.Add(cbox_mask_edit_3d_depth) + line6.Add(spin_mask_edit_3d_depth) + line6.Add(btn_clear_3d_poly) + + self.cbox_mask_edit_3d_depth = cbox_mask_edit_3d_depth + self.spin_mask_edit_3d_depth = spin_mask_edit_3d_depth + self.btn_clear_3d_poly = btn_clear_3d_poly + # Add lines into main sizer sizer = wx.BoxSizer(wx.VERTICAL) sizer.AddSpacer(7) @@ -768,6 +804,12 @@ def __init__(self, parent): sizer.Add(text_thresh, 0, wx.GROW | wx.EXPAND | wx.LEFT | wx.RIGHT, 5) sizer.AddSpacer(5) sizer.Add(gradient_thresh, 0, wx.EXPAND | wx.LEFT | wx.RIGHT, 5) + sizer.AddSpacer(5) + sizer.Add(m3ediv, 0, wx.EXPAND | wx.LEFT | wx.RIGHT, 5) + sizer.AddSpacer(3) + sizer.Add(line5, 0, wx.EXPAND | wx.LEFT | wx.RIGHT, 5) + sizer.AddSpacer(3) + sizer.Add(line6, 0, wx.EXPAND | wx.LEFT | wx.RIGHT, 5) sizer.AddSpacer(7) sizer.Fit(self) @@ -782,6 +824,12 @@ def __bind_events_wx(self): self.btn_brush_format.Bind(wx.EVT_MENU, self.OnMenu) self.Bind(grad.EVT_THRESHOLD_CHANGED, self.OnGradientChanged, self.gradient_thresh) self.combo_brush_op.Bind(wx.EVT_COMBOBOX, self.OnComboBrushOp) + self.cbox_mask_edit_3d.Bind(wx.EVT_CHECKBOX, self.OnCheckMaskEdit3D) + self.btn_do_3d_edit.Bind(wx.EVT_BUTTON, self.OnDoMaskEdit3D) + self.combo_mask_edit_3d_op.Bind(wx.EVT_COMBOBOX, self.OnComboMaskEdit3DMode) + self.cbox_mask_edit_3d_depth.Bind(wx.EVT_CHECKBOX, self.OnUseDepthMaskEdit3D) + self.spin_mask_edit_3d_depth.Bind(wx.EVT_SPINCTRLDOUBLE, self.OnSpinDepthMaskEdit3D) + self.btn_clear_3d_poly.Bind(wx.EVT_BUTTON, self.OnClearPolyMaskEdit3D) def __bind_events(self): Publisher.subscribe(self.SetThresholdBounds, "Update threshold limits") @@ -885,6 +933,34 @@ def OnSetUnit(self, evt): self.unit = self.txt_unit.GetLabel() Publisher.sendMessage("Set edition brush unit", unit=self.unit) + def OnCheckMaskEdit3D(self, evt): + style_id = const.STATE_MASK_3D_EDIT + + if self.cbox_mask_edit_3d.GetValue(): + Publisher.sendMessage("Enable style", style=style_id) + Publisher.sendMessage("Enable mask 3D preview") + else: + Publisher.sendMessage("Disable style", style=style_id) + + def OnDoMaskEdit3D(self, evt): + Publisher.sendMessage("M3E apply edit") + + def OnComboMaskEdit3DMode(self, evt): + op_id = evt.GetSelection() + Publisher.sendMessage("M3E set edit mode", mode=op_id) + + def OnUseDepthMaskEdit3D(self, evt): + cb_val = self.cbox_mask_edit_3d_depth.GetValue() + Publisher.sendMessage("M3E use depth", use=cb_val) + self.OnSpinDepthMaskEdit3D(evt) # To pass the set value + + def OnSpinDepthMaskEdit3D(self, evt): + spin_val = self.spin_mask_edit_3d_depth.GetValue() + Publisher.sendMessage("M3E depth value", value=spin_val) + + def OnClearPolyMaskEdit3D(self, evt): + Publisher.sendMessage("M3E clear polygons") + class WatershedTool(EditionTools): def __init__(self, parent): diff --git a/invesalius/segmentation/mask_3d_edit.py b/invesalius/segmentation/mask_3d_edit.py new file mode 100644 index 000000000..8733df756 --- /dev/null +++ b/invesalius/segmentation/mask_3d_edit.py @@ -0,0 +1,206 @@ +# -------------------------------------------------------------------------- +# Software: InVesalius - Software de Reconstrucao 3D de Imagens Medicas +# Copyright: (C) 2001 Centro de Pesquisas Renato Archer +# Homepage: http://www.softwarepublico.gov.br +# Contact: invesalius@cti.gov.br +# License: GNU - GPL 2 (LICENSE.txt/LICENCA.txt) +# -------------------------------------------------------------------------- +# Este programa e software livre; voce pode redistribui-lo e/ou +# modifica-lo sob os termos da Licenca Publica Geral GNU, conforme +# publicada pela Free Software Foundation; de acordo com a versao 2 +# da Licenca. +# +# Este programa eh distribuido na expectativa de ser util, mas SEM +# QUALQUER GARANTIA; sem mesmo a garantia implicita de +# COMERCIALIZACAO ou de ADEQUACAO A QUALQUER PROPOSITO EM +# PARTICULAR. Consulte a Licenca Publica Geral GNU para obter mais +# detalhes. +# -------------------------------------------------------------------------- + +import numpy as np +from skimage.draw import polygon + +import invesalius.constants as const +import invesalius.data.slice_ as slc +from invesalius.pubsub import pub as Publisher +from invesalius_cy.mask_cut import mask_cut, mask_cut_with_depth + + +class Mask3DEditException(Exception): + pass + + +class Mask3DEditor: + def __init__(self): + self.__bind_events() + self.polygons_to_operate = [] + self.spacing = None + self.edit_mode = const.MASK_3D_EDIT_INCLUDE + self.use_depth = False + self.depth_val = None + self.model_to_screen = None + self.model_view = None + self.resolution = None # (w, h) + self.clipping_range = None # (near, far) + + def __bind_events(self): + Publisher.subscribe(self.AddPolygon, "M3E add polygon") + Publisher.subscribe(self.DoMaskEdit, "M3E apply edit") + Publisher.subscribe(self.ClearPolygons, "M3E clear polygons") + Publisher.subscribe(self.SetCamParameters, "M3E set camera") + Publisher.subscribe(self.SetEditMode, "M3E set edit mode") + Publisher.subscribe(self.SetUseDepthForEdit, "M3E use depth") + Publisher.subscribe(self.SetDepthValue, "M3E depth value") + + def AddPolygon(self, points): + """ + Adds polygon to be used in the edit + """ + self.polygons_to_operate.append(points) + + def ClearPolygons(self): + """ + Discards all added polygons + """ + self.polygons_to_operate = [] + + def SetCamParameters(self, params): + """ + Sets the model-view-projection matrix used by the viewer + """ + if "model_to_screen" in params: + self.model_to_screen = params["model_to_screen"] + if "model_view" in params: + self.model_view = params["model_view"] + if "resolution" in params: + self.resolution = params["resolution"] + if "clipping_range" in params: + self.clipping_range = params["clipping_range"] + + def SetEditMode(self, mode): + """ + Sets edit mode to either discard points within or outside + the polygon. + """ + self.edit_mode = mode + + def SetUseDepthForEdit(self, use): + """ + Sets whether to perform a mask cut using depth or through all + """ + self.use_depth = use + + def SetDepthValue(self, value): + self.depth_val = value + + def __create_filter(self): + """ + Based on the polygons and screen resolution, + create a filter for the edit. + """ + w, h = self.resolution + _filter = np.zeros((h, w), dtype="uint8") + + # Include all selected polygons to create the cut + for poly_points in self.polygons_to_operate: + print(f"Poly Points: {poly_points}") + poly = np.array(poly_points) + rr, cc = polygon(poly[:, 1], poly[:, 0], _filter.shape) + _filter[rr, cc] = 1 + + if self.edit_mode == const.MASK_3D_EDIT_INCLUDE: + _filter = 1 - _filter + + return _filter + + # Unoptimized implementation + def __cut_mask(self, mask_data, spacing, _filter): + """ + Cuts an Invesalius Mask with a filter (pixel-wise filter) + mask_data: matrix data without metadata ([1:]) + spacing: the spacing used in the mask matrix + _filter: matrix the same size as the viewer that rasterizes + points to be edited + """ + + dz, dy, dx = mask_data.shape + sx, sy, sz = spacing + + M = self.model_to_screen + h, w = _filter.shape + + for z in range(dz): + for y in range(dy): + for x in range(dx): + # Voxel to world space + p0 = x * sx + p1 = y * sy + p2 = z * sz + p3 = 1.0 + + # _q = M * _p + _q0 = p0 * M[0, 0] + p1 * M[0, 1] + p2 * M[0, 2] + p3 * M[0, 3] + _q1 = p0 * M[1, 0] + p1 * M[1, 1] + p2 * M[1, 2] + p3 * M[1, 3] + _q2 = p0 * M[2, 0] + p1 * M[2, 1] + p2 * M[2, 2] + p3 * M[2, 3] + _q3 = p0 * M[3, 0] + p1 * M[3, 1] + p2 * M[3, 2] + p3 * M[3, 3] + + if _q3 > 0: + q0 = _q0 / _q3 + q1 = _q1 / _q3 + q2 = _q2 / _q3 + + # Normalized coordinates back to pixels + px = (q0 / 2.0 + 0.5) * (w - 1) + py = (q1 / 2.0 + 0.5) * (h - 1) + + if 0 <= px < w and 0 <= py < h: + if ( + _filter[int(py), int(px)] == 1 + ): # NOTE: The lack of round here might be a problem + mask_data[z, y, x] = 0 + + def DoMaskEdit(self): + if len(self.polygons_to_operate) == 0: + return + + s = slc.Slice() + _cur_mask = s.current_mask + + if _cur_mask is None: + raise Mask3DEditException("Attempted editing a non-existent mask") + + _filter = self.__create_filter() + + _prev_mat = _cur_mask.matrix.copy() + + # Unoptimized implementation + # self.__cut_mask(_cur_mask.matrix[1:, 1:, 1:], s.spacing, _filter) + + # Optimized implementation + _mat = _cur_mask.matrix[1:, 1:, 1:].copy() + sx, sy, sz = s.spacing + + print(s.spacing) + print(_mat.shape) + + if self.use_depth: + near, far = self.clipping_range + depth = near + (far - near) * self.depth_val + mask_cut_with_depth( + _mat, sx, sy, sz, depth, _filter, self.model_to_screen, self.model_view, _mat + ) + else: + mask_cut(_mat, sx, sy, sz, _filter, self.model_to_screen, _mat) + + _cur_mask.matrix[1:, 1:, 1:] = _mat + _cur_mask.modified(all_volume=True) + + # Discard all buffers to reupdate view + for ori in ["AXIAL", "CORONAL", "SAGITAL"]: + s.buffer_slices[ori].discard_buffer() + + # Save modification in the history + _cur_mask.save_history(0, "VOLUME", _cur_mask.matrix.copy(), _prev_mat) + + Publisher.sendMessage("Update mask 3D preview") + Publisher.sendMessage("Reload actual slice") diff --git a/invesalius_cy/mask_cut.pyx b/invesalius_cy/mask_cut.pyx new file mode 100644 index 000000000..bcf7d5077 --- /dev/null +++ b/invesalius_cy/mask_cut.pyx @@ -0,0 +1,131 @@ +import numpy as np +cimport numpy as np +cimport cython + +from libc.math cimport floor, ceil, sqrt, fabs, round +from cython.parallel import prange + +ctypedef np.float64_t DTYPEF64_t + +from .cy_my_types cimport image_t, mask_t + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +def mask_cut(np.ndarray[image_t, ndim=3] mask_data, + float sx, float sy, float sz, + np.ndarray[mask_t, ndim=2] filter, + np.ndarray[DTYPEF64_t, ndim=2] M, + np.ndarray[image_t, ndim=3] out): + + cdef int dz = mask_data.shape[0] + cdef int dy = mask_data.shape[1] + cdef int dx = mask_data.shape[2] + + cdef int x + cdef int y + cdef int z + + cdef int h = filter.shape[0] + cdef int w = filter.shape[1] + + cdef DTYPEF64_t px + cdef DTYPEF64_t py + + cdef DTYPEF64_t p0, p1, p2, p3 + cdef DTYPEF64_t _q0, _q1, _q2, _q3 + cdef DTYPEF64_t q0, q1, q2 + + for z in prange(dz, nogil=True): + for y in range(dy): + for x in range(dx): + p0 = (x*sx) + p1 = (y*sy) + p2 = (z*sz) + p3 = 1.0 + + _q0 = p0 * M[0, 0] + p1 * M[0, 1] + p2 * M[0, 2] + p3 * M[0, 3] + _q1 = p0 * M[1, 0] + p1 * M[1, 1] + p2 * M[1, 2] + p3 * M[1, 3] + _q2 = p0 * M[2, 0] + p1 * M[2, 1] + p2 * M[2, 2] + p3 * M[2, 3] + _q3 = p0 * M[3, 0] + p1 * M[3, 1] + p2 * M[3, 2] + p3 * M[3, 3] + + if _q3 > 0: + q0 = _q0/_q3 + q1 = _q1/_q3 + q2 = _q2/_q3 + + px = (q0/2.0 + 0.5) * (w - 1) + py = (q1/2.0 + 0.5) * (h - 1) + + if 0 <= px <= w and 0 <= py <= h: + if filter[(py), (px)]: + out[z, y, x] = 0 + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +def mask_cut_with_depth(np.ndarray[image_t, ndim=3] image, + float sx, float sy, float sz, + float max_depth, + np.ndarray[mask_t, ndim=2] mask, + np.ndarray[DTYPEF64_t, ndim=2] M, + np.ndarray[DTYPEF64_t, ndim=2] MV, + np.ndarray[image_t, ndim=3] out): + + cdef int dz = image.shape[0] + cdef int dy = image.shape[1] + cdef int dx = image.shape[2] + + cdef int x + cdef int y + cdef int z + + cdef int h = mask.shape[0] + cdef int w = mask.shape[1] + + cdef DTYPEF64_t px + cdef DTYPEF64_t py + + cdef DTYPEF64_t p0, p1, p2, p3 + cdef DTYPEF64_t _q0, _q1, _q2, _q3 + cdef DTYPEF64_t q0, q1, q2 + + cdef DTYPEF64_t _c0, _c1, _c2, _c3 + cdef DTYPEF64_t c0, c1, c2 + + cdef DTYPEF64_t dist + + for z in prange(dz, nogil=True): + for y in range(dy): + for x in range(dx): + p0 = (x*sx) + p1 = (y*sy) + p2 = (z*sz) + p3 = 1.0 + + _q0 = p0 * M[0, 0] + p1 * M[0, 1] + p2 * M[0, 2] + p3 * M[0, 3] + _q1 = p0 * M[1, 0] + p1 * M[1, 1] + p2 * M[1, 2] + p3 * M[1, 3] + _q2 = p0 * M[2, 0] + p1 * M[2, 1] + p2 * M[2, 2] + p3 * M[2, 3] + _q3 = p0 * M[3, 0] + p1 * M[3, 1] + p2 * M[3, 2] + p3 * M[3, 3] + + if _q3 > 0: + q0 = _q0/_q3 + q1 = _q1/_q3 + q2 = _q2/_q3 + + _c0 = p0 * MV[0, 0] + p1 * MV[0, 1] + p2 * MV[0, 2] + p3 * MV[0, 3] + _c1 = p0 * MV[1, 0] + p1 * MV[1, 1] + p2 * MV[1, 2] + p3 * MV[1, 3] + _c2 = p0 * MV[2, 0] + p1 * MV[2, 1] + p2 * MV[2, 2] + p3 * MV[2, 3] + _c3 = p0 * MV[3, 0] + p1 * MV[3, 1] + p2 * MV[3, 2] + p3 * MV[3, 3] + + c0 = _c0/_c3 + c1 = _c1/_c3 + c2 = _c2/_c3 + + dist = sqrt(c0*c0 + c1*c1 + c2*c2) + + if dist <= max_depth: + px = (q0/2.0 + 0.5) * (w - 1) + py = (q1/2.0 + 0.5) * (h - 1) + + if 0 <= px <= w and 0 <= py <= h: + if mask[round(py), round(px)]: + out[z, y, x] = 0 \ No newline at end of file diff --git a/setup.py b/setup.py index b183b7059..59780e67e 100644 --- a/setup.py +++ b/setup.py @@ -1,4 +1,3 @@ -import setuptools import logging import os import pathlib @@ -6,7 +5,8 @@ import sys import numpy -from Cython.Build import cythonize, build_ext +import setuptools +from Cython.Build import build_ext, cythonize if sys.platform == "darwin": unix_copt = ["-Xpreprocessor", "-fopenmp", "-lomp"] @@ -62,9 +62,7 @@ def run(self): plugin_folder = plugins_folder.joinpath(p) self.announce("Compiling plugin: {}".format(p)) os.chdir(plugin_folder) - subprocess.check_call( - [sys.executable, "setup.py", "build_ext", "--inplace"] - ) + subprocess.check_call([sys.executable, "setup.py", "build_ext", "--inplace"]) os.chdir(inv_folder) @@ -97,6 +95,11 @@ def run(self): ["invesalius_cy/cy_mesh.pyx"], language="c++", ), + setuptools.Extension( + "invesalius_cy.mask_cut", + ["invesalius_cy/mask_cut.pyx"], + language="c++", + ), ] ), )