From f13a4b79d6e6cc3f3ac885ccd20527fde9106719 Mon Sep 17 00:00:00 2001 From: Thiago Franco de Moraes Date: Wed, 9 Nov 2016 11:03:45 -0200 Subject: [PATCH] Casmoothing implemented in cython (#61) --- .gitignore | 6 ++++++ invesalius/data/cy_mesh.pyx | 466 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ invesalius/data/cy_my_types.pxd | 4 ++++ invesalius/data/surface.py | 29 +++++++++++++++++++++-------- invesalius/gui/dialogs.py | 2 +- setup.py | 18 ++++++++++++++++++ 6 files changed, 516 insertions(+), 9 deletions(-) create mode 100644 invesalius/data/cy_mesh.pyx diff --git a/.gitignore b/.gitignore index d1b6ccc..4ad10a2 100644 --- a/.gitignore +++ b/.gitignore @@ -24,3 +24,9 @@ tags build *.patch *.tgz + +*.pyd +*.cpp +*.diff + +*.directory diff --git a/invesalius/data/cy_mesh.pyx b/invesalius/data/cy_mesh.pyx new file mode 100644 index 0000000..c89da9d --- /dev/null +++ b/invesalius/data/cy_mesh.pyx @@ -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 = 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 diff --git a/invesalius/data/cy_my_types.pxd b/invesalius/data/cy_my_types.pxd index 36bd784..f87d0b0 100644 --- a/invesalius/data/cy_my_types.pxd +++ b/invesalius/data/cy_my_types.pxd @@ -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 diff --git a/invesalius/data/surface.py b/invesalius/data/surface.py index 0ec727e..d97e028 100644 --- a/invesalius/data/surface.py +++ b/invesalius/data/surface.py @@ -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() diff --git a/invesalius/gui/dialogs.py b/invesalius/gui/dialogs.py index d000e05..451a4c1 100644 --- a/invesalius/gui/dialogs.py +++ b/invesalius/gui/dialogs.py @@ -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) diff --git a/setup.py b/setup.py index c9223db..cc05315 100644 --- a/setup.py +++ b/setup.py @@ -29,6 +29,12 @@ if sys.platform == 'linux2': 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 == 'win32': 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++',), + ]) ) -- libgit2 0.21.2