diff --git a/invesalius_cy_/__init__.py b/invesalius_cy_/__init__.py deleted file mode 100644 index e69de29bb..000000000 diff --git a/invesalius_cy_/cy_mesh.pyx b/invesalius_cy_/cy_mesh.pyx deleted file mode 100644 index 00e0fdd6b..000000000 --- a/invesalius_cy_/cy_mesh.pyx +++ /dev/null @@ -1,478 +0,0 @@ -# distutils: language = c++ -# distutils: define_macros=NPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION -# cython: boundscheck=False -# cython: wraparound=False -# cython: initializedcheck=False -# cython: cdivision=True -# cython: nonecheck=False -# cython: language_level=3 - -import os -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 cimport prange -cimport openmp - -from .cy_my_types cimport vertex_t, normal_t, vertex_id_t - -import numpy as np - -from vtkmodules.util import numpy_support -from vtkmodules.vtkCommonCore import vtkPoints -from vtkmodules.vtkCommonDataModel import vtkCellArray, vtkPolyData - -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 range(_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 = vtkPoints() - points.SetData(numpy_support.numpy_to_vtk(vertices)) - - id_triangles = numpy_support.numpy_to_vtkIdTypeArray(faces) - triangles = vtkCellArray() - triangles.SetCells(faces.shape[0], id_triangles) - - pd = vtkPolyData() - pd.SetPoints(points) - pd.SetPolys(triangles) - - return pd - - cdef vector[vertex_id_t]* get_faces_by_vertex(self, int v_id) noexcept 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) noexcept 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) noexcept 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) noexcept 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) - dmax = dmax * dmax - 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 range(nf): - f_id = deref(idfaces)[nid] - if status_f.find(f_id) == status_f.end(): - status_f[f_id] = True - - for j in range(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 = (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) noexcept 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) - - cdef openmp.omp_lock_t lock - openmp.omp_init_lock(&lock) - - 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 range(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]: - openmp.omp_set_lock(&lock) - deref(weights)[vj_id] = value - openmp.omp_unset_lock(&lock) - - del near_vertices - - # for i in range(msize): - # if mesh.is_border(i): - # deref(weights)[i] = 0.0 - - # cdef vertex_id_t v0, v1, v2 - # for i in range(mesh.faces.shape[0]): - # for j in range(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 - - openmp.omp_destroy_lock(&lock) - return weights - - -cdef inline Point calc_d(Mesh mesh, vertex_id_t v_id) noexcept 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) noexcept 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 range(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 range(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) nogil: - """ - 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 range(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 deleted file mode 100644 index 13d6c6cf7..000000000 --- a/invesalius_cy_/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_cy_/floodfill.pyx b/invesalius_cy_/floodfill.pyx deleted file mode 100644 index 031b137e6..000000000 --- a/invesalius_cy_/floodfill.pyx +++ /dev/null @@ -1,264 +0,0 @@ -# distutils: define_macros=NPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION -# cython: boundscheck=False -# cython: wraparound=False -# cython: initializedcheck=False -# cython: cdivision=True -# cython: nonecheck=False -# cython: language_level=3 - -import numpy as np -cimport numpy as np -cimport cython - -from collections import deque - -from cython.parallel cimport 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 - - -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) - - -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) - - -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 - - -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 range(odz): - zo = z + k - offset_z - for j in range(ody): - yo = y + j - offset_y - for i in range(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 - - -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 - - -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 range(dz): - for y in range(dy): - for x in range(dx): - sizes[labels[z, y, x]] += 1 - - #Checking if any hole will be filled - for i in range(nlabels + 1): - if sizes[i] <= max_size: - modified = True - - if not modified: - return 0 - - for z in prange(dz, nogil=True): - for y in range(dy): - for x in range(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 deleted file mode 100644 index b63cc4155..000000000 --- a/invesalius_cy_/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_cy_/interpolation.pyx b/invesalius_cy_/interpolation.pyx deleted file mode 100644 index 4301ab489..000000000 --- a/invesalius_cy_/interpolation.pyx +++ /dev/null @@ -1,373 +0,0 @@ -# distutils: define_macros=NPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION -# cython: boundscheck=False -# cython: wraparound=False -# cython: initializedcheck=False -# cython: cdivision=True -# cython: nonecheck=False -# cython: language_level=3 - -# 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 cimport 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] -] - - -cdef inline image_t _G(const image_t[:, :, :] V, int x, int y, int z) noexcept 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] - - -cdef double nearest_neighbour_interp(const image_t[:, :, :] V, double x, double y, double z) noexcept nogil: - return V[(z), (y), (x)] - -cdef double interpolate(const image_t[:, :, :] V, double x, double y, double z) noexcept 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 - - -cdef inline double lanczos3_L(double x, int a) noexcept 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 - - -cdef double lanczos3(const image_t[:, :, :] V, double x, double y, double z) noexcept 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 range(zi, zf): - n = 0 - for j in range(yi, yf): - lx = 0 - for i in range(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 range(zi, zf): - n = 0 - ly = 0 - for j in range(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 range(zi, zf): - lz += temp_y[m] * lanczos3_L(z - k, a) - m += 1 - - return lz - - - - -cdef void calc_coef_tricub(image_t[:, :, :] V, double x, double y, double z, double [64] coef) noexcept 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, nogil=True): - coef[j] = 0.0 - for i in range(64): - coef[j] += (temp[j][i] * _x[i]) - - -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 range(4): - for j in range(4): - for k in range(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 - - -cdef double cubicInterpolate(double p[4], double x) noexcept 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]))) - - -cdef double bicubicInterpolate (double p[4][4], double x, double y) noexcept 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) - - -cdef double tricubicInterpolate(image_t[:, :, :] V, double x, double y, double z) noexcept 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 range(4): - for j in range(4): - for k in range(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 deleted file mode 100644 index 5055b91db..000000000 --- a/invesalius_cy_/mips.pyx +++ /dev/null @@ -1,374 +0,0 @@ -# distutils: define_macros=NPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION -# cython: boundscheck=False -# cython: wraparound=False -# cython: initializedcheck=False -# cython: cdivision=True -# cython: nonecheck=False -# cython: language_level=3 - -#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 cimport 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 - -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 range(sx): - for y in range(sy): - max = image[0, y, x] - if max >= tmin and max <= tmax: - start = 1 - else: - start = 0 - for z in range(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 range(sz): - for x in range(sx): - max = image[z, 0, x] - if max >= tmin and max <= tmax: - start = 1 - else: - start = 0 - for y in range(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 range(sz): - for y in range(sy): - max = image[z, y, 0] - if max >= tmin and max <= tmax: - start = 1 - else: - start = 0 - for x in range(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 - -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 - -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 - - -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 range(sy): - fmax = 0.0 - alpha_p = 0.0 - colour_p = 0.0 - for z in range(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 range(sx): - fmax = 0.0 - alpha_p = 0.0 - colour_p = 0.0 - for y in range(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 range(sy): - fmax = 0.0 - alpha_p = 0.0 - colour_p = 0.0 - for x in range(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) - - - -cdef inline void finite_difference(DTYPE16_t[:, :, :] image, - int x, int y, int z, float h, float *g) noexcept 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 - - - -cdef inline float calc_fcm_itensity(DTYPE16_t[:, :, :] image, - int x, int y, int z, float n, float* dir) noexcept 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 - -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 deleted file mode 100644 index c44462586..000000000 --- a/invesalius_cy_/transforms.pyx +++ /dev/null @@ -1,171 +0,0 @@ -# distutils: define_macros=NPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION -# cython: boundscheck=False -# cython: wraparound=False -# cython: initializedcheck=False -# cython: cdivision=True -# cython: nonecheck=False -# cython: language_level=3 - -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 cimport prange - -ctypedef double (*interp_function)(image_t[:, :, :], double, double, double) nogil - -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] - - -cdef image_t coord_transform(image_t[:, :, :] volume, double[:, :] M, int x, int y, int z, double sx, double sy, double sz, int minterpol, 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): - 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 - - -def apply_view_matrix_transform(image_t[:, :, :] volume, - tuple 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 double kkkkk = sx * sy * sz - - 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 range(n, n+odz): - for y in prange(dy, nogil=True): - for x in range(dx): - out[count, y, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, minterpol, cval) - count += 1 - - elif orientation == 'CORONAL': - for y in range(n, n+ody): - for z in prange(dz, nogil=True): - for x in range(dx): - out[z, count, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, minterpol, cval) - count += 1 - - elif orientation == 'SAGITAL': - for x in range(n, n+odx): - for z in prange(dz, nogil=True): - for y in range(dy): - out[z, y, count] = coord_transform(volume, M, x, y, z, sx, sy, sz, minterpol, cval) - count += 1 - - -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 range(sy): - for x in range(sx): - if volume[z, y, x] != 0: - for k in range(skz): - kz = z - skz // 2 + k - for j in range(sky): - ky = y - sky // 2 + j - for i in range(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)