From 18c2bcddd4fbb36ebbef13006c804f4ce661476a Mon Sep 17 00:00:00 2001 From: Thiago Franco de Moraes Date: Fri, 20 Mar 2020 14:35:27 -0300 Subject: [PATCH] Moved all cython codes to invesalius_cy --- invesalius/data/cy_mesh.pyx | 466 ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- invesalius/data/cy_my_types.pxd | 21 --------------------- invesalius/data/floodfill.pyx | 274 ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- invesalius/data/interpolation.pxd | 8 -------- invesalius/data/interpolation.pyx | 392 -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- invesalius/data/mask.py | 2 +- invesalius/data/mips.pyx | 379 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- invesalius/data/slice_.py | 4 ++-- invesalius/data/styles.py | 2 +- invesalius/data/surface.py | 4 ++-- invesalius/data/surface_process.py | 2 +- invesalius/data/transforms.pyx | 174 ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ invesalius_cy/cy_mesh.pyx | 466 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ invesalius_cy/cy_my_types.pxd | 21 +++++++++++++++++++++ invesalius_cy/floodfill.pyx | 274 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ invesalius_cy/interpolation.pxd | 8 ++++++++ invesalius_cy/interpolation.pyx | 392 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ invesalius_cy/mips.pyx | 379 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ invesalius_cy/transforms.pyx | 174 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ setup.py | 20 ++++++++++---------- 20 files changed, 1731 insertions(+), 1731 deletions(-) delete mode 100644 invesalius/data/cy_mesh.pyx delete mode 100644 invesalius/data/cy_my_types.pxd delete mode 100644 invesalius/data/floodfill.pyx delete mode 100644 invesalius/data/interpolation.pxd delete mode 100644 invesalius/data/interpolation.pyx delete mode 100644 invesalius/data/mips.pyx delete mode 100644 invesalius/data/transforms.pyx create mode 100644 invesalius_cy/cy_mesh.pyx create mode 100644 invesalius_cy/cy_my_types.pxd create mode 100644 invesalius_cy/floodfill.pyx create mode 100644 invesalius_cy/interpolation.pxd create mode 100644 invesalius_cy/interpolation.pyx create mode 100644 invesalius_cy/mips.pyx create mode 100644 invesalius_cy/transforms.pyx diff --git a/invesalius/data/cy_mesh.pyx b/invesalius/data/cy_mesh.pyx deleted file mode 100644 index 6d074f8..0000000 --- a/invesalius/data/cy_mesh.pyx +++ /dev/null @@ -1,466 +0,0 @@ -# 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 maximum 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 minimum 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 minimum 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 deleted file mode 100644 index 13d6c6c..0000000 --- a/invesalius/data/cy_my_types.pxd +++ /dev/null @@ -1,21 +0,0 @@ -import numpy as np -cimport numpy as np -cimport cython - -# ctypedef np.uint16_t image_t - -ctypedef fused image_t: - np.float64_t - np.int16_t - np.uint8_t - -ctypedef np.uint8_t mask_t - -ctypedef np.float32_t vertex_t -ctypedef np.float32_t normal_t - -# To compile in windows 64 -IF UNAME_MACHINE == 'AMD64': - ctypedef np.int64_t vertex_id_t -ELSE: - ctypedef np.int_t vertex_id_t diff --git a/invesalius/data/floodfill.pyx b/invesalius/data/floodfill.pyx deleted file mode 100644 index 6636cc9..0000000 --- a/invesalius/data/floodfill.pyx +++ /dev/null @@ -1,274 +0,0 @@ -import numpy as np -cimport numpy as np -cimport cython - -from collections import deque - -from cython.parallel import prange -from libc.math cimport floor, ceil -from libcpp cimport bool -from libcpp.deque cimport deque as cdeque -from libcpp.vector cimport vector - -from cy_my_types cimport image_t, mask_t - -cdef struct s_coord: - int x - int y - int z - -ctypedef s_coord coord - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.wraparound(False) -@cython.nonecheck(False) -@cython.cdivision(True) -cdef inline void append_queue(cdeque[int]& stack, int x, int y, int z, int d, int h, int w) nogil: - stack.push_back(z*h*w + y*w + x) - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.wraparound(False) -@cython.nonecheck(False) -@cython.cdivision(True) -cdef inline void pop_queue(cdeque[int]& stack, int* x, int* y, int* z, int d, int h, int w) nogil: - cdef int i = stack.front() - stack.pop_front() - x[0] = i % w - y[0] = (i / w) % h - z[0] = i / (h * w) - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -def floodfill(np.ndarray[image_t, ndim=3] data, int i, int j, int k, int v, int fill, np.ndarray[mask_t, ndim=3] out): - - cdef int to_return = 0 - if out is None: - out = np.zeros_like(data) - to_return = 1 - - cdef int x, y, z - cdef int w, h, d - - d = data.shape[0] - h = data.shape[1] - w = data.shape[2] - - stack = [(i, j, k), ] - out[k, j, i] = fill - - while stack: - x, y, z = stack.pop() - - if z + 1 < d and data[z + 1, y, x] == v and out[z + 1, y, x] != fill: - out[z + 1, y, x] = fill - stack.append((x, y, z + 1)) - - if z - 1 >= 0 and data[z - 1, y, x] == v and out[z - 1, y, x] != fill: - out[z - 1, y, x] = fill - stack.append((x, y, z - 1)) - - if y + 1 < h and data[z, y + 1, x] == v and out[z, y + 1, x] != fill: - out[z, y + 1, x] = fill - stack.append((x, y + 1, z)) - - if y - 1 >= 0 and data[z, y - 1, x] == v and out[z, y - 1, x] != fill: - out[z, y - 1, x] = fill - stack.append((x, y - 1, z)) - - if x + 1 < w and data[z, y, x + 1] == v and out[z, y, x + 1] != fill: - out[z, y, x + 1] = fill - stack.append((x + 1, y, z)) - - if x - 1 >= 0 and data[z, y, x - 1] == v and out[z, y, x - 1] != fill: - out[z, y, x - 1] = fill - stack.append((x - 1, y, z)) - - if to_return: - return out - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.wraparound(False) -@cython.nonecheck(False) -def floodfill_threshold(np.ndarray[image_t, ndim=3] data, list seeds, int t0, int t1, int fill, np.ndarray[mask_t, ndim=3] strct, np.ndarray[mask_t, ndim=3] out): - - cdef int to_return = 0 - if out is None: - out = np.zeros_like(data) - to_return = 1 - - cdef int x, y, z - cdef int dx, dy, dz - cdef int odx, ody, odz - cdef int xo, yo, zo - cdef int i, j, k - cdef int offset_x, offset_y, offset_z - - dz = data.shape[0] - dy = data.shape[1] - dx = data.shape[2] - - odz = strct.shape[0] - ody = strct.shape[1] - odx = strct.shape[2] - - cdef cdeque[coord] stack - cdef coord c - - offset_z = odz / 2 - offset_y = ody / 2 - offset_x = odx / 2 - - for i, j, k in seeds: - if data[k, j, i] >= t0 and data[k, j, i] <= t1: - c.x = i - c.y = j - c.z = k - stack.push_back(c) - out[k, j, i] = fill - - with nogil: - while stack.size(): - c = stack.back() - stack.pop_back() - - x = c.x - y = c.y - z = c.z - - out[z, y, x] = fill - - for k in xrange(odz): - zo = z + k - offset_z - for j in xrange(ody): - yo = y + j - offset_y - for i in xrange(odx): - if strct[k, j, i]: - xo = x + i - offset_x - if 0 <= xo < dx and 0 <= yo < dy and 0 <= zo < dz and out[zo, yo, xo] != fill and t0 <= data[zo, yo, xo] <= t1: - out[zo, yo, xo] = fill - c.x = xo - c.y = yo - c.z = zo - stack.push_back(c) - - if to_return: - return out - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.wraparound(False) -@cython.nonecheck(False) -def floodfill_auto_threshold(np.ndarray[image_t, ndim=3] data, list seeds, float p, int fill, np.ndarray[mask_t, ndim=3] out): - - cdef int to_return = 0 - if out is None: - out = np.zeros_like(data) - to_return = 1 - - cdef cdeque[int] stack - cdef int x, y, z - cdef int w, h, d - cdef int xo, yo, zo - cdef int t0, t1 - - cdef int i, j, k - - d = data.shape[0] - h = data.shape[1] - w = data.shape[2] - - - # stack = deque() - - x = 0 - y = 0 - z = 0 - - - for i, j, k in seeds: - append_queue(stack, i, j, k, d, h, w) - out[k, j, i] = fill - print i, j, k, d, h, w - - with nogil: - while stack.size(): - pop_queue(stack, &x, &y, &z, d, h, w) - - # print x, y, z, d, h, w - - xo = x - yo = y - zo = z - - t0 = ceil(data[z, y, x] * (1 - p)) - t1 = floor(data[z, y, x] * (1 + p)) - - if z + 1 < d and data[z + 1, y, x] >= t0 and data[z + 1, y, x] <= t1 and out[zo + 1, yo, xo] != fill: - out[zo + 1, yo, xo] = fill - append_queue(stack, x, y, z+1, d, h, w) - - if z - 1 >= 0 and data[z - 1, y, x] >= t0 and data[z - 1, y, x] <= t1 and out[zo - 1, yo, xo] != fill: - out[zo - 1, yo, xo] = fill - append_queue(stack, x, y, z-1, d, h, w) - - if y + 1 < h and data[z, y + 1, x] >= t0 and data[z, y + 1, x] <= t1 and out[zo, yo + 1, xo] != fill: - out[zo, yo + 1, xo] = fill - append_queue(stack, x, y+1, z, d, h, w) - - if y - 1 >= 0 and data[z, y - 1, x] >= t0 and data[z, y - 1, x] <= t1 and out[zo, yo - 1, xo] != fill: - out[zo, yo - 1, xo] = fill - append_queue(stack, x, y-1, z, d, h, w) - - if x + 1 < w and data[z, y, x + 1] >= t0 and data[z, y, x + 1] <= t1 and out[zo, yo, xo + 1] != fill: - out[zo, yo, xo + 1] = fill - append_queue(stack, x+1, y, z, d, h, w) - - if x - 1 >= 0 and data[z, y, x - 1] >= t0 and data[z, y, x - 1] <= t1 and out[zo, yo, xo - 1] != fill: - out[zo, yo, xo - 1] = fill - append_queue(stack, x-1, y, z, d, h, w) - - if to_return: - return out - - -@cython.boundscheck(False) -@cython.wraparound(False) -@cython.nonecheck(False) -def fill_holes_automatically(np.ndarray[mask_t, ndim=3] mask, np.ndarray[np.uint16_t, ndim=3] labels, unsigned int nlabels, unsigned int max_size): - """ - Fill mask holes automatically. The hole must <= max_size. Return True if any hole were filled. - """ - cdef np.ndarray[np.uint32_t, ndim=1] sizes = np.zeros(shape=(nlabels + 1), dtype=np.uint32) - cdef int x, y, z - cdef int dx, dy, dz - cdef int i - - cdef bool modified = False - - dz = mask.shape[0] - dy = mask.shape[1] - dx = mask.shape[2] - - for z in xrange(dz): - for y in xrange(dy): - for x in xrange(dx): - sizes[labels[z, y, x]] += 1 - - #Checking if any hole will be filled - for i in xrange(nlabels + 1): - if sizes[i] <= max_size: - modified = True - - if not modified: - return 0 - - for z in prange(dz, nogil=True): - for y in xrange(dy): - for x in xrange(dx): - if sizes[labels[z, y, x]] <= max_size: - mask[z, y, x] = 254 - - return modified diff --git a/invesalius/data/interpolation.pxd b/invesalius/data/interpolation.pxd deleted file mode 100644 index b63cc41..0000000 --- a/invesalius/data/interpolation.pxd +++ /dev/null @@ -1,8 +0,0 @@ -from .cy_my_types cimport image_t - -cdef double interpolate(image_t[:, :, :], double, double, double) nogil -cdef double tricub_interpolate(image_t[:, :, :], double, double, double) nogil -cdef double tricubicInterpolate (image_t[:, :, :], double, double, double) nogil -cdef double lanczos3 (image_t[:, :, :], double, double, double) nogil - -cdef double nearest_neighbour_interp(image_t[:, :, :], double, double, double) nogil diff --git a/invesalius/data/interpolation.pyx b/invesalius/data/interpolation.pyx deleted file mode 100644 index 9f44648..0000000 --- a/invesalius/data/interpolation.pyx +++ /dev/null @@ -1,392 +0,0 @@ -# from interpolation cimport interpolate - -import numpy as np -cimport numpy as np -cimport cython - -from libc.math cimport floor, ceil, sqrt, fabs, sin, M_PI -from cython.parallel import prange - -DEF LANCZOS_A = 4 -DEF SIZE_LANCZOS_TMP = LANCZOS_A * 2 - 1 - -cdef double[64][64] temp = [ - [ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 2, -2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 9, -9,-9, 9, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [-6, 6, 6,-6, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [-6, 6, 6,-6, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 4, -4,-4, 4, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0], - [-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 9, -9, 0, 0,-9, 9, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [-6, 6, 0, 0, 6,-6, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9, 0, 0,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0], - [ 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0], - [-27, 27,27,-27,27,-27,-27,27,-18,-9,18, 9,18, 9,-18,-9,-18,18,-9, 9,18,-18, 9,-9,-18,18,18,-18,-9, 9, 9,-9,-12,-6,-6,-3,12, 6, 6, 3,-12,-6,12, 6,-6,-3, 6, 3,-12,12,-6, 6,-6, 6,-3, 3,-8,-4,-4,-2,-4,-2,-2,-1], - [18, -18,-18,18,-18,18,18,-18, 9, 9,-9,-9,-9,-9, 9, 9,12,-12, 6,-6,-12,12,-6, 6,12,-12,-12,12, 6,-6,-6, 6, 6, 6, 3, 3,-6,-6,-3,-3, 6, 6,-6,-6, 3, 3,-3,-3, 8,-8, 4,-4, 4,-4, 2,-2, 4, 4, 2, 2, 2, 2, 1, 1], - [-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0], - [18, -18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6, 9,-9, 9,-9,-9, 9,-9, 9,12,-12,-12,12, 6,-6,-6, 6, 6, 3, 6, 3,-6,-3,-6,-3, 8, 4,-8,-4, 4, 2,-4,-2, 6,-6, 6,-6, 3,-3, 3,-3, 4, 2, 4, 2, 2, 1, 2, 1], - [-12, 12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-6, 6,-6, 6, 6,-6, 6,-6,-8, 8, 8,-8,-4, 4, 4,-4,-3,-3,-3,-3, 3, 3, 3, 3,-4,-4, 4, 4,-2,-2, 2, 2,-4, 4,-4, 4,-2, 2,-2, 2,-2,-2,-2,-2,-1,-1,-1,-1], - [ 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [-6, 6, 0, 0, 6,-6, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 4, -4, 0, 0,-4, 4, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4, 0, 0,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0], - [-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0], - [18, -18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6,12,-12, 6,-6,-12,12,-6, 6, 9,-9,-9, 9, 9,-9,-9, 9, 8, 4, 4, 2,-8,-4,-4,-2, 6, 3,-6,-3, 6, 3,-6,-3, 6,-6, 3,-3, 6,-6, 3,-3, 4, 2, 2, 1, 4, 2, 2, 1], - [-12, 12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-8, 8,-4, 4, 8,-8, 4,-4,-6, 6, 6,-6,-6, 6, 6,-6,-4,-4,-2,-2, 4, 4, 2, 2,-3,-3, 3, 3,-3,-3, 3, 3,-4, 4,-2, 2,-4, 4,-2, 2,-2,-2,-1,-1,-2,-2,-1,-1], - [ 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0], - [-12, 12,12,-12,12,-12,-12,12,-8,-4, 8, 4, 8, 4,-8,-4,-6, 6,-6, 6, 6,-6, 6,-6,-6, 6, 6,-6,-6, 6, 6,-6,-4,-2,-4,-2, 4, 2, 4, 2,-4,-2, 4, 2,-4,-2, 4, 2,-3, 3,-3, 3,-3, 3,-3, 3,-2,-1,-2,-1,-2,-1,-2,-1], - [ 8, -8,-8, 8,-8, 8, 8,-8, 4, 4,-4,-4,-4,-4, 4, 4, 4,-4, 4,-4,-4, 4,-4, 4, 4,-4,-4, 4, 4,-4,-4, 4, 2, 2, 2, 2,-2,-2,-2,-2, 2, 2,-2,-2, 2, 2,-2,-2, 2,-2, 2,-2, 2,-2, 2,-2, 1, 1, 1, 1, 1, 1, 1, 1] -] - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -cdef double nearest_neighbour_interp(image_t[:, :, :] V, double x, double y, double z) nogil: - return V[(z), (y), (x)] - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -cdef double interpolate(image_t[:, :, :] V, double x, double y, double z) nogil: - cdef double xd, yd, zd - cdef double c00, c10, c01, c11 - cdef double c0, c1 - cdef double c - - cdef int x0 = floor(x) - cdef int x1 = x0 + 1 - - cdef int y0 = floor(y) - cdef int y1 = y0 + 1 - - cdef int z0 = floor(z) - cdef int z1 = z0 + 1 - - if x0 == x1: - xd = 1.0 - else: - xd = (x - x0) / (x1 - x0) - - if y0 == y1: - yd = 1.0 - else: - yd = (y - y0) / (y1 - y0) - - if z0 == z1: - zd = 1.0 - else: - zd = (z - z0) / (z1 - z0) - - c00 = _G(V, x0, y0, z0)*(1 - xd) + _G(V, x1, y0, z0)*xd - c10 = _G(V, x0, y1, z0)*(1 - xd) + _G(V, x1, y1, z0)*xd - c01 = _G(V, x0, y0, z1)*(1 - xd) + _G(V, x1, y0, z1)*xd - c11 = _G(V, x0, y1, z1)*(1 - xd) + _G(V, x1, y1, z1)*xd - - c0 = c00*(1 - yd) + c10*yd - c1 = c01*(1 - yd) + c11*yd - - c = c0*(1 - zd) + c1*zd - - return c - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -cdef inline double lanczos3_L(double x, int a) nogil: - if x == 0: - return 1.0 - elif -a <= x < a: - return (a * sin(M_PI * x) * sin(M_PI * (x / a)))/(M_PI**2 * x**2) - else: - return 0.0 - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -cdef double lanczos3(image_t[:, :, :] V, double x, double y, double z) nogil: - cdef int a = LANCZOS_A - - cdef int xd = floor(x) - cdef int yd = floor(y) - cdef int zd = floor(z) - - cdef int xi = xd - a + 1 - cdef int xf = xd + a - - cdef int yi = yd - a + 1 - cdef int yf = yd + a - - cdef int zi = zd - a + 1 - cdef int zf = zd + a - - cdef double lx = 0.0 - cdef double ly = 0.0 - cdef double lz = 0.0 - - cdef double[SIZE_LANCZOS_TMP][SIZE_LANCZOS_TMP] temp_x - cdef double[SIZE_LANCZOS_TMP] temp_y - - cdef int i, j, k - cdef int m, n, o - - m = 0 - for k in xrange(zi, zf): - n = 0 - for j in xrange(yi, yf): - lx = 0 - for i in xrange(xi, xf): - lx += _G(V, i, j, k) * lanczos3_L(x - i, a) - temp_x[m][n] = lx - n += 1 - m += 1 - - m = 0 - for k in xrange(zi, zf): - n = 0 - ly = 0 - for j in xrange(yi, yf): - ly += temp_x[m][n] * lanczos3_L(y - j, a) - n += 1 - temp_y[m] = ly - m += 1 - - m = 0 - for k in xrange(zi, zf): - lz += temp_y[m] * lanczos3_L(z - k, a) - m += 1 - - return lz - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -cdef image_t _G(image_t[:, :, :] V, int x, int y, int z) nogil: - cdef int dz, dy, dx - dz = V.shape[0] - 1 - dy = V.shape[1] - 1 - dx = V.shape[2] - 1 - - if x < 0: - x = dx + x + 1 - elif x > dx: - x = x - dx - 1 - - if y < 0: - y = dy + y + 1 - elif y > dy: - y = y - dy - 1 - - if z < 0: - z = dz + z + 1 - elif z > dz: - z = z - dz - 1 - - return V[z, y, x] - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -cdef void calc_coef_tricub(image_t[:, :, :] V, double x, double y, double z, double [64] coef) nogil: - cdef int xi = floor(x) - cdef int yi = floor(y) - cdef int zi = floor(z) - - cdef double[64] _x - - cdef int i, j - - _x[0] = _G(V, xi, yi, zi) - _x[1] = _G(V, xi + 1, yi, zi) - _x[2] = _G(V, xi, yi + 1, zi) - _x[3] = _G(V, xi + 1, yi + 1, zi) - _x[4] = _G(V, xi, yi, zi + 1) - _x[5] = _G(V, xi + 1, yi, zi + 1) - _x[6] = _G(V, xi, yi + 1, zi + 1) - _x[7] = _G(V, xi + 1, yi + 1, zi + 1) - - _x[8] = 0.5*(_G(V, xi+1,yi,zi) - _G(V, xi-1, yi, zi)) - _x[9] = 0.5*(_G(V, xi+2,yi,zi) - _G(V, xi, yi, zi)) - _x[10] = 0.5*(_G(V, xi+1, yi+1,zi) - _G(V, xi-1, yi+1, zi)) - _x[11] = 0.5*(_G(V, xi+2, yi+1,zi) - _G(V, xi, yi+1, zi)) - _x[12] = 0.5*(_G(V, xi+1, yi,zi+1) - _G(V, xi-1, yi, zi+1)) - _x[13] = 0.5*(_G(V, xi+2, yi,zi+1) - _G(V, xi, yi, zi+1)) - _x[14] = 0.5*(_G(V, xi+1, yi+1,zi+1) - _G(V, xi-1, yi+1, zi+1)) - _x[15] = 0.5*(_G(V, xi+2, yi+1,zi+1) - _G(V, xi, yi+1, zi+1)) - _x[16] = 0.5*(_G(V, xi, yi+1,zi) - _G(V, xi, yi-1, zi)) - _x[17] = 0.5*(_G(V, xi+1, yi+1,zi) - _G(V, xi+1, yi-1, zi)) - _x[18] = 0.5*(_G(V, xi, yi+2,zi) - _G(V, xi, yi, zi)) - _x[19] = 0.5*(_G(V, xi+1, yi+2,zi) - _G(V, xi+1, yi, zi)) - _x[20] = 0.5*(_G(V, xi, yi+1,zi+1) - _G(V, xi, yi-1, zi+1)) - _x[21] = 0.5*(_G(V, xi+1, yi+1,zi+1) - _G(V, xi+1, yi-1, zi+1)) - _x[22] = 0.5*(_G(V, xi, yi+2,zi+1) - _G(V, xi, yi, zi+1)) - _x[23] = 0.5*(_G(V, xi+1, yi+2,zi+1) - _G(V, xi+1, yi, zi+1)) - _x[24] = 0.5*(_G(V, xi, yi,zi+1) - _G(V, xi, yi, zi-1)) - _x[25] = 0.5*(_G(V, xi+1, yi,zi+1) - _G(V, xi+1, yi, zi-1)) - _x[26] = 0.5*(_G(V, xi, yi+1,zi+1) - _G(V, xi, yi+1, zi-1)) - _x[27] = 0.5*(_G(V, xi+1, yi+1,zi+1) - _G(V, xi+1, yi+1, zi-1)) - _x[28] = 0.5*(_G(V, xi, yi,zi+2) - _G(V, xi, yi, zi)) - _x[29] = 0.5*(_G(V, xi+1, yi,zi+2) - _G(V, xi+1, yi, zi)) - _x[30] = 0.5*(_G(V, xi, yi+1,zi+2) - _G(V, xi, yi+1, zi)) - _x[31] = 0.5*(_G(V, xi+1, yi+1,zi+2) - _G(V, xi+1, yi+1, zi)) - - _x [32] = 0.25*(_G(V, xi+1, yi+1, zi) - _G(V, xi-1, yi+1, zi) - _G(V, xi+1, yi-1, zi) + _G(V, xi-1, yi-1, zi)) - _x [33] = 0.25*(_G(V, xi+2, yi+1, zi) - _G(V, xi, yi+1, zi) - _G(V, xi+2, yi-1, zi) + _G(V, xi, yi-1, zi)) - _x [34] = 0.25*(_G(V, xi+1, yi+2, zi) - _G(V, xi-1, yi+2, zi) - _G(V, xi+1, yi, zi) + _G(V, xi-1, yi, zi)) - _x [35] = 0.25*(_G(V, xi+2, yi+2, zi) - _G(V, xi, yi+2, zi) - _G(V, xi+2, yi, zi) + _G(V, xi, yi, zi)) - _x [36] = 0.25*(_G(V, xi+1, yi+1, zi+1) - _G(V, xi-1, yi+1, zi+1) - _G(V, xi+1, yi-1, zi+1) + _G(V, xi-1, yi-1, zi+1)) - _x [37] = 0.25*(_G(V, xi+2, yi+1, zi+1) - _G(V, xi, yi+1, zi+1) - _G(V, xi+2, yi-1, zi+1) + _G(V, xi, yi-1, zi+1)) - _x [38] = 0.25*(_G(V, xi+1, yi+2, zi+1) - _G(V, xi-1, yi+2, zi+1) - _G(V, xi+1, yi, zi+1) + _G(V, xi-1, yi, zi+1)) - _x [39] = 0.25*(_G(V, xi+2, yi+2, zi+1) - _G(V, xi, yi+2, zi+1) - _G(V, xi+2, yi, zi+1) + _G(V, xi, yi, zi+1)) - _x [40] = 0.25*(_G(V, xi+1, yi, zi+1) - _G(V, xi-1, yi, zi+1) - _G(V, xi+1, yi, zi-1) + _G(V, xi-1, yi, zi-1)) - _x [41] = 0.25*(_G(V, xi+2, yi, zi+1) - _G(V, xi, yi, zi+1) - _G(V, xi+2, yi, zi-1) + _G(V, xi, yi, zi-1)) - _x [42] = 0.25*(_G(V, xi+1, yi+1, zi+1) - _G(V, xi-1, yi+1, zi+1) - _G(V, xi+1, yi+1, zi-1) + _G(V, xi-1, yi+1, zi-1)) - _x [43] = 0.25*(_G(V, xi+2, yi+1, zi+1) - _G(V, xi, yi+1, zi+1) - _G(V, xi+2, yi+1, zi-1) + _G(V, xi, yi+1, zi-1)) - _x [44] = 0.25*(_G(V, xi+1, yi, zi+2) - _G(V, xi-1, yi, zi+2) - _G(V, xi+1, yi, zi) + _G(V, xi-1, yi, zi)) - _x [45] = 0.25*(_G(V, xi+2, yi, zi+2) - _G(V, xi, yi, zi+2) - _G(V, xi+2, yi, zi) + _G(V, xi, yi, zi)) - _x [46] = 0.25*(_G(V, xi+1, yi+1, zi+2) - _G(V, xi-1, yi+1, zi+2) - _G(V, xi+1, yi+1, zi) + _G(V, xi-1, yi+1, zi)) - _x [47] = 0.25*(_G(V, xi+2, yi+1, zi+2) - _G(V, xi, yi+1, zi+2) - _G(V, xi+2, yi+1, zi) + _G(V, xi, yi+1, zi)) - _x [48] = 0.25*(_G(V, xi, yi+1, zi+1) - _G(V, xi, yi-1, zi+1) - _G(V, xi, yi+1, zi-1) + _G(V, xi, yi-1, zi-1)) - _x [49] = 0.25*(_G(V, xi+1, yi+1, zi+1) - _G(V, xi+1, yi-1, zi+1) - _G(V, xi+1, yi+1, zi-1) + _G(V, xi+1, yi-1, zi-1)) - _x [50] = 0.25*(_G(V, xi, yi+2, zi+1) - _G(V, xi, yi, zi+1) - _G(V, xi, yi+2, zi-1) + _G(V, xi, yi, zi-1)) - _x [51] = 0.25*(_G(V, xi+1, yi+2, zi+1) - _G(V, xi+1, yi, zi+1) - _G(V, xi+1, yi+2, zi-1) + _G(V, xi+1, yi, zi-1)) - _x [52] = 0.25*(_G(V, xi, yi+1, zi+2) - _G(V, xi, yi-1, zi+2) - _G(V, xi, yi+1, zi) + _G(V, xi, yi-1, zi)) - _x [53] = 0.25*(_G(V, xi+1, yi+1, zi+2) - _G(V, xi+1, yi-1, zi+2) - _G(V, xi+1, yi+1, zi) + _G(V, xi+1, yi-1, zi)) - _x [54] = 0.25*(_G(V, xi, yi+2, zi+2) - _G(V, xi, yi, zi+2) - _G(V, xi, yi+2, zi) + _G(V, xi, yi, zi)) - _x [55] = 0.25*(_G(V, xi+1, yi+2, zi+2) - _G(V, xi+1, yi, zi+2) - _G(V, xi+1, yi+2, zi) + _G(V, xi+1, yi, zi)) - - _x[56] = 0.125*(_G(V, xi+1, yi+1, zi+1) - _G(V, xi-1, yi+1, zi+1) - _G(V, xi+1, yi-1, zi+1) + _G(V, xi-1, yi-1, zi+1) - _G(V, xi+1, yi+1, zi-1) + _G(V, xi-1,yi+1,zi-1)+_G(V, xi+1,yi-1,zi-1)-_G(V, xi-1,yi-1,zi-1)) - _x[57] = 0.125*(_G(V, xi+2, yi+1, zi+1) - _G(V, xi, yi+1, zi+1) - _G(V, xi+2, yi-1, zi+1) + _G(V, xi, yi-1, zi+1) - _G(V, xi+2, yi+1, zi-1) + _G(V, xi,yi+1,zi-1)+_G(V, xi+2,yi-1,zi-1)-_G(V, xi,yi-1,zi-1)) - _x[58] = 0.125*(_G(V, xi+1, yi+2, zi+1) - _G(V, xi-1, yi+2, zi+1) - _G(V, xi+1, yi, zi+1) + _G(V, xi-1, yi, zi+1) - _G(V, xi+1, yi+2, zi-1) + _G(V, xi-1,yi+2,zi-1)+_G(V, xi+1,yi,zi-1)-_G(V, xi-1,yi,zi-1)) - _x[59] = 0.125*(_G(V, xi+2, yi+2, zi+1) - _G(V, xi, yi+2, zi+1) - _G(V, xi+2, yi, zi+1) + _G(V, xi, yi, zi+1) - _G(V, xi+2, yi+2, zi-1) + _G(V, xi,yi+2,zi-1)+_G(V, xi+2,yi,zi-1)-_G(V, xi,yi,zi-1)) - _x[60] = 0.125*(_G(V, xi+1, yi+1, zi+2) - _G(V, xi-1, yi+1, zi+2) - _G(V, xi+1, yi-1, zi+2) + _G(V, xi-1, yi-1, zi+2) - _G(V, xi+1, yi+1, zi) + _G(V, xi-1,yi+1,zi)+_G(V, xi+1,yi-1,zi)-_G(V, xi-1,yi-1,zi)) - _x[61] = 0.125*(_G(V, xi+2, yi+1, zi+2) - _G(V, xi, yi+1, zi+2) - _G(V, xi+2, yi-1, zi+2) + _G(V, xi, yi-1, zi+2) - _G(V, xi+2, yi+1, zi) + _G(V, xi,yi+1,zi)+_G(V, xi+2,yi-1,zi)-_G(V, xi,yi-1,zi)) - _x[62] = 0.125*(_G(V, xi+1, yi+2, zi+2) - _G(V, xi-1, yi+2, zi+2) - _G(V, xi+1, yi, zi+2) + _G(V, xi-1, yi, zi+2) - _G(V, xi+1, yi+2, zi) + _G(V, xi-1,yi+2,zi)+_G(V, xi+1,yi,zi)-_G(V, xi-1,yi,zi)) - _x[63] = 0.125*(_G(V, xi+2, yi+2, zi+2) - _G(V, xi, yi+2, zi+2) - _G(V, xi+2, yi, zi+2) + _G(V, xi, yi, zi+2) - _G(V, xi+2, yi+2, zi) + _G(V, xi,yi+2,zi)+_G(V, xi+2,yi,zi)-_G(V, xi,yi,zi)) - - for j in prange(64): - coef[j] = 0.0 - for i in xrange(64): - coef[j] += (temp[j][i] * _x[i]) - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -cdef double tricub_interpolate(image_t[:, :, :] V, double x, double y, double z) nogil: - # From: Tricubic interpolation in three dimensions. Lekien and Marsden - cdef double[64] coef - cdef double result = 0.0 - calc_coef_tricub(V, x, y, z, coef) - - cdef int i, j, k - - cdef int xi = floor(x) - cdef int yi = floor(y) - cdef int zi = floor(z) - - for i in xrange(4): - for j in xrange(4): - for k in xrange(4): - result += (coef[i+4*j+16*k] * ((x-xi)**i) * ((y-yi)**j) * ((z-zi)**k)) - # return V[z, y, x] - # with gil: - # print result - return result - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -cdef double cubicInterpolate(double p[4], double x) nogil: - return p[1] + 0.5 * x*(p[2] - p[0] + x*(2.0*p[0] - 5.0*p[1] + 4.0*p[2] - p[3] + x*(3.0*(p[1] - p[2]) + p[3] - p[0]))) - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -cdef double bicubicInterpolate (double p[4][4], double x, double y) nogil: - cdef double arr[4] - arr[0] = cubicInterpolate(p[0], y) - arr[1] = cubicInterpolate(p[1], y) - arr[2] = cubicInterpolate(p[2], y) - arr[3] = cubicInterpolate(p[3], y) - return cubicInterpolate(arr, x) - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -cdef double tricubicInterpolate(image_t[:, :, :] V, double x, double y, double z) nogil: - # From http://www.paulinternet.nl/?page=bicubic - cdef double p[4][4][4] - - cdef int xi = floor(x) - cdef int yi = floor(y) - cdef int zi = floor(z) - - cdef int i, j, k - - for i in xrange(4): - for j in xrange(4): - for k in xrange(4): - p[i][j][k] = _G(V, xi + i -1, yi + j -1, zi + k - 1) - - cdef double arr[4] - arr[0] = bicubicInterpolate(p[0], y-yi, z-zi) - arr[1] = bicubicInterpolate(p[1], y-yi, z-zi) - arr[2] = bicubicInterpolate(p[2], y-yi, z-zi) - arr[3] = bicubicInterpolate(p[3], y-yi, z-zi) - return cubicInterpolate(arr, x-xi) - - -def tricub_interpolate_py(image_t[:, :, :] V, double x, double y, double z): - return tricub_interpolate(V, x, y, z) - -def tricub_interpolate2_py(image_t[:, :, :] V, double x, double y, double z): - return tricubicInterpolate(V, x, y, z) - -def trilin_interpolate_py(image_t[:, :, :] V, double x, double y, double z): - return interpolate(V, x, y, z) diff --git a/invesalius/data/mask.py b/invesalius/data/mask.py index 2bae7b4..40df971 100644 --- a/invesalius/data/mask.py +++ b/invesalius/data/mask.py @@ -30,7 +30,7 @@ import invesalius.constants as const import invesalius.data.imagedata_utils as iu import invesalius.session as ses -from . import floodfill +from invesalius_cy import floodfill from wx.lib.pubsub import pub as Publisher from scipy import ndimage diff --git a/invesalius/data/mips.pyx b/invesalius/data/mips.pyx deleted file mode 100644 index 6b0e0a2..0000000 --- a/invesalius/data/mips.pyx +++ /dev/null @@ -1,379 +0,0 @@ -#http://en.wikipedia.org/wiki/Local_maximum_intensity_projection -import numpy as np -cimport numpy as np -cimport cython - -from libc.math cimport floor, ceil, sqrt, fabs -from cython.parallel import prange - -DTYPE = np.uint8 -ctypedef np.uint8_t DTYPE_t - -DTYPE16 = np.int16 -ctypedef np.int16_t DTYPE16_t - -DTYPEF32 = np.float32 -ctypedef np.float32_t DTYPEF32_t - -@cython.boundscheck(False) # turn of bounds-checking for entire function -def lmip(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t tmin, - DTYPE16_t tmax, np.ndarray[DTYPE16_t, ndim=2] out): - cdef DTYPE16_t max - cdef int start - cdef int sz = image.shape[0] - cdef int sy = image.shape[1] - cdef int sx = image.shape[2] - - # AXIAL - if axis == 0: - for x in xrange(sx): - for y in xrange(sy): - max = image[0, y, x] - if max >= tmin and max <= tmax: - start = 1 - else: - start = 0 - for z in xrange(sz): - if image[z, y, x] > max: - max = image[z, y, x] - - elif image[z, y, x] < max and start: - break - - if image[z, y, x] >= tmin and image[z, y, x] <= tmax: - start = 1 - - out[y, x] = max - - #CORONAL - elif axis == 1: - for z in xrange(sz): - for x in xrange(sx): - max = image[z, 0, x] - if max >= tmin and max <= tmax: - start = 1 - else: - start = 0 - for y in xrange(sy): - if image[z, y, x] > max: - max = image[z, y, x] - - elif image[z, y, x] < max and start: - break - - if image[z, y, x] >= tmin and image[z, y, x] <= tmax: - start = 1 - - out[z, x] = max - - #CORONAL - elif axis == 2: - for z in xrange(sz): - for y in xrange(sy): - max = image[z, y, 0] - if max >= tmin and max <= tmax: - start = 1 - else: - start = 0 - for x in xrange(sx): - if image[z, y, x] > max: - max = image[z, y, x] - - elif image[z, y, x] < max and start: - break - - if image[z, y, x] >= tmin and image[z, y, x] <= tmax: - start = 1 - - out[z, y] = max - - -cdef DTYPE16_t get_colour(DTYPE16_t vl, DTYPE16_t wl, DTYPE16_t ww): - cdef DTYPE16_t out_colour - cdef DTYPE16_t min_value = wl - (ww / 2) - cdef DTYPE16_t max_value = wl + (ww / 2) - if vl < min_value: - out_colour = min_value - elif vl > max_value: - out_colour = max_value - else: - out_colour = vl - - return out_colour - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -cdef float get_opacity(DTYPE16_t vl, DTYPE16_t wl, DTYPE16_t ww) nogil: - cdef float out_opacity - cdef DTYPE16_t min_value = wl - (ww / 2) - cdef DTYPE16_t max_value = wl + (ww / 2) - if vl < min_value: - out_opacity = 0.0 - elif vl > max_value: - out_opacity = 1.0 - else: - out_opacity = 1.0/(max_value - min_value) * (vl - min_value) - - return out_opacity - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -cdef float get_opacity_f32(DTYPEF32_t vl, DTYPE16_t wl, DTYPE16_t ww) nogil: - cdef float out_opacity - cdef DTYPE16_t min_value = wl - (ww / 2) - cdef DTYPE16_t max_value = wl + (ww / 2) - if vl < min_value: - out_opacity = 0.0 - elif vl > max_value: - out_opacity = 1.0 - else: - out_opacity = 1.0/(max_value - min_value) * (vl - min_value) - - return out_opacity - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -def mida(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t wl, - DTYPE16_t ww, np.ndarray[DTYPE16_t, ndim=2] out): - cdef int sz = image.shape[0] - cdef int sy = image.shape[1] - cdef int sx = image.shape[2] - - cdef DTYPE16_t min = image.min() - cdef DTYPE16_t max = image.max() - cdef DTYPE16_t vl - - cdef DTYPE16_t min_value = wl - (ww / 2) - cdef DTYPE16_t max_value = wl + (ww / 2) - - cdef float fmax=0.0 - cdef float fpi - cdef float dl - cdef float bt - - cdef float alpha - cdef float alpha_p = 0.0 - cdef float colour - cdef float colour_p = 0 - - cdef int x, y, z - - # AXIAL - if axis == 0: - for x in prange(sx, nogil=True): - for y in xrange(sy): - fmax = 0.0 - alpha_p = 0.0 - colour_p = 0.0 - for z in xrange(sz): - vl = image[z, y, x] - fpi = 1.0/(max - min) * (vl - min) - if fpi > fmax: - dl = fpi - fmax - fmax = fpi - else: - dl = 0.0 - - bt = 1.0 - dl - - colour = fpi - alpha = get_opacity(vl, wl, ww) - colour = (bt * colour_p) + (1 - bt * alpha_p) * colour * alpha - alpha = (bt * alpha_p) + (1 - bt * alpha_p) * alpha - - colour_p = colour - alpha_p = alpha - - if alpha >= 1.0: - break - - - #out[y, x] = ((max_value - min_value) * colour + min_value) - out[y, x] = ((max - min) * colour + min) - - - #CORONAL - elif axis == 1: - for z in prange(sz, nogil=True): - for x in xrange(sx): - fmax = 0.0 - alpha_p = 0.0 - colour_p = 0.0 - for y in xrange(sy): - vl = image[z, y, x] - fpi = 1.0/(max - min) * (vl - min) - if fpi > fmax: - dl = fpi - fmax - fmax = fpi - else: - dl = 0.0 - - bt = 1.0 - dl - - colour = fpi - alpha = get_opacity(vl, wl, ww) - colour = (bt * colour_p) + (1 - bt * alpha_p) * colour * alpha - alpha = (bt * alpha_p) + (1 - bt * alpha_p) * alpha - - colour_p = colour - alpha_p = alpha - - if alpha >= 1.0: - break - - out[z, x] = ((max - min) * colour + min) - - #AXIAL - elif axis == 2: - for z in prange(sz, nogil=True): - for y in xrange(sy): - fmax = 0.0 - alpha_p = 0.0 - colour_p = 0.0 - for x in xrange(sx): - vl = image[z, y, x] - fpi = 1.0/(max - min) * (vl - min) - if fpi > fmax: - dl = fpi - fmax - fmax = fpi - else: - dl = 0.0 - - bt = 1.0 - dl - - colour = fpi - alpha = get_opacity(vl, wl, ww) - colour = (bt * colour_p) + (1 - bt * alpha_p) * colour * alpha - alpha = (bt * alpha_p) + (1 - bt * alpha_p) * alpha - - colour_p = colour - alpha_p = alpha - - if alpha >= 1.0: - break - - out[z, y] = ((max - min) * colour + min) - - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -cdef inline void finite_difference(DTYPE16_t[:, :, :] image, - int x, int y, int z, float h, float *g) nogil: - cdef int px, py, pz, fx, fy, fz - - cdef int sz = image.shape[0] - cdef int sy = image.shape[1] - cdef int sx = image.shape[2] - - cdef float gx, gy, gz - - if x == 0: - px = 0 - fx = 1 - elif x == sx - 1: - px = x - 1 - fx = x - else: - px = x - 1 - fx = x + 1 - - if y == 0: - py = 0 - fy = 1 - elif y == sy - 1: - py = y - 1 - fy = y - else: - py = y - 1 - fy = y + 1 - - if z == 0: - pz = 0 - fz = 1 - elif z == sz - 1: - pz = z - 1 - fz = z - else: - pz = z - 1 - fz = z + 1 - - gx = (image[z, y, fx] - image[z, y, px]) / (2*h) - gy = (image[z, fy, x] - image[z, py, x]) / (2*h) - gz = (image[fz, y, x] - image[pz, y, x]) / (2*h) - - g[0] = gx - g[1] = gy - g[2] = gz - - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -cdef inline float calc_fcm_itensity(DTYPE16_t[:, :, :] image, - int x, int y, int z, float n, float* dir) nogil: - cdef float g[3] - finite_difference(image, x, y, z, 1.0, g) - cdef float gm = sqrt(g[0]*g[0] + g[1]*g[1] + g[2]*g[2]) - cdef float d = g[0]*dir[0] + g[1]*dir[1] + g[2]*dir[2] - cdef float sf = (1.0 - fabs(d/gm))**n - #alpha = get_opacity_f32(gm, wl, ww) - cdef float vl = gm * sf - return vl - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -def fast_countour_mip(np.ndarray[DTYPE16_t, ndim=3] image, - float n, - int axis, - DTYPE16_t wl, DTYPE16_t ww, - int tmip, - np.ndarray[DTYPE16_t, ndim=2] out): - cdef int sz = image.shape[0] - cdef int sy = image.shape[1] - cdef int sx = image.shape[2] - cdef float gm - cdef float alpha - cdef float sf - cdef float d - - cdef float* g - cdef float* dir = [ 0, 0, 0 ] - - cdef DTYPE16_t[:, :, :] vimage = image - cdef np.ndarray[DTYPE16_t, ndim=3] tmp = np.empty_like(image) - - cdef DTYPE16_t min = image.min() - cdef DTYPE16_t max = image.max() - cdef float fmin = min - cdef float fmax = max - cdef float vl - cdef DTYPE16_t V - - cdef int x, y, z - - if axis == 0: - dir[2] = 1.0 - elif axis == 1: - dir[1] = 1.0 - elif axis == 2: - dir[0] = 1.0 - - for z in prange(sz, nogil=True): - for y in range(sy): - for x in range(sx): - vl = calc_fcm_itensity(vimage, x, y, z, n, dir) - tmp[z, y, x] = vl - - cdef DTYPE16_t tmin = tmp.min() - cdef DTYPE16_t tmax = tmp.max() - - #tmp = ((max - min)/(tmax - tmin)) * (tmp - tmin) + min - - if tmip == 0: - out[:] = tmp.max(axis) - elif tmip == 1: - lmip(tmp, axis, 700, 3033, out) - elif tmip == 2: - mida(tmp, axis, wl, ww, out) diff --git a/invesalius/data/slice_.py b/invesalius/data/slice_.py index 33096c8..e53235f 100644 --- a/invesalius/data/slice_.py +++ b/invesalius/data/slice_.py @@ -28,13 +28,13 @@ from wx.lib.pubsub import pub as Publisher import invesalius.constants as const import invesalius.data.converters as converters import invesalius.data.imagedata_utils as iu -import invesalius.data.transformations as transformations import invesalius.session as ses import invesalius.style as st import invesalius.utils as utils -from invesalius.data import mips, transforms +from invesalius.data import transformations from invesalius.data.mask import Mask from invesalius.project import Project +from invesalius_cy import mips, transforms OTHER = 0 PLIST = 1 diff --git a/invesalius/data/styles.py b/invesalius/data/styles.py index 5cfadf0..1fbd6c0 100644 --- a/invesalius/data/styles.py +++ b/invesalius/data/styles.py @@ -46,7 +46,7 @@ import invesalius.utils as utils from invesalius.data.measures import (CircleDensityMeasure, MeasureData, PolygonDensityMeasure) -from . import floodfill +from invesalius_cy import floodfill ORIENTATIONS = { "AXIAL": const.AXIAL, diff --git a/invesalius/data/surface.py b/invesalius/data/surface.py index 5e5fee6..d30fe90 100644 --- a/invesalius/data/surface.py +++ b/invesalius/data/surface.py @@ -59,9 +59,9 @@ import invesalius.utils as utl import invesalius.data.vtk_utils as vu from invesalius.gui import dialogs +from invesalius_cy import cy_mesh -from invesalius.data import cy_mesh -# TODO: Verificar ReleaseDataFlagOn and SetSource +# TODO: Verificar ReleaseDataFlagOn and SetSource class Surface(): diff --git a/invesalius/data/surface_process.py b/invesalius/data/surface_process.py index 19a0ed3..9ad896d 100644 --- a/invesalius/data/surface_process.py +++ b/invesalius/data/surface_process.py @@ -13,7 +13,7 @@ import vtk import invesalius.i18n as i18n import invesalius.data.converters as converters -from invesalius.data import cy_mesh +from invesalius_cy import cy_mesh import weakref from scipy import ndimage diff --git a/invesalius/data/transforms.pyx b/invesalius/data/transforms.pyx deleted file mode 100644 index 5c746d2..0000000 --- a/invesalius/data/transforms.pyx +++ /dev/null @@ -1,174 +0,0 @@ -import numpy as np -cimport numpy as np -cimport cython - -from .cy_my_types cimport image_t -from .interpolation cimport interpolate, tricub_interpolate, tricubicInterpolate, lanczos3, nearest_neighbour_interp - -from libc.math cimport floor, ceil, sqrt, fabs, round -from cython.parallel import prange - -ctypedef double (*interp_function)(image_t[:, :, :], double, double, double) nogil - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -cdef void mul_mat4_vec4(double[:, :] M, - double* coord, - double* out) nogil: - - out[0] = coord[0] * M[0, 0] + coord[1] * M[0, 1] + coord[2] * M[0, 2] + coord[3] * M[0, 3] - out[1] = coord[0] * M[1, 0] + coord[1] * M[1, 1] + coord[2] * M[1, 2] + coord[3] * M[1, 3] - out[2] = coord[0] * M[2, 0] + coord[1] * M[2, 1] + coord[2] * M[2, 2] + coord[3] * M[2, 3] - out[3] = coord[0] * M[3, 0] + coord[1] * M[3, 1] + coord[2] * M[3, 2] + coord[3] * M[3, 3] - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -cdef image_t coord_transform(image_t[:, :, :] volume, double[:, :] M, int x, int y, int z, double sx, double sy, double sz, interp_function f_interp, image_t cval) nogil: - - cdef double coord[4] - coord[0] = z*sz - coord[1] = y*sy - coord[2] = x*sx - coord[3] = 1.0 - - cdef double _ncoord[4] - _ncoord[3] = 1 - # _ncoord[:] = [0.0, 0.0, 0.0, 1.0] - - cdef unsigned int dz, dy, dx - dz = volume.shape[0] - dy = volume.shape[1] - dx = volume.shape[2] - - - mul_mat4_vec4(M, coord, _ncoord) - - cdef double nz, ny, nx - nz = (_ncoord[0]/_ncoord[3])/sz - ny = (_ncoord[1]/_ncoord[3])/sy - nx = (_ncoord[2]/_ncoord[3])/sx - - cdef double v - - if 0 <= nz <= (dz-1) and 0 <= ny <= (dy-1) and 0 <= nx <= (dx-1): - return f_interp(volume, nx, ny, nz) - # if minterpol == 0: - # return nearest_neighbour_interp(volume, nx, ny, nz) - # elif minterpol == 1: - # return interpolate(volume, nx, ny, nz) - # elif minterpol == 2: - # v = tricubicInterpolate(volume, nx, ny, nz) - # if (v < cval): - # v = cval - # return v - # else: - # v = lanczos3(volume, nx, ny, nz) - # if (v < cval): - # v = cval - # return v - else: - return cval - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -def apply_view_matrix_transform(image_t[:, :, :] volume, - spacing, - double[:, :] M, - unsigned int n, str orientation, - int minterpol, - image_t cval, - image_t[:, :, :] out): - - cdef int dz, dy, dx - cdef int z, y, x - dz = volume.shape[0] - dy = volume.shape[1] - dx = volume.shape[2] - - cdef unsigned int odz, ody, odx - odz = out.shape[0] - ody = out.shape[1] - odx = out.shape[2] - - cdef unsigned int count = 0 - - cdef double sx, sy, sz - sx = spacing[0] - sy = spacing[1] - sz = spacing[2] - - cdef interp_function f_interp - - if minterpol == 0: - f_interp = nearest_neighbour_interp - elif minterpol == 1: - f_interp = interpolate - elif minterpol == 2: - f_interp = tricubicInterpolate - else: - f_interp = lanczos3 - - if orientation == 'AXIAL': - for z in xrange(n, n+odz): - for y in prange(dy, nogil=True): - for x in xrange(dx): - out[count, y, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval) - count += 1 - - elif orientation == 'CORONAL': - for y in xrange(n, n+ody): - for z in prange(dz, nogil=True): - for x in xrange(dx): - out[z, count, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval) - count += 1 - - elif orientation == 'SAGITAL': - for x in xrange(n, n+odx): - for z in prange(dz, nogil=True): - for y in xrange(dy): - out[z, y, count] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval) - count += 1 - - -@cython.boundscheck(False) # turn of bounds-checking for entire function -@cython.cdivision(True) -@cython.wraparound(False) -def convolve_non_zero(image_t[:, :, :] volume, - image_t[:, :, :] kernel, - image_t cval): - cdef Py_ssize_t x, y, z, sx, sy, sz, kx, ky, kz, skx, sky, skz, i, j, k - cdef image_t v - - cdef image_t[:, :, :] out = np.zeros_like(volume) - - sz = volume.shape[0] - sy = volume.shape[1] - sx = volume.shape[2] - - skz = kernel.shape[0] - sky = kernel.shape[1] - skx = kernel.shape[2] - - for z in prange(sz, nogil=True): - for y in xrange(sy): - for x in xrange(sx): - if volume[z, y, x] != 0: - for k in xrange(skz): - kz = z - skz // 2 + k - for j in xrange(sky): - ky = y - sky // 2 + j - for i in xrange(skx): - kx = x - skx // 2 + i - - if 0 <= kz < sz and 0 <= ky < sy and 0 <= kx < sx: - v = volume[kz, ky, kx] - else: - v = cval - - out[z, y, x] += (v * kernel[k, j, i]) - return np.asarray(out) diff --git a/invesalius_cy/cy_mesh.pyx b/invesalius_cy/cy_mesh.pyx new file mode 100644 index 0000000..6d074f8 --- /dev/null +++ b/invesalius_cy/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 maximum 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 minimum 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 minimum 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_cy/cy_my_types.pxd b/invesalius_cy/cy_my_types.pxd new file mode 100644 index 0000000..13d6c6c --- /dev/null +++ b/invesalius_cy/cy_my_types.pxd @@ -0,0 +1,21 @@ +import numpy as np +cimport numpy as np +cimport cython + +# ctypedef np.uint16_t image_t + +ctypedef fused image_t: + np.float64_t + np.int16_t + np.uint8_t + +ctypedef np.uint8_t mask_t + +ctypedef np.float32_t vertex_t +ctypedef np.float32_t normal_t + +# To compile in windows 64 +IF UNAME_MACHINE == 'AMD64': + ctypedef np.int64_t vertex_id_t +ELSE: + ctypedef np.int_t vertex_id_t diff --git a/invesalius_cy/floodfill.pyx b/invesalius_cy/floodfill.pyx new file mode 100644 index 0000000..6636cc9 --- /dev/null +++ b/invesalius_cy/floodfill.pyx @@ -0,0 +1,274 @@ +import numpy as np +cimport numpy as np +cimport cython + +from collections import deque + +from cython.parallel import prange +from libc.math cimport floor, ceil +from libcpp cimport bool +from libcpp.deque cimport deque as cdeque +from libcpp.vector cimport vector + +from cy_my_types cimport image_t, mask_t + +cdef struct s_coord: + int x + int y + int z + +ctypedef s_coord coord + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.wraparound(False) +@cython.nonecheck(False) +@cython.cdivision(True) +cdef inline void append_queue(cdeque[int]& stack, int x, int y, int z, int d, int h, int w) nogil: + stack.push_back(z*h*w + y*w + x) + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.wraparound(False) +@cython.nonecheck(False) +@cython.cdivision(True) +cdef inline void pop_queue(cdeque[int]& stack, int* x, int* y, int* z, int d, int h, int w) nogil: + cdef int i = stack.front() + stack.pop_front() + x[0] = i % w + y[0] = (i / w) % h + z[0] = i / (h * w) + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +def floodfill(np.ndarray[image_t, ndim=3] data, int i, int j, int k, int v, int fill, np.ndarray[mask_t, ndim=3] out): + + cdef int to_return = 0 + if out is None: + out = np.zeros_like(data) + to_return = 1 + + cdef int x, y, z + cdef int w, h, d + + d = data.shape[0] + h = data.shape[1] + w = data.shape[2] + + stack = [(i, j, k), ] + out[k, j, i] = fill + + while stack: + x, y, z = stack.pop() + + if z + 1 < d and data[z + 1, y, x] == v and out[z + 1, y, x] != fill: + out[z + 1, y, x] = fill + stack.append((x, y, z + 1)) + + if z - 1 >= 0 and data[z - 1, y, x] == v and out[z - 1, y, x] != fill: + out[z - 1, y, x] = fill + stack.append((x, y, z - 1)) + + if y + 1 < h and data[z, y + 1, x] == v and out[z, y + 1, x] != fill: + out[z, y + 1, x] = fill + stack.append((x, y + 1, z)) + + if y - 1 >= 0 and data[z, y - 1, x] == v and out[z, y - 1, x] != fill: + out[z, y - 1, x] = fill + stack.append((x, y - 1, z)) + + if x + 1 < w and data[z, y, x + 1] == v and out[z, y, x + 1] != fill: + out[z, y, x + 1] = fill + stack.append((x + 1, y, z)) + + if x - 1 >= 0 and data[z, y, x - 1] == v and out[z, y, x - 1] != fill: + out[z, y, x - 1] = fill + stack.append((x - 1, y, z)) + + if to_return: + return out + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.wraparound(False) +@cython.nonecheck(False) +def floodfill_threshold(np.ndarray[image_t, ndim=3] data, list seeds, int t0, int t1, int fill, np.ndarray[mask_t, ndim=3] strct, np.ndarray[mask_t, ndim=3] out): + + cdef int to_return = 0 + if out is None: + out = np.zeros_like(data) + to_return = 1 + + cdef int x, y, z + cdef int dx, dy, dz + cdef int odx, ody, odz + cdef int xo, yo, zo + cdef int i, j, k + cdef int offset_x, offset_y, offset_z + + dz = data.shape[0] + dy = data.shape[1] + dx = data.shape[2] + + odz = strct.shape[0] + ody = strct.shape[1] + odx = strct.shape[2] + + cdef cdeque[coord] stack + cdef coord c + + offset_z = odz / 2 + offset_y = ody / 2 + offset_x = odx / 2 + + for i, j, k in seeds: + if data[k, j, i] >= t0 and data[k, j, i] <= t1: + c.x = i + c.y = j + c.z = k + stack.push_back(c) + out[k, j, i] = fill + + with nogil: + while stack.size(): + c = stack.back() + stack.pop_back() + + x = c.x + y = c.y + z = c.z + + out[z, y, x] = fill + + for k in xrange(odz): + zo = z + k - offset_z + for j in xrange(ody): + yo = y + j - offset_y + for i in xrange(odx): + if strct[k, j, i]: + xo = x + i - offset_x + if 0 <= xo < dx and 0 <= yo < dy and 0 <= zo < dz and out[zo, yo, xo] != fill and t0 <= data[zo, yo, xo] <= t1: + out[zo, yo, xo] = fill + c.x = xo + c.y = yo + c.z = zo + stack.push_back(c) + + if to_return: + return out + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.wraparound(False) +@cython.nonecheck(False) +def floodfill_auto_threshold(np.ndarray[image_t, ndim=3] data, list seeds, float p, int fill, np.ndarray[mask_t, ndim=3] out): + + cdef int to_return = 0 + if out is None: + out = np.zeros_like(data) + to_return = 1 + + cdef cdeque[int] stack + cdef int x, y, z + cdef int w, h, d + cdef int xo, yo, zo + cdef int t0, t1 + + cdef int i, j, k + + d = data.shape[0] + h = data.shape[1] + w = data.shape[2] + + + # stack = deque() + + x = 0 + y = 0 + z = 0 + + + for i, j, k in seeds: + append_queue(stack, i, j, k, d, h, w) + out[k, j, i] = fill + print i, j, k, d, h, w + + with nogil: + while stack.size(): + pop_queue(stack, &x, &y, &z, d, h, w) + + # print x, y, z, d, h, w + + xo = x + yo = y + zo = z + + t0 = ceil(data[z, y, x] * (1 - p)) + t1 = floor(data[z, y, x] * (1 + p)) + + if z + 1 < d and data[z + 1, y, x] >= t0 and data[z + 1, y, x] <= t1 and out[zo + 1, yo, xo] != fill: + out[zo + 1, yo, xo] = fill + append_queue(stack, x, y, z+1, d, h, w) + + if z - 1 >= 0 and data[z - 1, y, x] >= t0 and data[z - 1, y, x] <= t1 and out[zo - 1, yo, xo] != fill: + out[zo - 1, yo, xo] = fill + append_queue(stack, x, y, z-1, d, h, w) + + if y + 1 < h and data[z, y + 1, x] >= t0 and data[z, y + 1, x] <= t1 and out[zo, yo + 1, xo] != fill: + out[zo, yo + 1, xo] = fill + append_queue(stack, x, y+1, z, d, h, w) + + if y - 1 >= 0 and data[z, y - 1, x] >= t0 and data[z, y - 1, x] <= t1 and out[zo, yo - 1, xo] != fill: + out[zo, yo - 1, xo] = fill + append_queue(stack, x, y-1, z, d, h, w) + + if x + 1 < w and data[z, y, x + 1] >= t0 and data[z, y, x + 1] <= t1 and out[zo, yo, xo + 1] != fill: + out[zo, yo, xo + 1] = fill + append_queue(stack, x+1, y, z, d, h, w) + + if x - 1 >= 0 and data[z, y, x - 1] >= t0 and data[z, y, x - 1] <= t1 and out[zo, yo, xo - 1] != fill: + out[zo, yo, xo - 1] = fill + append_queue(stack, x-1, y, z, d, h, w) + + if to_return: + return out + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.nonecheck(False) +def fill_holes_automatically(np.ndarray[mask_t, ndim=3] mask, np.ndarray[np.uint16_t, ndim=3] labels, unsigned int nlabels, unsigned int max_size): + """ + Fill mask holes automatically. The hole must <= max_size. Return True if any hole were filled. + """ + cdef np.ndarray[np.uint32_t, ndim=1] sizes = np.zeros(shape=(nlabels + 1), dtype=np.uint32) + cdef int x, y, z + cdef int dx, dy, dz + cdef int i + + cdef bool modified = False + + dz = mask.shape[0] + dy = mask.shape[1] + dx = mask.shape[2] + + for z in xrange(dz): + for y in xrange(dy): + for x in xrange(dx): + sizes[labels[z, y, x]] += 1 + + #Checking if any hole will be filled + for i in xrange(nlabels + 1): + if sizes[i] <= max_size: + modified = True + + if not modified: + return 0 + + for z in prange(dz, nogil=True): + for y in xrange(dy): + for x in xrange(dx): + if sizes[labels[z, y, x]] <= max_size: + mask[z, y, x] = 254 + + return modified diff --git a/invesalius_cy/interpolation.pxd b/invesalius_cy/interpolation.pxd new file mode 100644 index 0000000..b63cc41 --- /dev/null +++ b/invesalius_cy/interpolation.pxd @@ -0,0 +1,8 @@ +from .cy_my_types cimport image_t + +cdef double interpolate(image_t[:, :, :], double, double, double) nogil +cdef double tricub_interpolate(image_t[:, :, :], double, double, double) nogil +cdef double tricubicInterpolate (image_t[:, :, :], double, double, double) nogil +cdef double lanczos3 (image_t[:, :, :], double, double, double) nogil + +cdef double nearest_neighbour_interp(image_t[:, :, :], double, double, double) nogil diff --git a/invesalius_cy/interpolation.pyx b/invesalius_cy/interpolation.pyx new file mode 100644 index 0000000..9f44648 --- /dev/null +++ b/invesalius_cy/interpolation.pyx @@ -0,0 +1,392 @@ +# from interpolation cimport interpolate + +import numpy as np +cimport numpy as np +cimport cython + +from libc.math cimport floor, ceil, sqrt, fabs, sin, M_PI +from cython.parallel import prange + +DEF LANCZOS_A = 4 +DEF SIZE_LANCZOS_TMP = LANCZOS_A * 2 - 1 + +cdef double[64][64] temp = [ + [ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 2, -2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 9, -9,-9, 9, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [-6, 6, 6,-6, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [-6, 6, 6,-6, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 4, -4,-4, 4, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0], + [-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 9, -9, 0, 0,-9, 9, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [-6, 6, 0, 0, 6,-6, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9, 0, 0,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0], + [ 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0], + [-27, 27,27,-27,27,-27,-27,27,-18,-9,18, 9,18, 9,-18,-9,-18,18,-9, 9,18,-18, 9,-9,-18,18,18,-18,-9, 9, 9,-9,-12,-6,-6,-3,12, 6, 6, 3,-12,-6,12, 6,-6,-3, 6, 3,-12,12,-6, 6,-6, 6,-3, 3,-8,-4,-4,-2,-4,-2,-2,-1], + [18, -18,-18,18,-18,18,18,-18, 9, 9,-9,-9,-9,-9, 9, 9,12,-12, 6,-6,-12,12,-6, 6,12,-12,-12,12, 6,-6,-6, 6, 6, 6, 3, 3,-6,-6,-3,-3, 6, 6,-6,-6, 3, 3,-3,-3, 8,-8, 4,-4, 4,-4, 2,-2, 4, 4, 2, 2, 2, 2, 1, 1], + [-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0], + [18, -18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6, 9,-9, 9,-9,-9, 9,-9, 9,12,-12,-12,12, 6,-6,-6, 6, 6, 3, 6, 3,-6,-3,-6,-3, 8, 4,-8,-4, 4, 2,-4,-2, 6,-6, 6,-6, 3,-3, 3,-3, 4, 2, 4, 2, 2, 1, 2, 1], + [-12, 12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-6, 6,-6, 6, 6,-6, 6,-6,-8, 8, 8,-8,-4, 4, 4,-4,-3,-3,-3,-3, 3, 3, 3, 3,-4,-4, 4, 4,-2,-2, 2, 2,-4, 4,-4, 4,-2, 2,-2, 2,-2,-2,-2,-2,-1,-1,-1,-1], + [ 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [-6, 6, 0, 0, 6,-6, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 4, -4, 0, 0,-4, 4, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4, 0, 0,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0], + [-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0], + [18, -18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6,12,-12, 6,-6,-12,12,-6, 6, 9,-9,-9, 9, 9,-9,-9, 9, 8, 4, 4, 2,-8,-4,-4,-2, 6, 3,-6,-3, 6, 3,-6,-3, 6,-6, 3,-3, 6,-6, 3,-3, 4, 2, 2, 1, 4, 2, 2, 1], + [-12, 12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-8, 8,-4, 4, 8,-8, 4,-4,-6, 6, 6,-6,-6, 6, 6,-6,-4,-4,-2,-2, 4, 4, 2, 2,-3,-3, 3, 3,-3,-3, 3, 3,-4, 4,-2, 2,-4, 4,-2, 2,-2,-2,-1,-1,-2,-2,-1,-1], + [ 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [ 0, 0, 0, 0, 0, 0, 0, 0, 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0], + [-12, 12,12,-12,12,-12,-12,12,-8,-4, 8, 4, 8, 4,-8,-4,-6, 6,-6, 6, 6,-6, 6,-6,-6, 6, 6,-6,-6, 6, 6,-6,-4,-2,-4,-2, 4, 2, 4, 2,-4,-2, 4, 2,-4,-2, 4, 2,-3, 3,-3, 3,-3, 3,-3, 3,-2,-1,-2,-1,-2,-1,-2,-1], + [ 8, -8,-8, 8,-8, 8, 8,-8, 4, 4,-4,-4,-4,-4, 4, 4, 4,-4, 4,-4,-4, 4,-4, 4, 4,-4,-4, 4, 4,-4,-4, 4, 2, 2, 2, 2,-2,-2,-2,-2, 2, 2,-2,-2, 2, 2,-2,-2, 2,-2, 2,-2, 2,-2, 2,-2, 1, 1, 1, 1, 1, 1, 1, 1] +] + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef double nearest_neighbour_interp(image_t[:, :, :] V, double x, double y, double z) nogil: + return V[(z), (y), (x)] + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef double interpolate(image_t[:, :, :] V, double x, double y, double z) nogil: + cdef double xd, yd, zd + cdef double c00, c10, c01, c11 + cdef double c0, c1 + cdef double c + + cdef int x0 = floor(x) + cdef int x1 = x0 + 1 + + cdef int y0 = floor(y) + cdef int y1 = y0 + 1 + + cdef int z0 = floor(z) + cdef int z1 = z0 + 1 + + if x0 == x1: + xd = 1.0 + else: + xd = (x - x0) / (x1 - x0) + + if y0 == y1: + yd = 1.0 + else: + yd = (y - y0) / (y1 - y0) + + if z0 == z1: + zd = 1.0 + else: + zd = (z - z0) / (z1 - z0) + + c00 = _G(V, x0, y0, z0)*(1 - xd) + _G(V, x1, y0, z0)*xd + c10 = _G(V, x0, y1, z0)*(1 - xd) + _G(V, x1, y1, z0)*xd + c01 = _G(V, x0, y0, z1)*(1 - xd) + _G(V, x1, y0, z1)*xd + c11 = _G(V, x0, y1, z1)*(1 - xd) + _G(V, x1, y1, z1)*xd + + c0 = c00*(1 - yd) + c10*yd + c1 = c01*(1 - yd) + c11*yd + + c = c0*(1 - zd) + c1*zd + + return c + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef inline double lanczos3_L(double x, int a) nogil: + if x == 0: + return 1.0 + elif -a <= x < a: + return (a * sin(M_PI * x) * sin(M_PI * (x / a)))/(M_PI**2 * x**2) + else: + return 0.0 + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef double lanczos3(image_t[:, :, :] V, double x, double y, double z) nogil: + cdef int a = LANCZOS_A + + cdef int xd = floor(x) + cdef int yd = floor(y) + cdef int zd = floor(z) + + cdef int xi = xd - a + 1 + cdef int xf = xd + a + + cdef int yi = yd - a + 1 + cdef int yf = yd + a + + cdef int zi = zd - a + 1 + cdef int zf = zd + a + + cdef double lx = 0.0 + cdef double ly = 0.0 + cdef double lz = 0.0 + + cdef double[SIZE_LANCZOS_TMP][SIZE_LANCZOS_TMP] temp_x + cdef double[SIZE_LANCZOS_TMP] temp_y + + cdef int i, j, k + cdef int m, n, o + + m = 0 + for k in xrange(zi, zf): + n = 0 + for j in xrange(yi, yf): + lx = 0 + for i in xrange(xi, xf): + lx += _G(V, i, j, k) * lanczos3_L(x - i, a) + temp_x[m][n] = lx + n += 1 + m += 1 + + m = 0 + for k in xrange(zi, zf): + n = 0 + ly = 0 + for j in xrange(yi, yf): + ly += temp_x[m][n] * lanczos3_L(y - j, a) + n += 1 + temp_y[m] = ly + m += 1 + + m = 0 + for k in xrange(zi, zf): + lz += temp_y[m] * lanczos3_L(z - k, a) + m += 1 + + return lz + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef image_t _G(image_t[:, :, :] V, int x, int y, int z) nogil: + cdef int dz, dy, dx + dz = V.shape[0] - 1 + dy = V.shape[1] - 1 + dx = V.shape[2] - 1 + + if x < 0: + x = dx + x + 1 + elif x > dx: + x = x - dx - 1 + + if y < 0: + y = dy + y + 1 + elif y > dy: + y = y - dy - 1 + + if z < 0: + z = dz + z + 1 + elif z > dz: + z = z - dz - 1 + + return V[z, y, x] + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef void calc_coef_tricub(image_t[:, :, :] V, double x, double y, double z, double [64] coef) nogil: + cdef int xi = floor(x) + cdef int yi = floor(y) + cdef int zi = floor(z) + + cdef double[64] _x + + cdef int i, j + + _x[0] = _G(V, xi, yi, zi) + _x[1] = _G(V, xi + 1, yi, zi) + _x[2] = _G(V, xi, yi + 1, zi) + _x[3] = _G(V, xi + 1, yi + 1, zi) + _x[4] = _G(V, xi, yi, zi + 1) + _x[5] = _G(V, xi + 1, yi, zi + 1) + _x[6] = _G(V, xi, yi + 1, zi + 1) + _x[7] = _G(V, xi + 1, yi + 1, zi + 1) + + _x[8] = 0.5*(_G(V, xi+1,yi,zi) - _G(V, xi-1, yi, zi)) + _x[9] = 0.5*(_G(V, xi+2,yi,zi) - _G(V, xi, yi, zi)) + _x[10] = 0.5*(_G(V, xi+1, yi+1,zi) - _G(V, xi-1, yi+1, zi)) + _x[11] = 0.5*(_G(V, xi+2, yi+1,zi) - _G(V, xi, yi+1, zi)) + _x[12] = 0.5*(_G(V, xi+1, yi,zi+1) - _G(V, xi-1, yi, zi+1)) + _x[13] = 0.5*(_G(V, xi+2, yi,zi+1) - _G(V, xi, yi, zi+1)) + _x[14] = 0.5*(_G(V, xi+1, yi+1,zi+1) - _G(V, xi-1, yi+1, zi+1)) + _x[15] = 0.5*(_G(V, xi+2, yi+1,zi+1) - _G(V, xi, yi+1, zi+1)) + _x[16] = 0.5*(_G(V, xi, yi+1,zi) - _G(V, xi, yi-1, zi)) + _x[17] = 0.5*(_G(V, xi+1, yi+1,zi) - _G(V, xi+1, yi-1, zi)) + _x[18] = 0.5*(_G(V, xi, yi+2,zi) - _G(V, xi, yi, zi)) + _x[19] = 0.5*(_G(V, xi+1, yi+2,zi) - _G(V, xi+1, yi, zi)) + _x[20] = 0.5*(_G(V, xi, yi+1,zi+1) - _G(V, xi, yi-1, zi+1)) + _x[21] = 0.5*(_G(V, xi+1, yi+1,zi+1) - _G(V, xi+1, yi-1, zi+1)) + _x[22] = 0.5*(_G(V, xi, yi+2,zi+1) - _G(V, xi, yi, zi+1)) + _x[23] = 0.5*(_G(V, xi+1, yi+2,zi+1) - _G(V, xi+1, yi, zi+1)) + _x[24] = 0.5*(_G(V, xi, yi,zi+1) - _G(V, xi, yi, zi-1)) + _x[25] = 0.5*(_G(V, xi+1, yi,zi+1) - _G(V, xi+1, yi, zi-1)) + _x[26] = 0.5*(_G(V, xi, yi+1,zi+1) - _G(V, xi, yi+1, zi-1)) + _x[27] = 0.5*(_G(V, xi+1, yi+1,zi+1) - _G(V, xi+1, yi+1, zi-1)) + _x[28] = 0.5*(_G(V, xi, yi,zi+2) - _G(V, xi, yi, zi)) + _x[29] = 0.5*(_G(V, xi+1, yi,zi+2) - _G(V, xi+1, yi, zi)) + _x[30] = 0.5*(_G(V, xi, yi+1,zi+2) - _G(V, xi, yi+1, zi)) + _x[31] = 0.5*(_G(V, xi+1, yi+1,zi+2) - _G(V, xi+1, yi+1, zi)) + + _x [32] = 0.25*(_G(V, xi+1, yi+1, zi) - _G(V, xi-1, yi+1, zi) - _G(V, xi+1, yi-1, zi) + _G(V, xi-1, yi-1, zi)) + _x [33] = 0.25*(_G(V, xi+2, yi+1, zi) - _G(V, xi, yi+1, zi) - _G(V, xi+2, yi-1, zi) + _G(V, xi, yi-1, zi)) + _x [34] = 0.25*(_G(V, xi+1, yi+2, zi) - _G(V, xi-1, yi+2, zi) - _G(V, xi+1, yi, zi) + _G(V, xi-1, yi, zi)) + _x [35] = 0.25*(_G(V, xi+2, yi+2, zi) - _G(V, xi, yi+2, zi) - _G(V, xi+2, yi, zi) + _G(V, xi, yi, zi)) + _x [36] = 0.25*(_G(V, xi+1, yi+1, zi+1) - _G(V, xi-1, yi+1, zi+1) - _G(V, xi+1, yi-1, zi+1) + _G(V, xi-1, yi-1, zi+1)) + _x [37] = 0.25*(_G(V, xi+2, yi+1, zi+1) - _G(V, xi, yi+1, zi+1) - _G(V, xi+2, yi-1, zi+1) + _G(V, xi, yi-1, zi+1)) + _x [38] = 0.25*(_G(V, xi+1, yi+2, zi+1) - _G(V, xi-1, yi+2, zi+1) - _G(V, xi+1, yi, zi+1) + _G(V, xi-1, yi, zi+1)) + _x [39] = 0.25*(_G(V, xi+2, yi+2, zi+1) - _G(V, xi, yi+2, zi+1) - _G(V, xi+2, yi, zi+1) + _G(V, xi, yi, zi+1)) + _x [40] = 0.25*(_G(V, xi+1, yi, zi+1) - _G(V, xi-1, yi, zi+1) - _G(V, xi+1, yi, zi-1) + _G(V, xi-1, yi, zi-1)) + _x [41] = 0.25*(_G(V, xi+2, yi, zi+1) - _G(V, xi, yi, zi+1) - _G(V, xi+2, yi, zi-1) + _G(V, xi, yi, zi-1)) + _x [42] = 0.25*(_G(V, xi+1, yi+1, zi+1) - _G(V, xi-1, yi+1, zi+1) - _G(V, xi+1, yi+1, zi-1) + _G(V, xi-1, yi+1, zi-1)) + _x [43] = 0.25*(_G(V, xi+2, yi+1, zi+1) - _G(V, xi, yi+1, zi+1) - _G(V, xi+2, yi+1, zi-1) + _G(V, xi, yi+1, zi-1)) + _x [44] = 0.25*(_G(V, xi+1, yi, zi+2) - _G(V, xi-1, yi, zi+2) - _G(V, xi+1, yi, zi) + _G(V, xi-1, yi, zi)) + _x [45] = 0.25*(_G(V, xi+2, yi, zi+2) - _G(V, xi, yi, zi+2) - _G(V, xi+2, yi, zi) + _G(V, xi, yi, zi)) + _x [46] = 0.25*(_G(V, xi+1, yi+1, zi+2) - _G(V, xi-1, yi+1, zi+2) - _G(V, xi+1, yi+1, zi) + _G(V, xi-1, yi+1, zi)) + _x [47] = 0.25*(_G(V, xi+2, yi+1, zi+2) - _G(V, xi, yi+1, zi+2) - _G(V, xi+2, yi+1, zi) + _G(V, xi, yi+1, zi)) + _x [48] = 0.25*(_G(V, xi, yi+1, zi+1) - _G(V, xi, yi-1, zi+1) - _G(V, xi, yi+1, zi-1) + _G(V, xi, yi-1, zi-1)) + _x [49] = 0.25*(_G(V, xi+1, yi+1, zi+1) - _G(V, xi+1, yi-1, zi+1) - _G(V, xi+1, yi+1, zi-1) + _G(V, xi+1, yi-1, zi-1)) + _x [50] = 0.25*(_G(V, xi, yi+2, zi+1) - _G(V, xi, yi, zi+1) - _G(V, xi, yi+2, zi-1) + _G(V, xi, yi, zi-1)) + _x [51] = 0.25*(_G(V, xi+1, yi+2, zi+1) - _G(V, xi+1, yi, zi+1) - _G(V, xi+1, yi+2, zi-1) + _G(V, xi+1, yi, zi-1)) + _x [52] = 0.25*(_G(V, xi, yi+1, zi+2) - _G(V, xi, yi-1, zi+2) - _G(V, xi, yi+1, zi) + _G(V, xi, yi-1, zi)) + _x [53] = 0.25*(_G(V, xi+1, yi+1, zi+2) - _G(V, xi+1, yi-1, zi+2) - _G(V, xi+1, yi+1, zi) + _G(V, xi+1, yi-1, zi)) + _x [54] = 0.25*(_G(V, xi, yi+2, zi+2) - _G(V, xi, yi, zi+2) - _G(V, xi, yi+2, zi) + _G(V, xi, yi, zi)) + _x [55] = 0.25*(_G(V, xi+1, yi+2, zi+2) - _G(V, xi+1, yi, zi+2) - _G(V, xi+1, yi+2, zi) + _G(V, xi+1, yi, zi)) + + _x[56] = 0.125*(_G(V, xi+1, yi+1, zi+1) - _G(V, xi-1, yi+1, zi+1) - _G(V, xi+1, yi-1, zi+1) + _G(V, xi-1, yi-1, zi+1) - _G(V, xi+1, yi+1, zi-1) + _G(V, xi-1,yi+1,zi-1)+_G(V, xi+1,yi-1,zi-1)-_G(V, xi-1,yi-1,zi-1)) + _x[57] = 0.125*(_G(V, xi+2, yi+1, zi+1) - _G(V, xi, yi+1, zi+1) - _G(V, xi+2, yi-1, zi+1) + _G(V, xi, yi-1, zi+1) - _G(V, xi+2, yi+1, zi-1) + _G(V, xi,yi+1,zi-1)+_G(V, xi+2,yi-1,zi-1)-_G(V, xi,yi-1,zi-1)) + _x[58] = 0.125*(_G(V, xi+1, yi+2, zi+1) - _G(V, xi-1, yi+2, zi+1) - _G(V, xi+1, yi, zi+1) + _G(V, xi-1, yi, zi+1) - _G(V, xi+1, yi+2, zi-1) + _G(V, xi-1,yi+2,zi-1)+_G(V, xi+1,yi,zi-1)-_G(V, xi-1,yi,zi-1)) + _x[59] = 0.125*(_G(V, xi+2, yi+2, zi+1) - _G(V, xi, yi+2, zi+1) - _G(V, xi+2, yi, zi+1) + _G(V, xi, yi, zi+1) - _G(V, xi+2, yi+2, zi-1) + _G(V, xi,yi+2,zi-1)+_G(V, xi+2,yi,zi-1)-_G(V, xi,yi,zi-1)) + _x[60] = 0.125*(_G(V, xi+1, yi+1, zi+2) - _G(V, xi-1, yi+1, zi+2) - _G(V, xi+1, yi-1, zi+2) + _G(V, xi-1, yi-1, zi+2) - _G(V, xi+1, yi+1, zi) + _G(V, xi-1,yi+1,zi)+_G(V, xi+1,yi-1,zi)-_G(V, xi-1,yi-1,zi)) + _x[61] = 0.125*(_G(V, xi+2, yi+1, zi+2) - _G(V, xi, yi+1, zi+2) - _G(V, xi+2, yi-1, zi+2) + _G(V, xi, yi-1, zi+2) - _G(V, xi+2, yi+1, zi) + _G(V, xi,yi+1,zi)+_G(V, xi+2,yi-1,zi)-_G(V, xi,yi-1,zi)) + _x[62] = 0.125*(_G(V, xi+1, yi+2, zi+2) - _G(V, xi-1, yi+2, zi+2) - _G(V, xi+1, yi, zi+2) + _G(V, xi-1, yi, zi+2) - _G(V, xi+1, yi+2, zi) + _G(V, xi-1,yi+2,zi)+_G(V, xi+1,yi,zi)-_G(V, xi-1,yi,zi)) + _x[63] = 0.125*(_G(V, xi+2, yi+2, zi+2) - _G(V, xi, yi+2, zi+2) - _G(V, xi+2, yi, zi+2) + _G(V, xi, yi, zi+2) - _G(V, xi+2, yi+2, zi) + _G(V, xi,yi+2,zi)+_G(V, xi+2,yi,zi)-_G(V, xi,yi,zi)) + + for j in prange(64): + coef[j] = 0.0 + for i in xrange(64): + coef[j] += (temp[j][i] * _x[i]) + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef double tricub_interpolate(image_t[:, :, :] V, double x, double y, double z) nogil: + # From: Tricubic interpolation in three dimensions. Lekien and Marsden + cdef double[64] coef + cdef double result = 0.0 + calc_coef_tricub(V, x, y, z, coef) + + cdef int i, j, k + + cdef int xi = floor(x) + cdef int yi = floor(y) + cdef int zi = floor(z) + + for i in xrange(4): + for j in xrange(4): + for k in xrange(4): + result += (coef[i+4*j+16*k] * ((x-xi)**i) * ((y-yi)**j) * ((z-zi)**k)) + # return V[z, y, x] + # with gil: + # print result + return result + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef double cubicInterpolate(double p[4], double x) nogil: + return p[1] + 0.5 * x*(p[2] - p[0] + x*(2.0*p[0] - 5.0*p[1] + 4.0*p[2] - p[3] + x*(3.0*(p[1] - p[2]) + p[3] - p[0]))) + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef double bicubicInterpolate (double p[4][4], double x, double y) nogil: + cdef double arr[4] + arr[0] = cubicInterpolate(p[0], y) + arr[1] = cubicInterpolate(p[1], y) + arr[2] = cubicInterpolate(p[2], y) + arr[3] = cubicInterpolate(p[3], y) + return cubicInterpolate(arr, x) + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef double tricubicInterpolate(image_t[:, :, :] V, double x, double y, double z) nogil: + # From http://www.paulinternet.nl/?page=bicubic + cdef double p[4][4][4] + + cdef int xi = floor(x) + cdef int yi = floor(y) + cdef int zi = floor(z) + + cdef int i, j, k + + for i in xrange(4): + for j in xrange(4): + for k in xrange(4): + p[i][j][k] = _G(V, xi + i -1, yi + j -1, zi + k - 1) + + cdef double arr[4] + arr[0] = bicubicInterpolate(p[0], y-yi, z-zi) + arr[1] = bicubicInterpolate(p[1], y-yi, z-zi) + arr[2] = bicubicInterpolate(p[2], y-yi, z-zi) + arr[3] = bicubicInterpolate(p[3], y-yi, z-zi) + return cubicInterpolate(arr, x-xi) + + +def tricub_interpolate_py(image_t[:, :, :] V, double x, double y, double z): + return tricub_interpolate(V, x, y, z) + +def tricub_interpolate2_py(image_t[:, :, :] V, double x, double y, double z): + return tricubicInterpolate(V, x, y, z) + +def trilin_interpolate_py(image_t[:, :, :] V, double x, double y, double z): + return interpolate(V, x, y, z) diff --git a/invesalius_cy/mips.pyx b/invesalius_cy/mips.pyx new file mode 100644 index 0000000..6b0e0a2 --- /dev/null +++ b/invesalius_cy/mips.pyx @@ -0,0 +1,379 @@ +#http://en.wikipedia.org/wiki/Local_maximum_intensity_projection +import numpy as np +cimport numpy as np +cimport cython + +from libc.math cimport floor, ceil, sqrt, fabs +from cython.parallel import prange + +DTYPE = np.uint8 +ctypedef np.uint8_t DTYPE_t + +DTYPE16 = np.int16 +ctypedef np.int16_t DTYPE16_t + +DTYPEF32 = np.float32 +ctypedef np.float32_t DTYPEF32_t + +@cython.boundscheck(False) # turn of bounds-checking for entire function +def lmip(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t tmin, + DTYPE16_t tmax, np.ndarray[DTYPE16_t, ndim=2] out): + cdef DTYPE16_t max + cdef int start + cdef int sz = image.shape[0] + cdef int sy = image.shape[1] + cdef int sx = image.shape[2] + + # AXIAL + if axis == 0: + for x in xrange(sx): + for y in xrange(sy): + max = image[0, y, x] + if max >= tmin and max <= tmax: + start = 1 + else: + start = 0 + for z in xrange(sz): + if image[z, y, x] > max: + max = image[z, y, x] + + elif image[z, y, x] < max and start: + break + + if image[z, y, x] >= tmin and image[z, y, x] <= tmax: + start = 1 + + out[y, x] = max + + #CORONAL + elif axis == 1: + for z in xrange(sz): + for x in xrange(sx): + max = image[z, 0, x] + if max >= tmin and max <= tmax: + start = 1 + else: + start = 0 + for y in xrange(sy): + if image[z, y, x] > max: + max = image[z, y, x] + + elif image[z, y, x] < max and start: + break + + if image[z, y, x] >= tmin and image[z, y, x] <= tmax: + start = 1 + + out[z, x] = max + + #CORONAL + elif axis == 2: + for z in xrange(sz): + for y in xrange(sy): + max = image[z, y, 0] + if max >= tmin and max <= tmax: + start = 1 + else: + start = 0 + for x in xrange(sx): + if image[z, y, x] > max: + max = image[z, y, x] + + elif image[z, y, x] < max and start: + break + + if image[z, y, x] >= tmin and image[z, y, x] <= tmax: + start = 1 + + out[z, y] = max + + +cdef DTYPE16_t get_colour(DTYPE16_t vl, DTYPE16_t wl, DTYPE16_t ww): + cdef DTYPE16_t out_colour + cdef DTYPE16_t min_value = wl - (ww / 2) + cdef DTYPE16_t max_value = wl + (ww / 2) + if vl < min_value: + out_colour = min_value + elif vl > max_value: + out_colour = max_value + else: + out_colour = vl + + return out_colour + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +cdef float get_opacity(DTYPE16_t vl, DTYPE16_t wl, DTYPE16_t ww) nogil: + cdef float out_opacity + cdef DTYPE16_t min_value = wl - (ww / 2) + cdef DTYPE16_t max_value = wl + (ww / 2) + if vl < min_value: + out_opacity = 0.0 + elif vl > max_value: + out_opacity = 1.0 + else: + out_opacity = 1.0/(max_value - min_value) * (vl - min_value) + + return out_opacity + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +cdef float get_opacity_f32(DTYPEF32_t vl, DTYPE16_t wl, DTYPE16_t ww) nogil: + cdef float out_opacity + cdef DTYPE16_t min_value = wl - (ww / 2) + cdef DTYPE16_t max_value = wl + (ww / 2) + if vl < min_value: + out_opacity = 0.0 + elif vl > max_value: + out_opacity = 1.0 + else: + out_opacity = 1.0/(max_value - min_value) * (vl - min_value) + + return out_opacity + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +def mida(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t wl, + DTYPE16_t ww, np.ndarray[DTYPE16_t, ndim=2] out): + cdef int sz = image.shape[0] + cdef int sy = image.shape[1] + cdef int sx = image.shape[2] + + cdef DTYPE16_t min = image.min() + cdef DTYPE16_t max = image.max() + cdef DTYPE16_t vl + + cdef DTYPE16_t min_value = wl - (ww / 2) + cdef DTYPE16_t max_value = wl + (ww / 2) + + cdef float fmax=0.0 + cdef float fpi + cdef float dl + cdef float bt + + cdef float alpha + cdef float alpha_p = 0.0 + cdef float colour + cdef float colour_p = 0 + + cdef int x, y, z + + # AXIAL + if axis == 0: + for x in prange(sx, nogil=True): + for y in xrange(sy): + fmax = 0.0 + alpha_p = 0.0 + colour_p = 0.0 + for z in xrange(sz): + vl = image[z, y, x] + fpi = 1.0/(max - min) * (vl - min) + if fpi > fmax: + dl = fpi - fmax + fmax = fpi + else: + dl = 0.0 + + bt = 1.0 - dl + + colour = fpi + alpha = get_opacity(vl, wl, ww) + colour = (bt * colour_p) + (1 - bt * alpha_p) * colour * alpha + alpha = (bt * alpha_p) + (1 - bt * alpha_p) * alpha + + colour_p = colour + alpha_p = alpha + + if alpha >= 1.0: + break + + + #out[y, x] = ((max_value - min_value) * colour + min_value) + out[y, x] = ((max - min) * colour + min) + + + #CORONAL + elif axis == 1: + for z in prange(sz, nogil=True): + for x in xrange(sx): + fmax = 0.0 + alpha_p = 0.0 + colour_p = 0.0 + for y in xrange(sy): + vl = image[z, y, x] + fpi = 1.0/(max - min) * (vl - min) + if fpi > fmax: + dl = fpi - fmax + fmax = fpi + else: + dl = 0.0 + + bt = 1.0 - dl + + colour = fpi + alpha = get_opacity(vl, wl, ww) + colour = (bt * colour_p) + (1 - bt * alpha_p) * colour * alpha + alpha = (bt * alpha_p) + (1 - bt * alpha_p) * alpha + + colour_p = colour + alpha_p = alpha + + if alpha >= 1.0: + break + + out[z, x] = ((max - min) * colour + min) + + #AXIAL + elif axis == 2: + for z in prange(sz, nogil=True): + for y in xrange(sy): + fmax = 0.0 + alpha_p = 0.0 + colour_p = 0.0 + for x in xrange(sx): + vl = image[z, y, x] + fpi = 1.0/(max - min) * (vl - min) + if fpi > fmax: + dl = fpi - fmax + fmax = fpi + else: + dl = 0.0 + + bt = 1.0 - dl + + colour = fpi + alpha = get_opacity(vl, wl, ww) + colour = (bt * colour_p) + (1 - bt * alpha_p) * colour * alpha + alpha = (bt * alpha_p) + (1 - bt * alpha_p) * alpha + + colour_p = colour + alpha_p = alpha + + if alpha >= 1.0: + break + + out[z, y] = ((max - min) * colour + min) + + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +cdef inline void finite_difference(DTYPE16_t[:, :, :] image, + int x, int y, int z, float h, float *g) nogil: + cdef int px, py, pz, fx, fy, fz + + cdef int sz = image.shape[0] + cdef int sy = image.shape[1] + cdef int sx = image.shape[2] + + cdef float gx, gy, gz + + if x == 0: + px = 0 + fx = 1 + elif x == sx - 1: + px = x - 1 + fx = x + else: + px = x - 1 + fx = x + 1 + + if y == 0: + py = 0 + fy = 1 + elif y == sy - 1: + py = y - 1 + fy = y + else: + py = y - 1 + fy = y + 1 + + if z == 0: + pz = 0 + fz = 1 + elif z == sz - 1: + pz = z - 1 + fz = z + else: + pz = z - 1 + fz = z + 1 + + gx = (image[z, y, fx] - image[z, y, px]) / (2*h) + gy = (image[z, fy, x] - image[z, py, x]) / (2*h) + gz = (image[fz, y, x] - image[pz, y, x]) / (2*h) + + g[0] = gx + g[1] = gy + g[2] = gz + + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +cdef inline float calc_fcm_itensity(DTYPE16_t[:, :, :] image, + int x, int y, int z, float n, float* dir) nogil: + cdef float g[3] + finite_difference(image, x, y, z, 1.0, g) + cdef float gm = sqrt(g[0]*g[0] + g[1]*g[1] + g[2]*g[2]) + cdef float d = g[0]*dir[0] + g[1]*dir[1] + g[2]*dir[2] + cdef float sf = (1.0 - fabs(d/gm))**n + #alpha = get_opacity_f32(gm, wl, ww) + cdef float vl = gm * sf + return vl + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +def fast_countour_mip(np.ndarray[DTYPE16_t, ndim=3] image, + float n, + int axis, + DTYPE16_t wl, DTYPE16_t ww, + int tmip, + np.ndarray[DTYPE16_t, ndim=2] out): + cdef int sz = image.shape[0] + cdef int sy = image.shape[1] + cdef int sx = image.shape[2] + cdef float gm + cdef float alpha + cdef float sf + cdef float d + + cdef float* g + cdef float* dir = [ 0, 0, 0 ] + + cdef DTYPE16_t[:, :, :] vimage = image + cdef np.ndarray[DTYPE16_t, ndim=3] tmp = np.empty_like(image) + + cdef DTYPE16_t min = image.min() + cdef DTYPE16_t max = image.max() + cdef float fmin = min + cdef float fmax = max + cdef float vl + cdef DTYPE16_t V + + cdef int x, y, z + + if axis == 0: + dir[2] = 1.0 + elif axis == 1: + dir[1] = 1.0 + elif axis == 2: + dir[0] = 1.0 + + for z in prange(sz, nogil=True): + for y in range(sy): + for x in range(sx): + vl = calc_fcm_itensity(vimage, x, y, z, n, dir) + tmp[z, y, x] = vl + + cdef DTYPE16_t tmin = tmp.min() + cdef DTYPE16_t tmax = tmp.max() + + #tmp = ((max - min)/(tmax - tmin)) * (tmp - tmin) + min + + if tmip == 0: + out[:] = tmp.max(axis) + elif tmip == 1: + lmip(tmp, axis, 700, 3033, out) + elif tmip == 2: + mida(tmp, axis, wl, ww, out) diff --git a/invesalius_cy/transforms.pyx b/invesalius_cy/transforms.pyx new file mode 100644 index 0000000..5c746d2 --- /dev/null +++ b/invesalius_cy/transforms.pyx @@ -0,0 +1,174 @@ +import numpy as np +cimport numpy as np +cimport cython + +from .cy_my_types cimport image_t +from .interpolation cimport interpolate, tricub_interpolate, tricubicInterpolate, lanczos3, nearest_neighbour_interp + +from libc.math cimport floor, ceil, sqrt, fabs, round +from cython.parallel import prange + +ctypedef double (*interp_function)(image_t[:, :, :], double, double, double) nogil + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef void mul_mat4_vec4(double[:, :] M, + double* coord, + double* out) nogil: + + out[0] = coord[0] * M[0, 0] + coord[1] * M[0, 1] + coord[2] * M[0, 2] + coord[3] * M[0, 3] + out[1] = coord[0] * M[1, 0] + coord[1] * M[1, 1] + coord[2] * M[1, 2] + coord[3] * M[1, 3] + out[2] = coord[0] * M[2, 0] + coord[1] * M[2, 1] + coord[2] * M[2, 2] + coord[3] * M[2, 3] + out[3] = coord[0] * M[3, 0] + coord[1] * M[3, 1] + coord[2] * M[3, 2] + coord[3] * M[3, 3] + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef image_t coord_transform(image_t[:, :, :] volume, double[:, :] M, int x, int y, int z, double sx, double sy, double sz, interp_function f_interp, image_t cval) nogil: + + cdef double coord[4] + coord[0] = z*sz + coord[1] = y*sy + coord[2] = x*sx + coord[3] = 1.0 + + cdef double _ncoord[4] + _ncoord[3] = 1 + # _ncoord[:] = [0.0, 0.0, 0.0, 1.0] + + cdef unsigned int dz, dy, dx + dz = volume.shape[0] + dy = volume.shape[1] + dx = volume.shape[2] + + + mul_mat4_vec4(M, coord, _ncoord) + + cdef double nz, ny, nx + nz = (_ncoord[0]/_ncoord[3])/sz + ny = (_ncoord[1]/_ncoord[3])/sy + nx = (_ncoord[2]/_ncoord[3])/sx + + cdef double v + + if 0 <= nz <= (dz-1) and 0 <= ny <= (dy-1) and 0 <= nx <= (dx-1): + return f_interp(volume, nx, ny, nz) + # if minterpol == 0: + # return nearest_neighbour_interp(volume, nx, ny, nz) + # elif minterpol == 1: + # return interpolate(volume, nx, ny, nz) + # elif minterpol == 2: + # v = tricubicInterpolate(volume, nx, ny, nz) + # if (v < cval): + # v = cval + # return v + # else: + # v = lanczos3(volume, nx, ny, nz) + # if (v < cval): + # v = cval + # return v + else: + return cval + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +def apply_view_matrix_transform(image_t[:, :, :] volume, + spacing, + double[:, :] M, + unsigned int n, str orientation, + int minterpol, + image_t cval, + image_t[:, :, :] out): + + cdef int dz, dy, dx + cdef int z, y, x + dz = volume.shape[0] + dy = volume.shape[1] + dx = volume.shape[2] + + cdef unsigned int odz, ody, odx + odz = out.shape[0] + ody = out.shape[1] + odx = out.shape[2] + + cdef unsigned int count = 0 + + cdef double sx, sy, sz + sx = spacing[0] + sy = spacing[1] + sz = spacing[2] + + cdef interp_function f_interp + + if minterpol == 0: + f_interp = nearest_neighbour_interp + elif minterpol == 1: + f_interp = interpolate + elif minterpol == 2: + f_interp = tricubicInterpolate + else: + f_interp = lanczos3 + + if orientation == 'AXIAL': + for z in xrange(n, n+odz): + for y in prange(dy, nogil=True): + for x in xrange(dx): + out[count, y, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval) + count += 1 + + elif orientation == 'CORONAL': + for y in xrange(n, n+ody): + for z in prange(dz, nogil=True): + for x in xrange(dx): + out[z, count, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval) + count += 1 + + elif orientation == 'SAGITAL': + for x in xrange(n, n+odx): + for z in prange(dz, nogil=True): + for y in xrange(dy): + out[z, y, count] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval) + count += 1 + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +def convolve_non_zero(image_t[:, :, :] volume, + image_t[:, :, :] kernel, + image_t cval): + cdef Py_ssize_t x, y, z, sx, sy, sz, kx, ky, kz, skx, sky, skz, i, j, k + cdef image_t v + + cdef image_t[:, :, :] out = np.zeros_like(volume) + + sz = volume.shape[0] + sy = volume.shape[1] + sx = volume.shape[2] + + skz = kernel.shape[0] + sky = kernel.shape[1] + skx = kernel.shape[2] + + for z in prange(sz, nogil=True): + for y in xrange(sy): + for x in xrange(sx): + if volume[z, y, x] != 0: + for k in xrange(skz): + kz = z - skz // 2 + k + for j in xrange(sky): + ky = y - sky // 2 + j + for i in xrange(skx): + kx = x - skx // 2 + i + + if 0 <= kz < sz and 0 <= ky < sy and 0 <= kx < sx: + v = volume[kz, ky, kx] + else: + v = cval + + out[z, y, x] += (v * kernel[k, j, i]) + return np.asarray(out) diff --git a/setup.py b/setup.py index 514815b..3187fa6 100644 --- a/setup.py +++ b/setup.py @@ -40,25 +40,25 @@ setup( ext_modules=cythonize( [ Extension( - "invesalius.data.mips", - ["invesalius/data/mips.pyx"], + "invesalius_cy.mips", + ["invesalius_cy/mips.pyx"], ), Extension( - "invesalius.data.interpolation", - ["invesalius/data/interpolation.pyx"], + "invesalius_cy.interpolation", + ["invesalius_cy/interpolation.pyx"], ), Extension( - "invesalius.data.transforms", - ["invesalius/data/transforms.pyx"], + "invesalius_cy.transforms", + ["invesalius_cy/transforms.pyx"], ), Extension( - "invesalius.data.floodfill", - ["invesalius/data/floodfill.pyx"], + "invesalius_cy.floodfill", + ["invesalius_cy/floodfill.pyx"], language="c++", ), Extension( - "invesalius.data.cy_mesh", - ["invesalius/data/cy_mesh.pyx"], + "invesalius_cy.cy_mesh", + ["invesalius_cy/cy_mesh.pyx"], language="c++", ), ] -- libgit2 0.21.2