Casmoothing implemented in cython (#61)

Implements casmoothing in Cython. - Implements in cython - Better handling of border vertices

Casmoothing implemented in cython (#61)
Implements casmoothing in Cython. - Implements in cython - Better handling of border vertices
Thiago Franco de Moraes · GitHub
1 parent 339908de
Showing 6 changed files with 516 additions and 9 deletions Show diff stats
.gitignore
invesalius/data/cy_mesh.pyx
invesalius/data/cy_my_types.pxd
invesalius/data/surface.py
invesalius/gui/dialogs.py
setup.py
@@ -24,3 +24,9 @@ tags
 build
 *.patch
 *.tgz
+
+*.pyd
+*.cpp
+*.diff
+
+*.directory
@@ -0,0 +1,466 @@
+# distutils: language = c++
+#cython: boundscheck=False
+#cython: wraparound=False
+#cython: initializedcheck=False
+#cython: cdivision=True
+#cython: nonecheck=False
+
+import sys
+import time
+cimport numpy as np
+
+from libc.math cimport sin, cos, acos, exp, sqrt, fabs, M_PI
+from libc.stdlib cimport abs as cabs
+from cython.operator cimport dereference as deref, preincrement as inc
+from libcpp.map cimport map
+from libcpp.unordered_map cimport unordered_map
+from libcpp.set cimport set
+from libcpp.vector cimport vector
+from libcpp.pair cimport pair
+from libcpp cimport bool
+from libcpp.deque cimport deque as cdeque
+from cython.parallel import prange
+
+from cy_my_types cimport vertex_t, normal_t, vertex_id_t
+
+import numpy as np
+import vtk
+
+from vtk.util import numpy_support
+
+ctypedef float weight_t
+
+cdef struct Point:
+    vertex_t x
+    vertex_t y
+    vertex_t z
+
+ctypedef pair[vertex_id_t, vertex_id_t] key
+
+
+cdef class Mesh:
+    cdef vertex_t[:, :] vertices
+    cdef vertex_id_t[:, :] faces
+    cdef normal_t[:, :] normals
+
+    cdef unordered_map[int, vector[vertex_id_t]] map_vface
+    cdef unordered_map[vertex_id_t, int] border_vertices
+
+    cdef bool _initialized
+
+    def __cinit__(self, pd=None, other=None):
+        cdef int i
+        cdef map[key, int] edge_nfaces
+        cdef map[key, int].iterator it
+        if pd:
+            self._initialized = True
+            _vertices = numpy_support.vtk_to_numpy(pd.GetPoints().GetData())
+            _vertices.shape = -1, 3
+
+            _faces = numpy_support.vtk_to_numpy(pd.GetPolys().GetData())
+            _faces.shape = -1, 4
+
+            _normals = numpy_support.vtk_to_numpy(pd.GetCellData().GetArray("Normals"))
+            _normals.shape = -1, 3
+
+            self.vertices = _vertices
+            self.faces = _faces
+            self.normals = _normals
+
+            for i in xrange(_faces.shape[0]):
+                self.map_vface[self.faces[i, 1]].push_back(i)
+                self.map_vface[self.faces[i, 2]].push_back(i)
+                self.map_vface[self.faces[i, 3]].push_back(i)
+
+                edge_nfaces[key(min(self.faces[i, 1], self.faces[i, 2]), max(self.faces[i, 1], self.faces[i, 2]))] += 1
+                edge_nfaces[key(min(self.faces[i, 2], self.faces[i, 3]), max(self.faces[i, 2], self.faces[i, 3]))] += 1
+                edge_nfaces[key(min(self.faces[i, 1], self.faces[i, 3]), max(self.faces[i, 1], self.faces[i, 3]))] += 1
+
+            it = edge_nfaces.begin()
+
+            while it != edge_nfaces.end():
+                if deref(it).second == 1:
+                    self.border_vertices[deref(it).first.first] = 1
+                    self.border_vertices[deref(it).first.second] = 1
+
+                inc(it)
+
+        elif other:
+            _other = <Mesh>other
+            self._initialized = True
+            self.vertices = _other.vertices.copy()
+            self.faces = _other.faces.copy()
+            self.normals = _other.normals.copy()
+            self.map_vface = unordered_map[int, vector[vertex_id_t]](_other.map_vface)
+            self.border_vertices = unordered_map[vertex_id_t, int](_other.border_vertices)
+        else:
+            self._initialized = False
+
+    cdef void copy_to(self, Mesh other):
+        """
+        Copies self content to other.
+        """
+        if self._initialized:
+            other.vertices[:] = self.vertices
+            other.faces[:] = self.faces
+            other.normals[:] = self.normals
+            other.map_vface = unordered_map[int, vector[vertex_id_t]](self.map_vface)
+            other.border_vertices = unordered_map[vertex_id_t, int](self.border_vertices)
+        else:
+            other.vertices = self.vertices.copy()
+            other.faces = self.faces.copy()
+            other.normals = self.normals.copy()
+
+            other.map_vface = self.map_vface
+            other.border_vertices = self.border_vertices
+
+    def to_vtk(self):
+        """
+        Converts Mesh to vtkPolyData.
+        """
+        vertices = np.asarray(self.vertices)
+        faces = np.asarray(self.faces)
+        normals = np.asarray(self.normals)
+
+        points = vtk.vtkPoints()
+        points.SetData(numpy_support.numpy_to_vtk(vertices))
+
+        id_triangles = numpy_support.numpy_to_vtkIdTypeArray(faces)
+        triangles = vtk.vtkCellArray()
+        triangles.SetCells(faces.shape[0], id_triangles)
+
+        pd = vtk.vtkPolyData()
+        pd.SetPoints(points)
+        pd.SetPolys(triangles)
+
+        return pd
+
+    cdef vector[vertex_id_t]* get_faces_by_vertex(self, int v_id) nogil:
+        """
+        Returns the faces whose vertex `v_id' is part.
+        """
+        return &self.map_vface[v_id]
+
+    cdef set[vertex_id_t]* get_ring1(self, vertex_id_t v_id) nogil:
+        """
+        Returns the ring1 of vertex `v_id'
+        """
+        cdef vertex_id_t f_id
+        cdef set[vertex_id_t]* ring1 = new set[vertex_id_t]()
+        cdef vector[vertex_id_t].iterator it = self.map_vface[v_id].begin()
+
+        while it != self.map_vface[v_id].end():
+            f_id = deref(it)
+            inc(it)
+            if self.faces[f_id, 1] != v_id:
+                ring1.insert(self.faces[f_id, 1])
+            if self.faces[f_id, 2] != v_id:
+                ring1.insert(self.faces[f_id, 2])
+            if self.faces[f_id, 3] != v_id:
+                ring1.insert(self.faces[f_id, 3])
+
+        return ring1
+
+    cdef bool is_border(self, vertex_id_t v_id) nogil:
+        """
+        Check if vertex `v_id' is a vertex border.
+        """
+        return self.border_vertices.find(v_id) != self.border_vertices.end()
+
+    cdef vector[vertex_id_t]* get_near_vertices_to_v(self, vertex_id_t v_id, float dmax) nogil:
+        """
+        Returns all vertices with distance at most `d' to the vertex `v_id'
+
+        Params:
+            v_id: id of the vertex
+            dmax: the maximun distance.
+        """
+        cdef vector[vertex_id_t]* idfaces
+        cdef vector[vertex_id_t]* near_vertices = new vector[vertex_id_t]()
+
+        cdef cdeque[vertex_id_t] to_visit
+        cdef unordered_map[vertex_id_t, bool] status_v
+        cdef unordered_map[vertex_id_t, bool] status_f
+
+        cdef vertex_t *vip
+        cdef vertex_t *vjp
+
+        cdef float distance
+        cdef int nf, nid, j
+        cdef vertex_id_t f_id, vj
+
+        vip = &self.vertices[v_id, 0]
+        to_visit.push_back(v_id)
+        while(not to_visit.empty()):
+            v_id = to_visit.front()
+            to_visit.pop_front()
+
+            status_v[v_id] = True
+
+            idfaces = self.get_faces_by_vertex(v_id)
+            nf = idfaces.size()
+
+            for nid in xrange(nf):
+                f_id = deref(idfaces)[nid]
+                if status_f.find(f_id) == status_f.end():
+                    status_f[f_id] = True
+
+                    for j in xrange(3):
+                        vj = self.faces[f_id, j+1]
+                        if status_v.find(vj) == status_v.end():
+                            status_v[vj] = True
+                            vjp = &self.vertices[vj, 0]
+                            distance = sqrt((vip[0] - vjp[0]) * (vip[0] - vjp[0]) \
+                                          + (vip[1] - vjp[1]) * (vip[1] - vjp[1]) \
+                                          + (vip[2] - vjp[2]) * (vip[2] - vjp[2]))
+                            if distance <= dmax:
+                                near_vertices.push_back(vj)
+                                to_visit.push_back(vj)
+
+        return near_vertices
+
+
+cdef vector[weight_t]* calc_artifacts_weight(Mesh mesh, vector[vertex_id_t]& vertices_staircase, float tmax, float bmin) nogil:
+    """
+    Calculate the artifact weight based on distance of each vertex to its
+    nearest staircase artifact vertex.
+
+    Params:
+        mesh: Mesh
+        vertices_staircase: the identified staircase artifact vertices
+        tmax: max distance the vertex must be to its nearest artifact vertex
+              to considered to calculate the weight
+        bmin: The minimun weight.
+    """
+    cdef int vi_id, vj_id, nnv, n_ids, i, j
+    cdef vector[vertex_id_t]* near_vertices
+    cdef weight_t value
+    cdef float d
+    n_ids = vertices_staircase.size()
+
+    cdef vertex_t* vi
+    cdef vertex_t* vj
+    cdef size_t msize
+
+    msize = mesh.vertices.shape[0]
+    cdef vector[weight_t]* weights = new vector[weight_t](msize)
+    weights.assign(msize, bmin)
+
+    for i in prange(n_ids, nogil=True):
+        vi_id = vertices_staircase[i]
+        deref(weights)[vi_id] = 1.0
+
+        vi = &mesh.vertices[vi_id, 0]
+        near_vertices = mesh.get_near_vertices_to_v(vi_id, tmax)
+        nnv = near_vertices.size()
+
+        for j in xrange(nnv):
+            vj_id = deref(near_vertices)[j]
+            vj = &mesh.vertices[vj_id, 0]
+
+            d = sqrt((vi[0] - vj[0]) * (vi[0] - vj[0])\
+                   + (vi[1] - vj[1]) * (vi[1] - vj[1])\
+                   + (vi[2] - vj[2]) * (vi[2] - vj[2]))
+            value = (1.0 - d/tmax) * (1.0 - bmin) + bmin
+
+            if value > deref(weights)[vj_id]:
+                deref(weights)[vj_id] = value
+
+        del near_vertices
+
+    #  for i in xrange(msize):
+        #  if mesh.is_border(i):
+            #  deref(weights)[i] = 0.0
+
+    #  cdef vertex_id_t v0, v1, v2
+    #  for i in xrange(mesh.faces.shape[0]):
+        #  for j in xrange(1, 4):
+            #  v0 = mesh.faces[i, j]
+            #  vi = &mesh.vertices[v0, 0]
+            #  if mesh.is_border(v0):
+                #  deref(weights)[v0] = 0.0
+                #  v1 = mesh.faces[i, (j + 1) % 3 + 1]
+                #  if mesh.is_border(v1):
+                    #  vi = &mesh.vertices[v1, 0]
+                    #  deref(weights)[v0] = 0.0
+
+    return weights
+
+
+cdef inline Point calc_d(Mesh mesh, vertex_id_t v_id) nogil:
+    cdef Point D
+    cdef int nf, f_id, nid
+    cdef float n=0
+    cdef int i
+    cdef vertex_t* vi
+    cdef vertex_t* vj
+    cdef set[vertex_id_t]* vertices
+    cdef set[vertex_id_t].iterator it
+    cdef vertex_id_t vj_id
+
+    D.x = 0.0
+    D.y = 0.0
+    D.z = 0.0
+
+    vertices = mesh.get_ring1(v_id)
+    vi = &mesh.vertices[v_id, 0]
+
+    if mesh.is_border(v_id):
+        it = vertices.begin()
+        while it != vertices.end():
+            vj_id = deref(it)
+            if mesh.is_border(vj_id):
+                vj = &mesh.vertices[vj_id, 0]
+
+                D.x = D.x + (vi[0] - vj[0])
+                D.y = D.y + (vi[1] - vj[1])
+                D.z = D.z + (vi[2] - vj[2])
+                n += 1.0
+
+            inc(it)
+    else:
+        it = vertices.begin()
+        while it != vertices.end():
+            vj_id = deref(it)
+            vj = &mesh.vertices[vj_id, 0]
+
+            D.x = D.x + (vi[0] - vj[0])
+            D.y = D.y + (vi[1] - vj[1])
+            D.z = D.z + (vi[2] - vj[2])
+            n += 1.0
+
+            inc(it)
+
+    del vertices
+
+    D.x = D.x / n
+    D.y = D.y / n
+    D.z = D.z / n
+    return D
+
+cdef vector[vertex_id_t]* find_staircase_artifacts(Mesh mesh, double[3] stack_orientation, double T) nogil:
+    """
+    This function is used to find vertices at staircase artifacts, which are
+    those vertices whose incident faces' orientation differences are
+    greater than T.
+
+    Params:
+        mesh: Mesh
+        stack_orientation: orientation of slice stacking
+        T: Min angle (between vertex faces and stack_orientation) to consider a
+           vertex a staircase artifact.
+    """
+    cdef int nv, nf, f_id, v_id
+    cdef double of_z, of_y, of_x, min_z, max_z, min_y, max_y, min_x, max_x;
+    cdef vector[vertex_id_t]* f_ids
+    cdef normal_t* normal
+
+    cdef vector[vertex_id_t]* output = new vector[vertex_id_t]()
+    cdef int i
+
+    nv = mesh.vertices.shape[0]
+
+    for v_id in xrange(nv):
+        max_z = -10000
+        min_z = 10000
+        max_y = -10000
+        min_y = 10000
+        max_x = -10000
+        min_x = 10000
+
+        f_ids = mesh.get_faces_by_vertex(v_id)
+        nf = deref(f_ids).size()
+
+        for i in xrange(nf):
+            f_id = deref(f_ids)[i]
+            normal = &mesh.normals[f_id, 0]
+
+            of_z = 1 - fabs(normal[0]*stack_orientation[0] + normal[1]*stack_orientation[1] + normal[2]*stack_orientation[2]);
+            of_y = 1 - fabs(normal[0]*0 + normal[1]*1 + normal[2]*0);
+            of_x = 1 - fabs(normal[0]*1 + normal[1]*0 + normal[2]*0);
+
+            if (of_z > max_z):
+                max_z = of_z
+
+            if (of_z < min_z):
+                min_z = of_z
+
+            if (of_y > max_y):
+                max_y = of_y
+
+            if (of_y < min_y):
+                min_y = of_y
+
+            if (of_x > max_x):
+                max_x = of_x
+
+            if (of_x < min_x):
+                min_x = of_x
+
+
+            if ((fabs(max_z - min_z) >= T) or (fabs(max_y - min_y) >= T) or (fabs(max_x - min_x) >= T)):
+                output.push_back(v_id)
+                break
+    return output
+
+
+cdef void taubin_smooth(Mesh mesh, vector[weight_t]& weights, float l, float m, int steps):
+    """
+    Implementation of Taubin's smooth algorithm described in the paper "A
+    Signal Processing Approach To Fair Surface Design". His benefeat is it
+    avoids surface shrinking.
+    """
+    cdef int s, i, nvertices
+    nvertices = mesh.vertices.shape[0]
+    cdef vector[Point] D = vector[Point](nvertices)
+    cdef vertex_t* vi
+    for s in xrange(steps):
+        for i in prange(nvertices, nogil=True):
+            D[i] = calc_d(mesh, i)
+
+        for i in prange(nvertices, nogil=True):
+            mesh.vertices[i, 0] += weights[i]*l*D[i].x;
+            mesh.vertices[i, 1] += weights[i]*l*D[i].y;
+            mesh.vertices[i, 2] += weights[i]*l*D[i].z;
+
+        for i in prange(nvertices, nogil=True):
+            D[i] = calc_d(mesh, i)
+
+        for i in prange(nvertices, nogil=True):
+            mesh.vertices[i, 0] += weights[i]*m*D[i].x;
+            mesh.vertices[i, 1] += weights[i]*m*D[i].y;
+            mesh.vertices[i, 2] += weights[i]*m*D[i].z;
+
+
+def ca_smoothing(Mesh mesh, double T, double tmax, double bmin, int n_iters):
+    """
+    This is a implementation of the paper "Context-aware mesh smoothing for
+    biomedical applications". It can be used to smooth meshes generated by
+    binary images to remove its staircase artifacts and keep the fine features.
+
+    Params:
+        mesh: Mesh
+        T: Min angle (between vertex faces and stack_orientation) to consider a
+           vertex a staircase artifact
+        tmax: max distance the vertex must be to its nearest artifact vertex
+              to considered to calculate the weight
+        bmin: The minimun weight
+        n_iters: Number of iterations.
+    """
+    cdef double[3] stack_orientation = [0.0, 0.0, 1.0]
+
+    t0 = time.time()
+    cdef vector[vertex_id_t]* vertices_staircase =  find_staircase_artifacts(mesh, stack_orientation, T)
+    print "vertices staircase", time.time() - t0
+
+    t0 = time.time()
+    cdef vector[weight_t]* weights = calc_artifacts_weight(mesh, deref(vertices_staircase), tmax, bmin)
+    print "Weights", time.time() - t0
+
+    del vertices_staircase
+
+    t0 = time.time()
+    taubin_smooth(mesh, deref(weights), 0.5, -0.53, n_iters)
+    print "taubin", time.time() - t0
+
+    del weights
@@ -10,3 +10,7 @@ ctypedef fused image_t:
     np.uint8_t
 ctypedef np.uint8_t mask_t
+
+ctypedef np.float32_t vertex_t
+ctypedef np.float32_t normal_t
+ctypedef np.int64_t vertex_id_t
@@ -42,6 +42,7 @@ try:
 except ImportError:
     import data.ca_smoothing as ca_smoothing
+import cy_mesh
 # TODO: Verificar ReleaseDataFlagOn and SetSource 
 class Surface():
@@ -564,16 +565,28 @@ class SurfaceManager():
             #  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'])
+            #  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()
+            w = vtk.vtkPLYWriter()
+            w.SetInputData(polydata)
+            w.SetFileName('/tmp/ca_smoothing_inv.ply')
+            w.Write()
         else:
             #smoother = vtk.vtkWindowedSincPolyDataFilter()
@@ -1276,7 +1276,7 @@ class CAOptions(wx.Panel):
                                          max_val=100.0, increment=0.1,
                                          digits=2)
-        self.min_weight = floatspin.FloatSpin(self, -1, value=0.2, min_val=0.0,
+        self.min_weight = floatspin.FloatSpin(self, -1, value=0.5, min_val=0.0,
                                          max_val=1.0, increment=0.1,
                                          digits=1)
@@ -29,6 +29,12 @@ if sys.platform == &#39;linux2&#39;:
                        Extension("invesalius.data.floodfill", ["invesalius/data/floodfill.pyx"],
                                  include_dirs=[numpy.get_include()],
                                  language='c++',),
+
+                       Extension("invesalius.data.cy_mesh", ["invesalius/data/cy_mesh.pyx"],
+                                 include_dirs=[numpy.get_include()],
+                                 extra_compile_args=['-fopenmp', '-std=c++11'],
+                                 extra_link_args=['-fopenmp', '-std=c++11'],
+                                 language='c++',),
                        ])
          )
@@ -50,6 +56,11 @@ elif sys.platform == &#39;win32&#39;:
                                  Extension("invesalius.data.floodfill", ["invesalius/data/floodfill.pyx"],
                                            include_dirs=[numpy.get_include()],
                                            language='c++',),
+
+                                 Extension("invesalius.data.cy_mesh", ["invesalius/data/cy_mesh.pyx"],
+                                           include_dirs=[numpy.get_include()],
+                                           extra_compile_args=['/openmp',],
+                                           language='c++',),
                                  ])
     )
@@ -75,5 +86,12 @@ else:
                                  Extension("invesalius.data.floodfill", ["invesalius/data/floodfill.pyx"],
                                            include_dirs=[numpy.get_include()],
                                            language='c++',),
+
+                                 Extension("invesalius.data.cy_mesh", ["invesalius/data/cy_mesh.pyx"],
+                                           include_dirs=[numpy.get_include()],
+                                           extra_compile_args=['-fopenmp', '-std=c++11'],
+                                           extra_link_args=['-fopenmp', '-std=c++11'],
+                                           language='c++',),
+
                                  ])
     )
	@@ -24,3 +24,9 @@ tags		@@ -24,3 +24,9 @@ tags
24	build	24	build
25	*.patch	25	*.patch
26	*.tgz	26	*.tgz
		27	+
		28	+*.pyd
		29	+*.cpp
		30	+*.diff
		31	+
		32	+*.directory
@@ -0,0 +1,466 @@		@@ -0,0 +1,466 @@
	1	+# distutils: language = c++
	2	+#cython: boundscheck=False
	3	+#cython: wraparound=False
	4	+#cython: initializedcheck=False
	5	+#cython: cdivision=True
	6	+#cython: nonecheck=False
	7	+
	8	+import sys
	9	+import time
	10	+cimport numpy as np
	11	+
	12	+from libc.math cimport sin, cos, acos, exp, sqrt, fabs, M_PI
	13	+from libc.stdlib cimport abs as cabs
	14	+from cython.operator cimport dereference as deref, preincrement as inc
	15	+from libcpp.map cimport map
	16	+from libcpp.unordered_map cimport unordered_map
	17	+from libcpp.set cimport set
	18	+from libcpp.vector cimport vector
	19	+from libcpp.pair cimport pair
	20	+from libcpp cimport bool
	21	+from libcpp.deque cimport deque as cdeque
	22	+from cython.parallel import prange
	23	+
	24	+from cy_my_types cimport vertex_t, normal_t, vertex_id_t
	25	+
	26	+import numpy as np
	27	+import vtk
	28	+
	29	+from vtk.util import numpy_support
	30	+
	31	+ctypedef float weight_t
	32	+
	33	+cdef struct Point:
	34	+ vertex_t x
	35	+ vertex_t y
	36	+ vertex_t z
	37	+
	38	+ctypedef pair[vertex_id_t, vertex_id_t] key
	39	+
	40	+
	41	+cdef class Mesh:
	42	+ cdef vertex_t[:, :] vertices
	43	+ cdef vertex_id_t[:, :] faces
	44	+ cdef normal_t[:, :] normals
	45	+
	46	+ cdef unordered_map[int, vector[vertex_id_t]] map_vface
	47	+ cdef unordered_map[vertex_id_t, int] border_vertices
	48	+
	49	+ cdef bool _initialized
	50	+
	51	+ def __cinit__(self, pd=None, other=None):
	52	+ cdef int i
	53	+ cdef map[key, int] edge_nfaces
	54	+ cdef map[key, int].iterator it
	55	+ if pd:
	56	+ self._initialized = True
	57	+ _vertices = numpy_support.vtk_to_numpy(pd.GetPoints().GetData())
	58	+ _vertices.shape = -1, 3
	59	+
	60	+ _faces = numpy_support.vtk_to_numpy(pd.GetPolys().GetData())
	61	+ _faces.shape = -1, 4
	62	+
	63	+ _normals = numpy_support.vtk_to_numpy(pd.GetCellData().GetArray("Normals"))
	64	+ _normals.shape = -1, 3
	65	+
	66	+ self.vertices = _vertices
	67	+ self.faces = _faces
	68	+ self.normals = _normals
	69	+
	70	+ for i in xrange(_faces.shape[0]):
	71	+ self.map_vface[self.faces[i, 1]].push_back(i)
	72	+ self.map_vface[self.faces[i, 2]].push_back(i)
	73	+ self.map_vface[self.faces[i, 3]].push_back(i)
	74	+
	75	+ edge_nfaces[key(min(self.faces[i, 1], self.faces[i, 2]), max(self.faces[i, 1], self.faces[i, 2]))] += 1
	76	+ edge_nfaces[key(min(self.faces[i, 2], self.faces[i, 3]), max(self.faces[i, 2], self.faces[i, 3]))] += 1
	77	+ edge_nfaces[key(min(self.faces[i, 1], self.faces[i, 3]), max(self.faces[i, 1], self.faces[i, 3]))] += 1
	78	+
	79	+ it = edge_nfaces.begin()
	80	+
	81	+ while it != edge_nfaces.end():
	82	+ if deref(it).second == 1:
	83	+ self.border_vertices[deref(it).first.first] = 1
	84	+ self.border_vertices[deref(it).first.second] = 1
	85	+
	86	+ inc(it)
	87	+
	88	+ elif other:
	89	+ _other = <Mesh>other
	90	+ self._initialized = True
	91	+ self.vertices = _other.vertices.copy()
	92	+ self.faces = _other.faces.copy()
	93	+ self.normals = _other.normals.copy()
	94	+ self.map_vface = unordered_map[int, vector[vertex_id_t]](_other.map_vface)
	95	+ self.border_vertices = unordered_map[vertex_id_t, int](_other.border_vertices)
	96	+ else:
	97	+ self._initialized = False
	98	+
	99	+ cdef void copy_to(self, Mesh other):
	100	+ """
	101	+ Copies self content to other.
	102	+ """
	103	+ if self._initialized:
	104	+ other.vertices[:] = self.vertices
	105	+ other.faces[:] = self.faces
	106	+ other.normals[:] = self.normals
	107	+ other.map_vface = unordered_map[int, vector[vertex_id_t]](self.map_vface)
	108	+ other.border_vertices = unordered_map[vertex_id_t, int](self.border_vertices)
	109	+ else:
	110	+ other.vertices = self.vertices.copy()
	111	+ other.faces = self.faces.copy()
	112	+ other.normals = self.normals.copy()
	113	+
	114	+ other.map_vface = self.map_vface
	115	+ other.border_vertices = self.border_vertices
	116	+
	117	+ def to_vtk(self):
	118	+ """
	119	+ Converts Mesh to vtkPolyData.
	120	+ """
	121	+ vertices = np.asarray(self.vertices)
	122	+ faces = np.asarray(self.faces)
	123	+ normals = np.asarray(self.normals)
	124	+
	125	+ points = vtk.vtkPoints()
	126	+ points.SetData(numpy_support.numpy_to_vtk(vertices))
	127	+
	128	+ id_triangles = numpy_support.numpy_to_vtkIdTypeArray(faces)
	129	+ triangles = vtk.vtkCellArray()
	130	+ triangles.SetCells(faces.shape[0], id_triangles)
	131	+
	132	+ pd = vtk.vtkPolyData()
	133	+ pd.SetPoints(points)
	134	+ pd.SetPolys(triangles)
	135	+
	136	+ return pd
	137	+
	138	+ cdef vector[vertex_id_t]* get_faces_by_vertex(self, int v_id) nogil:
	139	+ """
	140	+ Returns the faces whose vertex `v_id' is part.
	141	+ """
	142	+ return &self.map_vface[v_id]
	143	+
	144	+ cdef set[vertex_id_t]* get_ring1(self, vertex_id_t v_id) nogil:
	145	+ """
	146	+ Returns the ring1 of vertex `v_id'
	147	+ """
	148	+ cdef vertex_id_t f_id
	149	+ cdef set[vertex_id_t]* ring1 = new set[vertex_id_t]()
	150	+ cdef vector[vertex_id_t].iterator it = self.map_vface[v_id].begin()
	151	+
	152	+ while it != self.map_vface[v_id].end():
	153	+ f_id = deref(it)
	154	+ inc(it)
	155	+ if self.faces[f_id, 1] != v_id:
	156	+ ring1.insert(self.faces[f_id, 1])
	157	+ if self.faces[f_id, 2] != v_id:
	158	+ ring1.insert(self.faces[f_id, 2])
	159	+ if self.faces[f_id, 3] != v_id:
	160	+ ring1.insert(self.faces[f_id, 3])
	161	+
	162	+ return ring1
	163	+
	164	+ cdef bool is_border(self, vertex_id_t v_id) nogil:
	165	+ """
	166	+ Check if vertex `v_id' is a vertex border.
	167	+ """
	168	+ return self.border_vertices.find(v_id) != self.border_vertices.end()
	169	+
	170	+ cdef vector[vertex_id_t]* get_near_vertices_to_v(self, vertex_id_t v_id, float dmax) nogil:
	171	+ """
	172	+ Returns all vertices with distance at most `d' to the vertex `v_id'
	173	+
	174	+ Params:
	175	+ v_id: id of the vertex
	176	+ dmax: the maximun distance.
	177	+ """
	178	+ cdef vector[vertex_id_t]* idfaces
	179	+ cdef vector[vertex_id_t]* near_vertices = new vector[vertex_id_t]()
	180	+
	181	+ cdef cdeque[vertex_id_t] to_visit
	182	+ cdef unordered_map[vertex_id_t, bool] status_v
	183	+ cdef unordered_map[vertex_id_t, bool] status_f
	184	+
	185	+ cdef vertex_t *vip
	186	+ cdef vertex_t *vjp
	187	+
	188	+ cdef float distance
	189	+ cdef int nf, nid, j
	190	+ cdef vertex_id_t f_id, vj
	191	+
	192	+ vip = &self.vertices[v_id, 0]
	193	+ to_visit.push_back(v_id)
	194	+ while(not to_visit.empty()):
	195	+ v_id = to_visit.front()
	196	+ to_visit.pop_front()
	197	+
	198	+ status_v[v_id] = True
	199	+
	200	+ idfaces = self.get_faces_by_vertex(v_id)
	201	+ nf = idfaces.size()
	202	+
	203	+ for nid in xrange(nf):
	204	+ f_id = deref(idfaces)[nid]
	205	+ if status_f.find(f_id) == status_f.end():
	206	+ status_f[f_id] = True
	207	+
	208	+ for j in xrange(3):
	209	+ vj = self.faces[f_id, j+1]
	210	+ if status_v.find(vj) == status_v.end():
	211	+ status_v[vj] = True
	212	+ vjp = &self.vertices[vj, 0]
	213	+ distance = sqrt((vip[0] - vjp[0]) * (vip[0] - vjp[0]) \
	214	+ + (vip[1] - vjp[1]) * (vip[1] - vjp[1]) \
	215	+ + (vip[2] - vjp[2]) * (vip[2] - vjp[2]))
	216	+ if distance <= dmax:
	217	+ near_vertices.push_back(vj)
	218	+ to_visit.push_back(vj)
	219	+
	220	+ return near_vertices
	221	+
	222	+
	223	+cdef vector[weight_t]* calc_artifacts_weight(Mesh mesh, vector[vertex_id_t]& vertices_staircase, float tmax, float bmin) nogil:
	224	+ """
	225	+ Calculate the artifact weight based on distance of each vertex to its
	226	+ nearest staircase artifact vertex.
	227	+
	228	+ Params:
	229	+ mesh: Mesh
	230	+ vertices_staircase: the identified staircase artifact vertices
	231	+ tmax: max distance the vertex must be to its nearest artifact vertex
	232	+ to considered to calculate the weight
	233	+ bmin: The minimun weight.
	234	+ """
	235	+ cdef int vi_id, vj_id, nnv, n_ids, i, j
	236	+ cdef vector[vertex_id_t]* near_vertices
	237	+ cdef weight_t value
	238	+ cdef float d
	239	+ n_ids = vertices_staircase.size()
	240	+
	241	+ cdef vertex_t* vi
	242	+ cdef vertex_t* vj
	243	+ cdef size_t msize
	244	+
	245	+ msize = mesh.vertices.shape[0]
	246	+ cdef vector[weight_t]* weights = new vector[weight_t](msize)
	247	+ weights.assign(msize, bmin)
	248	+
	249	+ for i in prange(n_ids, nogil=True):
	250	+ vi_id = vertices_staircase[i]
	251	+ deref(weights)[vi_id] = 1.0
	252	+
	253	+ vi = &mesh.vertices[vi_id, 0]
	254	+ near_vertices = mesh.get_near_vertices_to_v(vi_id, tmax)
	255	+ nnv = near_vertices.size()
	256	+
	257	+ for j in xrange(nnv):
	258	+ vj_id = deref(near_vertices)[j]
	259	+ vj = &mesh.vertices[vj_id, 0]
	260	+
	261	+ d = sqrt((vi[0] - vj[0]) * (vi[0] - vj[0])\
	262	+ + (vi[1] - vj[1]) * (vi[1] - vj[1])\
	263	+ + (vi[2] - vj[2]) * (vi[2] - vj[2]))
	264	+ value = (1.0 - d/tmax) * (1.0 - bmin) + bmin
	265	+
	266	+ if value > deref(weights)[vj_id]:
	267	+ deref(weights)[vj_id] = value
	268	+
	269	+ del near_vertices
	270	+
	271	+ # for i in xrange(msize):
	272	+ # if mesh.is_border(i):
	273	+ # deref(weights)[i] = 0.0
	274	+
	275	+ # cdef vertex_id_t v0, v1, v2
	276	+ # for i in xrange(mesh.faces.shape[0]):
	277	+ # for j in xrange(1, 4):
	278	+ # v0 = mesh.faces[i, j]
	279	+ # vi = &mesh.vertices[v0, 0]
	280	+ # if mesh.is_border(v0):
	281	+ # deref(weights)[v0] = 0.0
	282	+ # v1 = mesh.faces[i, (j + 1) % 3 + 1]
	283	+ # if mesh.is_border(v1):
	284	+ # vi = &mesh.vertices[v1, 0]
	285	+ # deref(weights)[v0] = 0.0
	286	+
	287	+ return weights
	288	+
	289	+
	290	+cdef inline Point calc_d(Mesh mesh, vertex_id_t v_id) nogil:
	291	+ cdef Point D
	292	+ cdef int nf, f_id, nid
	293	+ cdef float n=0
	294	+ cdef int i
	295	+ cdef vertex_t* vi
	296	+ cdef vertex_t* vj
	297	+ cdef set[vertex_id_t]* vertices
	298	+ cdef set[vertex_id_t].iterator it
	299	+ cdef vertex_id_t vj_id
	300	+
	301	+ D.x = 0.0
	302	+ D.y = 0.0
	303	+ D.z = 0.0
	304	+
	305	+ vertices = mesh.get_ring1(v_id)
	306	+ vi = &mesh.vertices[v_id, 0]
	307	+
	308	+ if mesh.is_border(v_id):
	309	+ it = vertices.begin()
	310	+ while it != vertices.end():
	311	+ vj_id = deref(it)
	312	+ if mesh.is_border(vj_id):
	313	+ vj = &mesh.vertices[vj_id, 0]
	314	+
	315	+ D.x = D.x + (vi[0] - vj[0])
	316	+ D.y = D.y + (vi[1] - vj[1])
	317	+ D.z = D.z + (vi[2] - vj[2])
	318	+ n += 1.0
	319	+
	320	+ inc(it)
	321	+ else:
	322	+ it = vertices.begin()
	323	+ while it != vertices.end():
	324	+ vj_id = deref(it)
	325	+ vj = &mesh.vertices[vj_id, 0]
	326	+
	327	+ D.x = D.x + (vi[0] - vj[0])
	328	+ D.y = D.y + (vi[1] - vj[1])
	329	+ D.z = D.z + (vi[2] - vj[2])
	330	+ n += 1.0
	331	+
	332	+ inc(it)
	333	+
	334	+ del vertices
	335	+
	336	+ D.x = D.x / n
	337	+ D.y = D.y / n
	338	+ D.z = D.z / n
	339	+ return D
	340	+
	341	+cdef vector[vertex_id_t]* find_staircase_artifacts(Mesh mesh, double[3] stack_orientation, double T) nogil:
	342	+ """
	343	+ This function is used to find vertices at staircase artifacts, which are
	344	+ those vertices whose incident faces' orientation differences are
	345	+ greater than T.
	346	+
	347	+ Params:
	348	+ mesh: Mesh
	349	+ stack_orientation: orientation of slice stacking
	350	+ T: Min angle (between vertex faces and stack_orientation) to consider a
	351	+ vertex a staircase artifact.
	352	+ """
	353	+ cdef int nv, nf, f_id, v_id
	354	+ cdef double of_z, of_y, of_x, min_z, max_z, min_y, max_y, min_x, max_x;
	355	+ cdef vector[vertex_id_t]* f_ids
	356	+ cdef normal_t* normal
	357	+
	358	+ cdef vector[vertex_id_t]* output = new vector[vertex_id_t]()
	359	+ cdef int i
	360	+
	361	+ nv = mesh.vertices.shape[0]
	362	+
	363	+ for v_id in xrange(nv):
	364	+ max_z = -10000
	365	+ min_z = 10000
	366	+ max_y = -10000
	367	+ min_y = 10000
	368	+ max_x = -10000
	369	+ min_x = 10000
	370	+
	371	+ f_ids = mesh.get_faces_by_vertex(v_id)
	372	+ nf = deref(f_ids).size()
	373	+
	374	+ for i in xrange(nf):
	375	+ f_id = deref(f_ids)[i]
	376	+ normal = &mesh.normals[f_id, 0]
	377	+
	378	+ of_z = 1 - fabs(normal[0]stack_orientation[0] + normal[1]stack_orientation[1] + normal[2]*stack_orientation[2]);
	379	+ of_y = 1 - fabs(normal[0]0 + normal[1]1 + normal[2]*0);
	380	+ of_x = 1 - fabs(normal[0]1 + normal[1]0 + normal[2]*0);
	381	+
	382	+ if (of_z > max_z):
	383	+ max_z = of_z
	384	+
	385	+ if (of_z < min_z):
	386	+ min_z = of_z
	387	+
	388	+ if (of_y > max_y):
	389	+ max_y = of_y
	390	+
	391	+ if (of_y < min_y):
	392	+ min_y = of_y
	393	+
	394	+ if (of_x > max_x):
	395	+ max_x = of_x
	396	+
	397	+ if (of_x < min_x):
	398	+ min_x = of_x
	399	+
	400	+
	401	+ if ((fabs(max_z - min_z) >= T) or (fabs(max_y - min_y) >= T) or (fabs(max_x - min_x) >= T)):
	402	+ output.push_back(v_id)
	403	+ break
	404	+ return output
	405	+
	406	+
	407	+cdef void taubin_smooth(Mesh mesh, vector[weight_t]& weights, float l, float m, int steps):
	408	+ """
	409	+ Implementation of Taubin's smooth algorithm described in the paper "A
	410	+ Signal Processing Approach To Fair Surface Design". His benefeat is it
	411	+ avoids surface shrinking.
	412	+ """
	413	+ cdef int s, i, nvertices
	414	+ nvertices = mesh.vertices.shape[0]
	415	+ cdef vector[Point] D = vector[Point](nvertices)
	416	+ cdef vertex_t* vi
	417	+ for s in xrange(steps):
	418	+ for i in prange(nvertices, nogil=True):
	419	+ D[i] = calc_d(mesh, i)
	420	+
	421	+ for i in prange(nvertices, nogil=True):
	422	+ mesh.vertices[i, 0] += weights[i]lD[i].x;
	423	+ mesh.vertices[i, 1] += weights[i]lD[i].y;
	424	+ mesh.vertices[i, 2] += weights[i]lD[i].z;
	425	+
	426	+ for i in prange(nvertices, nogil=True):
	427	+ D[i] = calc_d(mesh, i)
	428	+
	429	+ for i in prange(nvertices, nogil=True):
	430	+ mesh.vertices[i, 0] += weights[i]mD[i].x;
	431	+ mesh.vertices[i, 1] += weights[i]mD[i].y;
	432	+ mesh.vertices[i, 2] += weights[i]mD[i].z;
	433	+
	434	+
	435	+def ca_smoothing(Mesh mesh, double T, double tmax, double bmin, int n_iters):
	436	+ """
	437	+ This is a implementation of the paper "Context-aware mesh smoothing for
	438	+ biomedical applications". It can be used to smooth meshes generated by
	439	+ binary images to remove its staircase artifacts and keep the fine features.
	440	+
	441	+ Params:
	442	+ mesh: Mesh
	443	+ T: Min angle (between vertex faces and stack_orientation) to consider a
	444	+ vertex a staircase artifact
	445	+ tmax: max distance the vertex must be to its nearest artifact vertex
	446	+ to considered to calculate the weight
	447	+ bmin: The minimun weight
	448	+ n_iters: Number of iterations.
	449	+ """
	450	+ cdef double[3] stack_orientation = [0.0, 0.0, 1.0]
	451	+
	452	+ t0 = time.time()
	453	+ cdef vector[vertex_id_t]* vertices_staircase = find_staircase_artifacts(mesh, stack_orientation, T)
	454	+ print "vertices staircase", time.time() - t0
	455	+
	456	+ t0 = time.time()
	457	+ cdef vector[weight_t]* weights = calc_artifacts_weight(mesh, deref(vertices_staircase), tmax, bmin)
	458	+ print "Weights", time.time() - t0
	459	+
	460	+ del vertices_staircase
	461	+
	462	+ t0 = time.time()
	463	+ taubin_smooth(mesh, deref(weights), 0.5, -0.53, n_iters)
	464	+ print "taubin", time.time() - t0
	465	+
	466	+ del weights
	@@ -10,3 +10,7 @@ ctypedef fused image_t:		@@ -10,3 +10,7 @@ ctypedef fused image_t:
10	np.uint8_t	10	np.uint8_t
11		11
12	ctypedef np.uint8_t mask_t	12	ctypedef np.uint8_t mask_t
		13	+
		14	+ctypedef np.float32_t vertex_t
		15	+ctypedef np.float32_t normal_t
		16	+ctypedef np.int64_t vertex_id_t
	@@ -42,6 +42,7 @@ try:		@@ -42,6 +42,7 @@ try:
42	except ImportError:	42	except ImportError:
43	import data.ca_smoothing as ca_smoothing	43	import data.ca_smoothing as ca_smoothing
44		44
		45	+import cy_mesh
45	# TODO: Verificar ReleaseDataFlagOn and SetSource	46	# TODO: Verificar ReleaseDataFlagOn and SetSource
46		47
47	class Surface():	48	class Surface():
	@@ -564,16 +565,28 @@ class SurfaceManager():		@@ -564,16 +565,28 @@ class SurfaceManager():
564	# polydata.SetSource(None)	565	# polydata.SetSource(None)
565	del clean	566	del clean
566		567
567	- try:
568	- polydata.BuildLinks()
569	- except TypeError:
570	- polydata.BuildLinks(0)
571	- polydata = ca_smoothing.ca_smoothing(polydata, options['angle'],
572	- options['max distance'],
573	- options['min weight'],
574	- options['steps'])	568	+ # try:
		569	+ # polydata.BuildLinks()
		570	+ # except TypeError:
		571	+ # polydata.BuildLinks(0)
		572	+ # polydata = ca_smoothing.ca_smoothing(polydata, options['angle'],
		573	+ # options['max distance'],
		574	+ # options['min weight'],
		575	+ # options['steps'])
		576	+
		577	+ mesh = cy_mesh.Mesh(polydata)
		578	+ cy_mesh.ca_smoothing(mesh, options['angle'],
		579	+ options['max distance'],
		580	+ options['min weight'],
		581	+ options['steps'])
		582	+ # polydata = mesh.to_vtk()
		583	+
575	# polydata.SetSource(None)	584	# polydata.SetSource(None)
576	# polydata.DebugOn()	585	# polydata.DebugOn()
		586	+ w = vtk.vtkPLYWriter()
		587	+ w.SetInputData(polydata)
		588	+ w.SetFileName('/tmp/ca_smoothing_inv.ply')
		589	+ w.Write()
577		590
578	else:	591	else:
579	#smoother = vtk.vtkWindowedSincPolyDataFilter()	592	#smoother = vtk.vtkWindowedSincPolyDataFilter()
	@@ -1276,7 +1276,7 @@ class CAOptions(wx.Panel):		@@ -1276,7 +1276,7 @@ class CAOptions(wx.Panel):
1276	max_val=100.0, increment=0.1,	1276	max_val=100.0, increment=0.1,
1277	digits=2)	1277	digits=2)
1278		1278
1279	- self.min_weight = floatspin.FloatSpin(self, -1, value=0.2, min_val=0.0,	1279	+ self.min_weight = floatspin.FloatSpin(self, -1, value=0.5, min_val=0.0,
1280	max_val=1.0, increment=0.1,	1280	max_val=1.0, increment=0.1,
1281	digits=1)	1281	digits=1)
1282		1282
	@@ -29,6 +29,12 @@ if sys.platform == 'linux2':		@@ -29,6 +29,12 @@ if sys.platform == 'linux2':
29	Extension("invesalius.data.floodfill", ["invesalius/data/floodfill.pyx"],	29	Extension("invesalius.data.floodfill", ["invesalius/data/floodfill.pyx"],
30	include_dirs=[numpy.get_include()],	30	include_dirs=[numpy.get_include()],
31	language='c++',),	31	language='c++',),
		32	+
		33	+ Extension("invesalius.data.cy_mesh", ["invesalius/data/cy_mesh.pyx"],
		34	+ include_dirs=[numpy.get_include()],
		35	+ extra_compile_args=['-fopenmp', '-std=c++11'],
		36	+ extra_link_args=['-fopenmp', '-std=c++11'],
		37	+ language='c++',),
32	])	38	])
33	)	39	)
34		40
	@@ -50,6 +56,11 @@ elif sys.platform == 'win32':		@@ -50,6 +56,11 @@ elif sys.platform == 'win32':
50	Extension("invesalius.data.floodfill", ["invesalius/data/floodfill.pyx"],	56	Extension("invesalius.data.floodfill", ["invesalius/data/floodfill.pyx"],
51	include_dirs=[numpy.get_include()],	57	include_dirs=[numpy.get_include()],
52	language='c++',),	58	language='c++',),
		59	+
		60	+ Extension("invesalius.data.cy_mesh", ["invesalius/data/cy_mesh.pyx"],
		61	+ include_dirs=[numpy.get_include()],
		62	+ extra_compile_args=['/openmp',],
		63	+ language='c++',),
53	])	64	])
54	)	65	)
55		66
	@@ -75,5 +86,12 @@ else:		@@ -75,5 +86,12 @@ else:
75	Extension("invesalius.data.floodfill", ["invesalius/data/floodfill.pyx"],	86	Extension("invesalius.data.floodfill", ["invesalius/data/floodfill.pyx"],
76	include_dirs=[numpy.get_include()],	87	include_dirs=[numpy.get_include()],
77	language='c++',),	88	language='c++',),
		89	+
		90	+ Extension("invesalius.data.cy_mesh", ["invesalius/data/cy_mesh.pyx"],
		91	+ include_dirs=[numpy.get_include()],
		92	+ extra_compile_args=['-fopenmp', '-std=c++11'],
		93	+ extra_link_args=['-fopenmp', '-std=c++11'],
		94	+ language='c++',),
		95	+
78	])	96	])
79	)	97	)