Download as
Browse Code »

Commit 18c2bcddd4fbb36ebbef13006c804f4ce661476a

Authored by Thiago Franco de Moraes
1 parent 73ddbb08
Exists in master

Moved all cython codes to invesalius_cy

Showing 20 changed files with 1731 additions and 1731 deletions   Show diff stats
invesalius/data/cy_mesh.pyx
View file @73ddbb0
@@ -1,466 +0,0 @@ @@ -1,466 +0,0 @@
1 -# distutils: language = c++  
2 -#cython: boundscheck=False  
3 -#cython: wraparound=False  
4 -#cython: initializedcheck=False  
5 -#cython: cdivision=True  
6 -#cython: nonecheck=False  
7 -  
8 -import sys  
9 -import time  
10 -cimport numpy as np  
11 -  
12 -from libc.math cimport sin, cos, acos, exp, sqrt, fabs, M_PI  
13 -from libc.stdlib cimport abs as cabs  
14 -from cython.operator cimport dereference as deref, preincrement as inc  
15 -from libcpp.map cimport map  
16 -from libcpp.unordered_map cimport unordered_map  
17 -from libcpp.set cimport set  
18 -from libcpp.vector cimport vector  
19 -from libcpp.pair cimport pair  
20 -from libcpp cimport bool  
21 -from libcpp.deque cimport deque as cdeque  
22 -from cython.parallel import prange  
23 -  
24 -from cy_my_types cimport vertex_t, normal_t, vertex_id_t  
25 -  
26 -import numpy as np  
27 -import vtk  
28 -  
29 -from vtk.util import numpy_support  
30 -  
31 -ctypedef float weight_t  
32 -  
33 -cdef struct Point:  
34 - vertex_t x  
35 - vertex_t y  
36 - vertex_t z  
37 -  
38 -ctypedef pair[vertex_id_t, vertex_id_t] key  
39 -  
40 -  
41 -cdef class Mesh:  
42 - cdef vertex_t[:, :] vertices  
43 - cdef vertex_id_t[:, :] faces  
44 - cdef normal_t[:, :] normals  
45 -  
46 - cdef unordered_map[int, vector[vertex_id_t]] map_vface  
47 - cdef unordered_map[vertex_id_t, int] border_vertices  
48 -  
49 - cdef bool _initialized  
50 -  
51 - def __cinit__(self, pd=None, other=None):  
52 - cdef int i  
53 - cdef map[key, int] edge_nfaces  
54 - cdef map[key, int].iterator it  
55 - if pd:  
56 - self._initialized = True  
57 - _vertices = numpy_support.vtk_to_numpy(pd.GetPoints().GetData())  
58 - _vertices.shape = -1, 3  
59 -  
60 - _faces = numpy_support.vtk_to_numpy(pd.GetPolys().GetData())  
61 - _faces.shape = -1, 4  
62 -  
63 - _normals = numpy_support.vtk_to_numpy(pd.GetCellData().GetArray("Normals"))  
64 - _normals.shape = -1, 3  
65 -  
66 - self.vertices = _vertices  
67 - self.faces = _faces  
68 - self.normals = _normals  
69 -  
70 - for i in xrange(_faces.shape[0]):  
71 - self.map_vface[self.faces[i, 1]].push_back(i)  
72 - self.map_vface[self.faces[i, 2]].push_back(i)  
73 - self.map_vface[self.faces[i, 3]].push_back(i)  
74 -  
75 - edge_nfaces[key(min(self.faces[i, 1], self.faces[i, 2]), max(self.faces[i, 1], self.faces[i, 2]))] += 1  
76 - edge_nfaces[key(min(self.faces[i, 2], self.faces[i, 3]), max(self.faces[i, 2], self.faces[i, 3]))] += 1  
77 - edge_nfaces[key(min(self.faces[i, 1], self.faces[i, 3]), max(self.faces[i, 1], self.faces[i, 3]))] += 1  
78 -  
79 - it = edge_nfaces.begin()  
80 -  
81 - while it != edge_nfaces.end():  
82 - if deref(it).second == 1:  
83 - self.border_vertices[deref(it).first.first] = 1  
84 - self.border_vertices[deref(it).first.second] = 1  
85 -  
86 - inc(it)  
87 -  
88 - elif other:  
89 - _other = <Mesh>other  
90 - self._initialized = True  
91 - self.vertices = _other.vertices.copy()  
92 - self.faces = _other.faces.copy()  
93 - self.normals = _other.normals.copy()  
94 - self.map_vface = unordered_map[int, vector[vertex_id_t]](_other.map_vface)  
95 - self.border_vertices = unordered_map[vertex_id_t, int](_other.border_vertices)  
96 - else:  
97 - self._initialized = False  
98 -  
99 - cdef void copy_to(self, Mesh other):  
100 - """  
101 - Copies self content to other.  
102 - """  
103 - if self._initialized:  
104 - other.vertices[:] = self.vertices  
105 - other.faces[:] = self.faces  
106 - other.normals[:] = self.normals  
107 - other.map_vface = unordered_map[int, vector[vertex_id_t]](self.map_vface)  
108 - other.border_vertices = unordered_map[vertex_id_t, int](self.border_vertices)  
109 - else:  
110 - other.vertices = self.vertices.copy()  
111 - other.faces = self.faces.copy()  
112 - other.normals = self.normals.copy()  
113 -  
114 - other.map_vface = self.map_vface  
115 - other.border_vertices = self.border_vertices  
116 -  
117 - def to_vtk(self):  
118 - """  
119 - Converts Mesh to vtkPolyData.  
120 - """  
121 - vertices = np.asarray(self.vertices)  
122 - faces = np.asarray(self.faces)  
123 - normals = np.asarray(self.normals)  
124 -  
125 - points = vtk.vtkPoints()  
126 - points.SetData(numpy_support.numpy_to_vtk(vertices))  
127 -  
128 - id_triangles = numpy_support.numpy_to_vtkIdTypeArray(faces)  
129 - triangles = vtk.vtkCellArray()  
130 - triangles.SetCells(faces.shape[0], id_triangles)  
131 -  
132 - pd = vtk.vtkPolyData()  
133 - pd.SetPoints(points)  
134 - pd.SetPolys(triangles)  
135 -  
136 - return pd  
137 -  
138 - cdef vector[vertex_id_t]* get_faces_by_vertex(self, int v_id) nogil:  
139 - """  
140 - Returns the faces whose vertex `v_id' is part.  
141 - """  
142 - return &self.map_vface[v_id]  
143 -  
144 - cdef set[vertex_id_t]* get_ring1(self, vertex_id_t v_id) nogil:  
145 - """  
146 - Returns the ring1 of vertex `v_id'  
147 - """  
148 - cdef vertex_id_t f_id  
149 - cdef set[vertex_id_t]* ring1 = new set[vertex_id_t]()  
150 - cdef vector[vertex_id_t].iterator it = self.map_vface[v_id].begin()  
151 -  
152 - while it != self.map_vface[v_id].end():  
153 - f_id = deref(it)  
154 - inc(it)  
155 - if self.faces[f_id, 1] != v_id:  
156 - ring1.insert(self.faces[f_id, 1])  
157 - if self.faces[f_id, 2] != v_id:  
158 - ring1.insert(self.faces[f_id, 2])  
159 - if self.faces[f_id, 3] != v_id:  
160 - ring1.insert(self.faces[f_id, 3])  
161 -  
162 - return ring1  
163 -  
164 - cdef bool is_border(self, vertex_id_t v_id) nogil:  
165 - """  
166 - Check if vertex `v_id' is a vertex border.  
167 - """  
168 - return self.border_vertices.find(v_id) != self.border_vertices.end()  
169 -  
170 - cdef vector[vertex_id_t]* get_near_vertices_to_v(self, vertex_id_t v_id, float dmax) nogil:  
171 - """  
172 - Returns all vertices with distance at most `d' to the vertex `v_id'  
173 -  
174 - Params:  
175 - v_id: id of the vertex  
176 - dmax: the maximum distance.  
177 - """  
178 - cdef vector[vertex_id_t]* idfaces  
179 - cdef vector[vertex_id_t]* near_vertices = new vector[vertex_id_t]()  
180 -  
181 - cdef cdeque[vertex_id_t] to_visit  
182 - cdef unordered_map[vertex_id_t, bool] status_v  
183 - cdef unordered_map[vertex_id_t, bool] status_f  
184 -  
185 - cdef vertex_t *vip  
186 - cdef vertex_t *vjp  
187 -  
188 - cdef float distance  
189 - cdef int nf, nid, j  
190 - cdef vertex_id_t f_id, vj  
191 -  
192 - vip = &self.vertices[v_id, 0]  
193 - to_visit.push_back(v_id)  
194 - while(not to_visit.empty()):  
195 - v_id = to_visit.front()  
196 - to_visit.pop_front()  
197 -  
198 - status_v[v_id] = True  
199 -  
200 - idfaces = self.get_faces_by_vertex(v_id)  
201 - nf = idfaces.size()  
202 -  
203 - for nid in xrange(nf):  
204 - f_id = deref(idfaces)[nid]  
205 - if status_f.find(f_id) == status_f.end():  
206 - status_f[f_id] = True  
207 -  
208 - for j in xrange(3):  
209 - vj = self.faces[f_id, j+1]  
210 - if status_v.find(vj) == status_v.end():  
211 - status_v[vj] = True  
212 - vjp = &self.vertices[vj, 0]  
213 - distance = sqrt((vip[0] - vjp[0]) * (vip[0] - vjp[0]) \  
214 - + (vip[1] - vjp[1]) * (vip[1] - vjp[1]) \  
215 - + (vip[2] - vjp[2]) * (vip[2] - vjp[2]))  
216 - if distance <= dmax:  
217 - near_vertices.push_back(vj)  
218 - to_visit.push_back(vj)  
219 -  
220 - return near_vertices  
221 -  
222 -  
223 -cdef vector[weight_t]* calc_artifacts_weight(Mesh mesh, vector[vertex_id_t]& vertices_staircase, float tmax, float bmin) nogil:  
224 - """  
225 - Calculate the artifact weight based on distance of each vertex to its  
226 - nearest staircase artifact vertex.  
227 -  
228 - Params:  
229 - mesh: Mesh  
230 - vertices_staircase: the identified staircase artifact vertices  
231 - tmax: max distance the vertex must be to its nearest artifact vertex  
232 - to considered to calculate the weight  
233 - bmin: The minimum weight.  
234 - """  
235 - cdef int vi_id, vj_id, nnv, n_ids, i, j  
236 - cdef vector[vertex_id_t]* near_vertices  
237 - cdef weight_t value  
238 - cdef float d  
239 - n_ids = vertices_staircase.size()  
240 -  
241 - cdef vertex_t* vi  
242 - cdef vertex_t* vj  
243 - cdef size_t msize  
244 -  
245 - msize = mesh.vertices.shape[0]  
246 - cdef vector[weight_t]* weights = new vector[weight_t](msize)  
247 - weights.assign(msize, bmin)  
248 -  
249 - for i in prange(n_ids, nogil=True):  
250 - vi_id = vertices_staircase[i]  
251 - deref(weights)[vi_id] = 1.0  
252 -  
253 - vi = &mesh.vertices[vi_id, 0]  
254 - near_vertices = mesh.get_near_vertices_to_v(vi_id, tmax)  
255 - nnv = near_vertices.size()  
256 -  
257 - for j in xrange(nnv):  
258 - vj_id = deref(near_vertices)[j]  
259 - vj = &mesh.vertices[vj_id, 0]  
260 -  
261 - d = sqrt((vi[0] - vj[0]) * (vi[0] - vj[0])\  
262 - + (vi[1] - vj[1]) * (vi[1] - vj[1])\  
263 - + (vi[2] - vj[2]) * (vi[2] - vj[2]))  
264 - value = (1.0 - d/tmax) * (1.0 - bmin) + bmin  
265 -  
266 - if value > deref(weights)[vj_id]:  
267 - deref(weights)[vj_id] = value  
268 -  
269 - del near_vertices  
270 -  
271 - # for i in xrange(msize):  
272 - # if mesh.is_border(i):  
273 - # deref(weights)[i] = 0.0  
274 -  
275 - # cdef vertex_id_t v0, v1, v2  
276 - # for i in xrange(mesh.faces.shape[0]):  
277 - # for j in xrange(1, 4):  
278 - # v0 = mesh.faces[i, j]  
279 - # vi = &mesh.vertices[v0, 0]  
280 - # if mesh.is_border(v0):  
281 - # deref(weights)[v0] = 0.0  
282 - # v1 = mesh.faces[i, (j + 1) % 3 + 1]  
283 - # if mesh.is_border(v1):  
284 - # vi = &mesh.vertices[v1, 0]  
285 - # deref(weights)[v0] = 0.0  
286 -  
287 - return weights  
288 -  
289 -  
290 -cdef inline Point calc_d(Mesh mesh, vertex_id_t v_id) nogil:  
291 - cdef Point D  
292 - cdef int nf, f_id, nid  
293 - cdef float n=0  
294 - cdef int i  
295 - cdef vertex_t* vi  
296 - cdef vertex_t* vj  
297 - cdef set[vertex_id_t]* vertices  
298 - cdef set[vertex_id_t].iterator it  
299 - cdef vertex_id_t vj_id  
300 -  
301 - D.x = 0.0  
302 - D.y = 0.0  
303 - D.z = 0.0  
304 -  
305 - vertices = mesh.get_ring1(v_id)  
306 - vi = &mesh.vertices[v_id, 0]  
307 -  
308 - if mesh.is_border(v_id):  
309 - it = vertices.begin()  
310 - while it != vertices.end():  
311 - vj_id = deref(it)  
312 - if mesh.is_border(vj_id):  
313 - vj = &mesh.vertices[vj_id, 0]  
314 -  
315 - D.x = D.x + (vi[0] - vj[0])  
316 - D.y = D.y + (vi[1] - vj[1])  
317 - D.z = D.z + (vi[2] - vj[2])  
318 - n += 1.0  
319 -  
320 - inc(it)  
321 - else:  
322 - it = vertices.begin()  
323 - while it != vertices.end():  
324 - vj_id = deref(it)  
325 - vj = &mesh.vertices[vj_id, 0]  
326 -  
327 - D.x = D.x + (vi[0] - vj[0])  
328 - D.y = D.y + (vi[1] - vj[1])  
329 - D.z = D.z + (vi[2] - vj[2])  
330 - n += 1.0  
331 -  
332 - inc(it)  
333 -  
334 - del vertices  
335 -  
336 - D.x = D.x / n  
337 - D.y = D.y / n  
338 - D.z = D.z / n  
339 - return D  
340 -  
341 -cdef vector[vertex_id_t]* find_staircase_artifacts(Mesh mesh, double[3] stack_orientation, double T) nogil:  
342 - """  
343 - This function is used to find vertices at staircase artifacts, which are  
344 - those vertices whose incident faces' orientation differences are  
345 - greater than T.  
346 -  
347 - Params:  
348 - mesh: Mesh  
349 - stack_orientation: orientation of slice stacking  
350 - T: Min angle (between vertex faces and stack_orientation) to consider a  
351 - vertex a staircase artifact.  
352 - """  
353 - cdef int nv, nf, f_id, v_id  
354 - cdef double of_z, of_y, of_x, min_z, max_z, min_y, max_y, min_x, max_x;  
355 - cdef vector[vertex_id_t]* f_ids  
356 - cdef normal_t* normal  
357 -  
358 - cdef vector[vertex_id_t]* output = new vector[vertex_id_t]()  
359 - cdef int i  
360 -  
361 - nv = mesh.vertices.shape[0]  
362 -  
363 - for v_id in xrange(nv):  
364 - max_z = -10000  
365 - min_z = 10000  
366 - max_y = -10000  
367 - min_y = 10000  
368 - max_x = -10000  
369 - min_x = 10000  
370 -  
371 - f_ids = mesh.get_faces_by_vertex(v_id)  
372 - nf = deref(f_ids).size()  
373 -  
374 - for i in xrange(nf):  
375 - f_id = deref(f_ids)[i]  
376 - normal = &mesh.normals[f_id, 0]  
377 -  
378 - of_z = 1 - fabs(normal[0]*stack_orientation[0] + normal[1]*stack_orientation[1] + normal[2]*stack_orientation[2]);  
379 - of_y = 1 - fabs(normal[0]*0 + normal[1]*1 + normal[2]*0);  
380 - of_x = 1 - fabs(normal[0]*1 + normal[1]*0 + normal[2]*0);  
381 -  
382 - if (of_z > max_z):  
383 - max_z = of_z  
384 -  
385 - if (of_z < min_z):  
386 - min_z = of_z  
387 -  
388 - if (of_y > max_y):  
389 - max_y = of_y  
390 -  
391 - if (of_y < min_y):  
392 - min_y = of_y  
393 -  
394 - if (of_x > max_x):  
395 - max_x = of_x  
396 -  
397 - if (of_x < min_x):  
398 - min_x = of_x  
399 -  
400 -  
401 - if ((fabs(max_z - min_z) >= T) or (fabs(max_y - min_y) >= T) or (fabs(max_x - min_x) >= T)):  
402 - output.push_back(v_id)  
403 - break  
404 - return output  
405 -  
406 -  
407 -cdef void taubin_smooth(Mesh mesh, vector[weight_t]& weights, float l, float m, int steps):  
408 - """  
409 - Implementation of Taubin's smooth algorithm described in the paper "A  
410 - Signal Processing Approach To Fair Surface Design". His benefeat is it  
411 - avoids surface shrinking.  
412 - """  
413 - cdef int s, i, nvertices  
414 - nvertices = mesh.vertices.shape[0]  
415 - cdef vector[Point] D = vector[Point](nvertices)  
416 - cdef vertex_t* vi  
417 - for s in xrange(steps):  
418 - for i in prange(nvertices, nogil=True):  
419 - D[i] = calc_d(mesh, i)  
420 -  
421 - for i in prange(nvertices, nogil=True):  
422 - mesh.vertices[i, 0] += weights[i]*l*D[i].x;  
423 - mesh.vertices[i, 1] += weights[i]*l*D[i].y;  
424 - mesh.vertices[i, 2] += weights[i]*l*D[i].z;  
425 -  
426 - for i in prange(nvertices, nogil=True):  
427 - D[i] = calc_d(mesh, i)  
428 -  
429 - for i in prange(nvertices, nogil=True):  
430 - mesh.vertices[i, 0] += weights[i]*m*D[i].x;  
431 - mesh.vertices[i, 1] += weights[i]*m*D[i].y;  
432 - mesh.vertices[i, 2] += weights[i]*m*D[i].z;  
433 -  
434 -  
435 -def ca_smoothing(Mesh mesh, double T, double tmax, double bmin, int n_iters):  
436 - """  
437 - This is a implementation of the paper "Context-aware mesh smoothing for  
438 - biomedical applications". It can be used to smooth meshes generated by  
439 - binary images to remove its staircase artifacts and keep the fine features.  
440 -  
441 - Params:  
442 - mesh: Mesh  
443 - T: Min angle (between vertex faces and stack_orientation) to consider a  
444 - vertex a staircase artifact  
445 - tmax: max distance the vertex must be to its nearest artifact vertex  
446 - to considered to calculate the weight  
447 - bmin: The minimum weight  
448 - n_iters: Number of iterations.  
449 - """  
450 - cdef double[3] stack_orientation = [0.0, 0.0, 1.0]  
451 -  
452 - t0 = time.time()  
453 - cdef vector[vertex_id_t]* vertices_staircase = find_staircase_artifacts(mesh, stack_orientation, T)  
454 - print "vertices staircase", time.time() - t0  
455 -  
456 - t0 = time.time()  
457 - cdef vector[weight_t]* weights = calc_artifacts_weight(mesh, deref(vertices_staircase), tmax, bmin)  
458 - print "Weights", time.time() - t0  
459 -  
460 - del vertices_staircase  
461 -  
462 - t0 = time.time()  
463 - taubin_smooth(mesh, deref(weights), 0.5, -0.53, n_iters)  
464 - print "taubin", time.time() - t0  
465 -  
466 - del weights  
invesalius/data/cy_my_types.pxd
View file @73ddbb0
@@ -1,21 +0,0 @@ @@ -1,21 +0,0 @@
1 -import numpy as np  
2 -cimport numpy as np  
3 -cimport cython  
4 -  
5 -# ctypedef np.uint16_t image_t  
6 -  
7 -ctypedef fused image_t:  
8 - np.float64_t  
9 - np.int16_t  
10 - np.uint8_t  
11 -  
12 -ctypedef np.uint8_t mask_t  
13 -  
14 -ctypedef np.float32_t vertex_t  
15 -ctypedef np.float32_t normal_t  
16 -  
17 -# To compile in windows 64  
18 -IF UNAME_MACHINE == 'AMD64':  
19 - ctypedef np.int64_t vertex_id_t  
20 -ELSE:  
21 - ctypedef np.int_t vertex_id_t  
invesalius/data/floodfill.pyx
View file @73ddbb0
@@ -1,274 +0,0 @@ @@ -1,274 +0,0 @@
1 -import numpy as np  
2 -cimport numpy as np  
3 -cimport cython  
4 -  
5 -from collections import deque  
6 -  
7 -from cython.parallel import prange  
8 -from libc.math cimport floor, ceil  
9 -from libcpp cimport bool  
10 -from libcpp.deque cimport deque as cdeque  
11 -from libcpp.vector cimport vector  
12 -  
13 -from cy_my_types cimport image_t, mask_t  
14 -  
15 -cdef struct s_coord:  
16 - int x  
17 - int y  
18 - int z  
19 -  
20 -ctypedef s_coord coord  
21 -  
22 -  
23 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
24 -@cython.wraparound(False)  
25 -@cython.nonecheck(False)  
26 -@cython.cdivision(True)  
27 -cdef inline void append_queue(cdeque[int]& stack, int x, int y, int z, int d, int h, int w) nogil:  
28 - stack.push_back(z*h*w + y*w + x)  
29 -  
30 -  
31 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
32 -@cython.wraparound(False)  
33 -@cython.nonecheck(False)  
34 -@cython.cdivision(True)  
35 -cdef inline void pop_queue(cdeque[int]& stack, int* x, int* y, int* z, int d, int h, int w) nogil:  
36 - cdef int i = stack.front()  
37 - stack.pop_front()  
38 - x[0] = i % w  
39 - y[0] = (i / w) % h  
40 - z[0] = i / (h * w)  
41 -  
42 -  
43 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
44 -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):  
45 -  
46 - cdef int to_return = 0  
47 - if out is None:  
48 - out = np.zeros_like(data)  
49 - to_return = 1  
50 -  
51 - cdef int x, y, z  
52 - cdef int w, h, d  
53 -  
54 - d = data.shape[0]  
55 - h = data.shape[1]  
56 - w = data.shape[2]  
57 -  
58 - stack = [(i, j, k), ]  
59 - out[k, j, i] = fill  
60 -  
61 - while stack:  
62 - x, y, z = stack.pop()  
63 -  
64 - if z + 1 < d and data[z + 1, y, x] == v and out[z + 1, y, x] != fill:  
65 - out[z + 1, y, x] = fill  
66 - stack.append((x, y, z + 1))  
67 -  
68 - if z - 1 >= 0 and data[z - 1, y, x] == v and out[z - 1, y, x] != fill:  
69 - out[z - 1, y, x] = fill  
70 - stack.append((x, y, z - 1))  
71 -  
72 - if y + 1 < h and data[z, y + 1, x] == v and out[z, y + 1, x] != fill:  
73 - out[z, y + 1, x] = fill  
74 - stack.append((x, y + 1, z))  
75 -  
76 - if y - 1 >= 0 and data[z, y - 1, x] == v and out[z, y - 1, x] != fill:  
77 - out[z, y - 1, x] = fill  
78 - stack.append((x, y - 1, z))  
79 -  
80 - if x + 1 < w and data[z, y, x + 1] == v and out[z, y, x + 1] != fill:  
81 - out[z, y, x + 1] = fill  
82 - stack.append((x + 1, y, z))  
83 -  
84 - if x - 1 >= 0 and data[z, y, x - 1] == v and out[z, y, x - 1] != fill:  
85 - out[z, y, x - 1] = fill  
86 - stack.append((x - 1, y, z))  
87 -  
88 - if to_return:  
89 - return out  
90 -  
91 -  
92 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
93 -@cython.wraparound(False)  
94 -@cython.nonecheck(False)  
95 -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):  
96 -  
97 - cdef int to_return = 0  
98 - if out is None:  
99 - out = np.zeros_like(data)  
100 - to_return = 1  
101 -  
102 - cdef int x, y, z  
103 - cdef int dx, dy, dz  
104 - cdef int odx, ody, odz  
105 - cdef int xo, yo, zo  
106 - cdef int i, j, k  
107 - cdef int offset_x, offset_y, offset_z  
108 -  
109 - dz = data.shape[0]  
110 - dy = data.shape[1]  
111 - dx = data.shape[2]  
112 -  
113 - odz = strct.shape[0]  
114 - ody = strct.shape[1]  
115 - odx = strct.shape[2]  
116 -  
117 - cdef cdeque[coord] stack  
118 - cdef coord c  
119 -  
120 - offset_z = odz / 2  
121 - offset_y = ody / 2  
122 - offset_x = odx / 2  
123 -  
124 - for i, j, k in seeds:  
125 - if data[k, j, i] >= t0 and data[k, j, i] <= t1:  
126 - c.x = i  
127 - c.y = j  
128 - c.z = k  
129 - stack.push_back(c)  
130 - out[k, j, i] = fill  
131 -  
132 - with nogil:  
133 - while stack.size():  
134 - c = stack.back()  
135 - stack.pop_back()  
136 -  
137 - x = c.x  
138 - y = c.y  
139 - z = c.z  
140 -  
141 - out[z, y, x] = fill  
142 -  
143 - for k in xrange(odz):  
144 - zo = z + k - offset_z  
145 - for j in xrange(ody):  
146 - yo = y + j - offset_y  
147 - for i in xrange(odx):  
148 - if strct[k, j, i]:  
149 - xo = x + i - offset_x  
150 - 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:  
151 - out[zo, yo, xo] = fill  
152 - c.x = xo  
153 - c.y = yo  
154 - c.z = zo  
155 - stack.push_back(c)  
156 -  
157 - if to_return:  
158 - return out  
159 -  
160 -  
161 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
162 -@cython.wraparound(False)  
163 -@cython.nonecheck(False)  
164 -def floodfill_auto_threshold(np.ndarray[image_t, ndim=3] data, list seeds, float p, int fill, np.ndarray[mask_t, ndim=3] out):  
165 -  
166 - cdef int to_return = 0  
167 - if out is None:  
168 - out = np.zeros_like(data)  
169 - to_return = 1  
170 -  
171 - cdef cdeque[int] stack  
172 - cdef int x, y, z  
173 - cdef int w, h, d  
174 - cdef int xo, yo, zo  
175 - cdef int t0, t1  
176 -  
177 - cdef int i, j, k  
178 -  
179 - d = data.shape[0]  
180 - h = data.shape[1]  
181 - w = data.shape[2]  
182 -  
183 -  
184 - # stack = deque()  
185 -  
186 - x = 0  
187 - y = 0  
188 - z = 0  
189 -  
190 -  
191 - for i, j, k in seeds:  
192 - append_queue(stack, i, j, k, d, h, w)  
193 - out[k, j, i] = fill  
194 - print i, j, k, d, h, w  
195 -  
196 - with nogil:  
197 - while stack.size():  
198 - pop_queue(stack, &x, &y, &z, d, h, w)  
199 -  
200 - # print x, y, z, d, h, w  
201 -  
202 - xo = x  
203 - yo = y  
204 - zo = z  
205 -  
206 - t0 = <int>ceil(data[z, y, x] * (1 - p))  
207 - t1 = <int>floor(data[z, y, x] * (1 + p))  
208 -  
209 - 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:  
210 - out[zo + 1, yo, xo] = fill  
211 - append_queue(stack, x, y, z+1, d, h, w)  
212 -  
213 - 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:  
214 - out[zo - 1, yo, xo] = fill  
215 - append_queue(stack, x, y, z-1, d, h, w)  
216 -  
217 - 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:  
218 - out[zo, yo + 1, xo] = fill  
219 - append_queue(stack, x, y+1, z, d, h, w)  
220 -  
221 - 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:  
222 - out[zo, yo - 1, xo] = fill  
223 - append_queue(stack, x, y-1, z, d, h, w)  
224 -  
225 - 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:  
226 - out[zo, yo, xo + 1] = fill  
227 - append_queue(stack, x+1, y, z, d, h, w)  
228 -  
229 - 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:  
230 - out[zo, yo, xo - 1] = fill  
231 - append_queue(stack, x-1, y, z, d, h, w)  
232 -  
233 - if to_return:  
234 - return out  
235 -  
236 -  
237 -@cython.boundscheck(False)  
238 -@cython.wraparound(False)  
239 -@cython.nonecheck(False)  
240 -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):  
241 - """  
242 - Fill mask holes automatically. The hole must <= max_size. Return True if any hole were filled.  
243 - """  
244 - cdef np.ndarray[np.uint32_t, ndim=1] sizes = np.zeros(shape=(nlabels + 1), dtype=np.uint32)  
245 - cdef int x, y, z  
246 - cdef int dx, dy, dz  
247 - cdef int i  
248 -  
249 - cdef bool modified = False  
250 -  
251 - dz = mask.shape[0]  
252 - dy = mask.shape[1]  
253 - dx = mask.shape[2]  
254 -  
255 - for z in xrange(dz):  
256 - for y in xrange(dy):  
257 - for x in xrange(dx):  
258 - sizes[labels[z, y, x]] += 1  
259 -  
260 - #Checking if any hole will be filled  
261 - for i in xrange(nlabels + 1):  
262 - if sizes[i] <= max_size:  
263 - modified = True  
264 -  
265 - if not modified:  
266 - return 0  
267 -  
268 - for z in prange(dz, nogil=True):  
269 - for y in xrange(dy):  
270 - for x in xrange(dx):  
271 - if sizes[labels[z, y, x]] <= max_size:  
272 - mask[z, y, x] = 254  
273 -  
274 - return modified  
invesalius/data/interpolation.pxd
View file @73ddbb0
@@ -1,8 +0,0 @@ @@ -1,8 +0,0 @@
1 -from .cy_my_types cimport image_t  
2 -  
3 -cdef double interpolate(image_t[:, :, :], double, double, double) nogil  
4 -cdef double tricub_interpolate(image_t[:, :, :], double, double, double) nogil  
5 -cdef double tricubicInterpolate (image_t[:, :, :], double, double, double) nogil  
6 -cdef double lanczos3 (image_t[:, :, :], double, double, double) nogil  
7 -  
8 -cdef double nearest_neighbour_interp(image_t[:, :, :], double, double, double) nogil  
invesalius/data/interpolation.pyx
View file @73ddbb0
@@ -1,392 +0,0 @@ @@ -1,392 +0,0 @@
1 -# from interpolation cimport interpolate  
2 -  
3 -import numpy as np  
4 -cimport numpy as np  
5 -cimport cython  
6 -  
7 -from libc.math cimport floor, ceil, sqrt, fabs, sin, M_PI  
8 -from cython.parallel import prange  
9 -  
10 -DEF LANCZOS_A = 4  
11 -DEF SIZE_LANCZOS_TMP = LANCZOS_A * 2 - 1  
12 -  
13 -cdef double[64][64] temp = [  
14 - [ 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],  
15 - [ 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],  
16 - [-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],  
17 - [ 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],  
18 - [ 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],  
19 - [ 0, 0, 0, 0, 0, 0, 0, 0, 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],  
20 - [ 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],  
21 - [ 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],  
22 - [-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],  
23 - [ 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],  
24 - [ 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],  
25 - [-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],  
26 - [ 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],  
27 - [ 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],  
28 - [-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],  
29 - [ 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],  
30 - [ 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],  
31 - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],  
32 - [ 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],  
33 - [ 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],  
34 - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],  
35 - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],  
36 - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],  
37 - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],  
38 - [ 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],  
39 - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],  
40 - [ 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],  
41 - [ 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],  
42 - [ 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],  
43 - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],  
44 - [ 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],  
45 - [ 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],  
46 - [-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],  
47 - [ 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],  
48 - [ 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],  
49 - [-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],  
50 - [ 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],  
51 - [ 0, 0, 0, 0, 0, 0, 0, 0, 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],  
52 - [ 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],  
53 - [ 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],  
54 - [ 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],  
55 - [ 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],  
56 - [-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],  
57 - [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],  
58 - [-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],  
59 - [ 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],  
60 - [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],  
61 - [-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],  
62 - [ 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],  
63 - [ 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],  
64 - [-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],  
65 - [ 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],  
66 - [ 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],  
67 - [ 0, 0, 0, 0, 0, 0, 0, 0, 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],  
68 - [ 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],  
69 - [ 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],  
70 - [-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],  
71 - [ 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],  
72 - [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],  
73 - [-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],  
74 - [ 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],  
75 - [ 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],  
76 - [-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],  
77 - [ 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]  
78 -]  
79 -  
80 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
81 -@cython.cdivision(True)  
82 -@cython.wraparound(False)  
83 -cdef double nearest_neighbour_interp(image_t[:, :, :] V, double x, double y, double z) nogil:  
84 - return V[<int>(z), <int>(y), <int>(x)]  
85 -  
86 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
87 -@cython.cdivision(True)  
88 -@cython.wraparound(False)  
89 -cdef double interpolate(image_t[:, :, :] V, double x, double y, double z) nogil:  
90 - cdef double xd, yd, zd  
91 - cdef double c00, c10, c01, c11  
92 - cdef double c0, c1  
93 - cdef double c  
94 -  
95 - cdef int x0 = <int>floor(x)  
96 - cdef int x1 = x0 + 1  
97 -  
98 - cdef int y0 = <int>floor(y)  
99 - cdef int y1 = y0 + 1  
100 -  
101 - cdef int z0 = <int>floor(z)  
102 - cdef int z1 = z0 + 1  
103 -  
104 - if x0 == x1:  
105 - xd = 1.0  
106 - else:  
107 - xd = (x - x0) / (x1 - x0)  
108 -  
109 - if y0 == y1:  
110 - yd = 1.0  
111 - else:  
112 - yd = (y - y0) / (y1 - y0)  
113 -  
114 - if z0 == z1:  
115 - zd = 1.0  
116 - else:  
117 - zd = (z - z0) / (z1 - z0)  
118 -  
119 - c00 = _G(V, x0, y0, z0)*(1 - xd) + _G(V, x1, y0, z0)*xd  
120 - c10 = _G(V, x0, y1, z0)*(1 - xd) + _G(V, x1, y1, z0)*xd  
121 - c01 = _G(V, x0, y0, z1)*(1 - xd) + _G(V, x1, y0, z1)*xd  
122 - c11 = _G(V, x0, y1, z1)*(1 - xd) + _G(V, x1, y1, z1)*xd  
123 -  
124 - c0 = c00*(1 - yd) + c10*yd  
125 - c1 = c01*(1 - yd) + c11*yd  
126 -  
127 - c = c0*(1 - zd) + c1*zd  
128 -  
129 - return c  
130 -  
131 -  
132 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
133 -@cython.cdivision(True)  
134 -@cython.wraparound(False)  
135 -cdef inline double lanczos3_L(double x, int a) nogil:  
136 - if x == 0:  
137 - return 1.0  
138 - elif -a <= x < a:  
139 - return (a * sin(M_PI * x) * sin(M_PI * (x / a)))/(M_PI**2 * x**2)  
140 - else:  
141 - return 0.0  
142 -  
143 -  
144 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
145 -@cython.cdivision(True)  
146 -@cython.wraparound(False)  
147 -cdef double lanczos3(image_t[:, :, :] V, double x, double y, double z) nogil:  
148 - cdef int a = LANCZOS_A  
149 -  
150 - cdef int xd = <int>floor(x)  
151 - cdef int yd = <int>floor(y)  
152 - cdef int zd = <int>floor(z)  
153 -  
154 - cdef int xi = xd - a + 1  
155 - cdef int xf = xd + a  
156 -  
157 - cdef int yi = yd - a + 1  
158 - cdef int yf = yd + a  
159 -  
160 - cdef int zi = zd - a + 1  
161 - cdef int zf = zd + a  
162 -  
163 - cdef double lx = 0.0  
164 - cdef double ly = 0.0  
165 - cdef double lz = 0.0  
166 -  
167 - cdef double[SIZE_LANCZOS_TMP][SIZE_LANCZOS_TMP] temp_x  
168 - cdef double[SIZE_LANCZOS_TMP] temp_y  
169 -  
170 - cdef int i, j, k  
171 - cdef int m, n, o  
172 -  
173 - m = 0  
174 - for k in xrange(zi, zf):  
175 - n = 0  
176 - for j in xrange(yi, yf):  
177 - lx = 0  
178 - for i in xrange(xi, xf):  
179 - lx += _G(V, i, j, k) * lanczos3_L(x - i, a)  
180 - temp_x[m][n] = lx  
181 - n += 1  
182 - m += 1  
183 -  
184 - m = 0  
185 - for k in xrange(zi, zf):  
186 - n = 0  
187 - ly = 0  
188 - for j in xrange(yi, yf):  
189 - ly += temp_x[m][n] * lanczos3_L(y - j, a)  
190 - n += 1  
191 - temp_y[m] = ly  
192 - m += 1  
193 -  
194 - m = 0  
195 - for k in xrange(zi, zf):  
196 - lz += temp_y[m] * lanczos3_L(z - k, a)  
197 - m += 1  
198 -  
199 - return lz  
200 -  
201 -  
202 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
203 -@cython.cdivision(True)  
204 -@cython.wraparound(False)  
205 -cdef image_t _G(image_t[:, :, :] V, int x, int y, int z) nogil:  
206 - cdef int dz, dy, dx  
207 - dz = V.shape[0] - 1  
208 - dy = V.shape[1] - 1  
209 - dx = V.shape[2] - 1  
210 -  
211 - if x < 0:  
212 - x = dx + x + 1  
213 - elif x > dx:  
214 - x = x - dx - 1  
215 -  
216 - if y < 0:  
217 - y = dy + y + 1  
218 - elif y > dy:  
219 - y = y - dy - 1  
220 -  
221 - if z < 0:  
222 - z = dz + z + 1  
223 - elif z > dz:  
224 - z = z - dz - 1  
225 -  
226 - return V[z, y, x]  
227 -  
228 -  
229 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
230 -@cython.cdivision(True)  
231 -@cython.wraparound(False)  
232 -cdef void calc_coef_tricub(image_t[:, :, :] V, double x, double y, double z, double [64] coef) nogil:  
233 - cdef int xi = <int>floor(x)  
234 - cdef int yi = <int>floor(y)  
235 - cdef int zi = <int>floor(z)  
236 -  
237 - cdef double[64] _x  
238 -  
239 - cdef int i, j  
240 -  
241 - _x[0] = _G(V, xi, yi, zi)  
242 - _x[1] = _G(V, xi + 1, yi, zi)  
243 - _x[2] = _G(V, xi, yi + 1, zi)  
244 - _x[3] = _G(V, xi + 1, yi + 1, zi)  
245 - _x[4] = _G(V, xi, yi, zi + 1)  
246 - _x[5] = _G(V, xi + 1, yi, zi + 1)  
247 - _x[6] = _G(V, xi, yi + 1, zi + 1)  
248 - _x[7] = _G(V, xi + 1, yi + 1, zi + 1)  
249 -  
250 - _x[8] = 0.5*(_G(V, xi+1,yi,zi) - _G(V, xi-1, yi, zi))  
251 - _x[9] = 0.5*(_G(V, xi+2,yi,zi) - _G(V, xi, yi, zi))  
252 - _x[10] = 0.5*(_G(V, xi+1, yi+1,zi) - _G(V, xi-1, yi+1, zi))  
253 - _x[11] = 0.5*(_G(V, xi+2, yi+1,zi) - _G(V, xi, yi+1, zi))  
254 - _x[12] = 0.5*(_G(V, xi+1, yi,zi+1) - _G(V, xi-1, yi, zi+1))  
255 - _x[13] = 0.5*(_G(V, xi+2, yi,zi+1) - _G(V, xi, yi, zi+1))  
256 - _x[14] = 0.5*(_G(V, xi+1, yi+1,zi+1) - _G(V, xi-1, yi+1, zi+1))  
257 - _x[15] = 0.5*(_G(V, xi+2, yi+1,zi+1) - _G(V, xi, yi+1, zi+1))  
258 - _x[16] = 0.5*(_G(V, xi, yi+1,zi) - _G(V, xi, yi-1, zi))  
259 - _x[17] = 0.5*(_G(V, xi+1, yi+1,zi) - _G(V, xi+1, yi-1, zi))  
260 - _x[18] = 0.5*(_G(V, xi, yi+2,zi) - _G(V, xi, yi, zi))  
261 - _x[19] = 0.5*(_G(V, xi+1, yi+2,zi) - _G(V, xi+1, yi, zi))  
262 - _x[20] = 0.5*(_G(V, xi, yi+1,zi+1) - _G(V, xi, yi-1, zi+1))  
263 - _x[21] = 0.5*(_G(V, xi+1, yi+1,zi+1) - _G(V, xi+1, yi-1, zi+1))  
264 - _x[22] = 0.5*(_G(V, xi, yi+2,zi+1) - _G(V, xi, yi, zi+1))  
265 - _x[23] = 0.5*(_G(V, xi+1, yi+2,zi+1) - _G(V, xi+1, yi, zi+1))  
266 - _x[24] = 0.5*(_G(V, xi, yi,zi+1) - _G(V, xi, yi, zi-1))  
267 - _x[25] = 0.5*(_G(V, xi+1, yi,zi+1) - _G(V, xi+1, yi, zi-1))  
268 - _x[26] = 0.5*(_G(V, xi, yi+1,zi+1) - _G(V, xi, yi+1, zi-1))  
269 - _x[27] = 0.5*(_G(V, xi+1, yi+1,zi+1) - _G(V, xi+1, yi+1, zi-1))  
270 - _x[28] = 0.5*(_G(V, xi, yi,zi+2) - _G(V, xi, yi, zi))  
271 - _x[29] = 0.5*(_G(V, xi+1, yi,zi+2) - _G(V, xi+1, yi, zi))  
272 - _x[30] = 0.5*(_G(V, xi, yi+1,zi+2) - _G(V, xi, yi+1, zi))  
273 - _x[31] = 0.5*(_G(V, xi+1, yi+1,zi+2) - _G(V, xi+1, yi+1, zi))  
274 -  
275 - _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))  
276 - _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))  
277 - _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))  
278 - _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))  
279 - _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))  
280 - _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))  
281 - _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))  
282 - _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))  
283 - _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))  
284 - _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))  
285 - _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))  
286 - _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))  
287 - _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))  
288 - _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))  
289 - _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))  
290 - _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))  
291 - _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))  
292 - _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))  
293 - _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))  
294 - _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))  
295 - _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))  
296 - _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))  
297 - _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))  
298 - _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))  
299 -  
300 - _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))  
301 - _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))  
302 - _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))  
303 - _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))  
304 - _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))  
305 - _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))  
306 - _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))  
307 - _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))  
308 -  
309 - for j in prange(64):  
310 - coef[j] = 0.0  
311 - for i in xrange(64):  
312 - coef[j] += (temp[j][i] * _x[i])  
313 -  
314 -  
315 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
316 -@cython.cdivision(True)  
317 -@cython.wraparound(False)  
318 -cdef double tricub_interpolate(image_t[:, :, :] V, double x, double y, double z) nogil:  
319 - # From: Tricubic interpolation in three dimensions. Lekien and Marsden  
320 - cdef double[64] coef  
321 - cdef double result = 0.0  
322 - calc_coef_tricub(V, x, y, z, coef)  
323 -  
324 - cdef int i, j, k  
325 -  
326 - cdef int xi = <int>floor(x)  
327 - cdef int yi = <int>floor(y)  
328 - cdef int zi = <int>floor(z)  
329 -  
330 - for i in xrange(4):  
331 - for j in xrange(4):  
332 - for k in xrange(4):  
333 - result += (coef[i+4*j+16*k] * ((x-xi)**i) * ((y-yi)**j) * ((z-zi)**k))  
334 - # return V[<int>z, <int>y, <int>x]  
335 - # with gil:  
336 - # print result  
337 - return result  
338 -  
339 -  
340 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
341 -@cython.cdivision(True)  
342 -@cython.wraparound(False)  
343 -cdef double cubicInterpolate(double p[4], double x) nogil:  
344 - 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])))  
345 -  
346 -  
347 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
348 -@cython.cdivision(True)  
349 -@cython.wraparound(False)  
350 -cdef double bicubicInterpolate (double p[4][4], double x, double y) nogil:  
351 - cdef double arr[4]  
352 - arr[0] = cubicInterpolate(p[0], y)  
353 - arr[1] = cubicInterpolate(p[1], y)  
354 - arr[2] = cubicInterpolate(p[2], y)  
355 - arr[3] = cubicInterpolate(p[3], y)  
356 - return cubicInterpolate(arr, x)  
357 -  
358 -  
359 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
360 -@cython.cdivision(True)  
361 -@cython.wraparound(False)  
362 -cdef double tricubicInterpolate(image_t[:, :, :] V, double x, double y, double z) nogil:  
363 - # From http://www.paulinternet.nl/?page=bicubic  
364 - cdef double p[4][4][4]  
365 -  
366 - cdef int xi = <int>floor(x)  
367 - cdef int yi = <int>floor(y)  
368 - cdef int zi = <int>floor(z)  
369 -  
370 - cdef int i, j, k  
371 -  
372 - for i in xrange(4):  
373 - for j in xrange(4):  
374 - for k in xrange(4):  
375 - p[i][j][k] = _G(V, xi + i -1, yi + j -1, zi + k - 1)  
376 -  
377 - cdef double arr[4]  
378 - arr[0] = bicubicInterpolate(p[0], y-yi, z-zi)  
379 - arr[1] = bicubicInterpolate(p[1], y-yi, z-zi)  
380 - arr[2] = bicubicInterpolate(p[2], y-yi, z-zi)  
381 - arr[3] = bicubicInterpolate(p[3], y-yi, z-zi)  
382 - return cubicInterpolate(arr, x-xi)  
383 -  
384 -  
385 -def tricub_interpolate_py(image_t[:, :, :] V, double x, double y, double z):  
386 - return tricub_interpolate(V, x, y, z)  
387 -  
388 -def tricub_interpolate2_py(image_t[:, :, :] V, double x, double y, double z):  
389 - return tricubicInterpolate(V, x, y, z)  
390 -  
391 -def trilin_interpolate_py(image_t[:, :, :] V, double x, double y, double z):  
392 - return interpolate(V, x, y, z)  
invesalius/data/mask.py
  Diff comments   View file @18c2bcd
@@ -30,7 +30,7 @@ import invesalius.constants as const @@ -30,7 +30,7 @@ import invesalius.constants as const
30 import invesalius.data.imagedata_utils as iu 30 import invesalius.data.imagedata_utils as iu
31 import invesalius.session as ses 31 import invesalius.session as ses
32 32
33 -from . import floodfill 33 +from invesalius_cy import floodfill
34 34
35 from wx.lib.pubsub import pub as Publisher 35 from wx.lib.pubsub import pub as Publisher
36 from scipy import ndimage 36 from scipy import ndimage
invesalius/data/mips.pyx
View file @73ddbb0
@@ -1,379 +0,0 @@ @@ -1,379 +0,0 @@
1 -#http://en.wikipedia.org/wiki/Local_maximum_intensity_projection  
2 -import numpy as np  
3 -cimport numpy as np  
4 -cimport cython  
5 -  
6 -from libc.math cimport floor, ceil, sqrt, fabs  
7 -from cython.parallel import prange  
8 -  
9 -DTYPE = np.uint8  
10 -ctypedef np.uint8_t DTYPE_t  
11 -  
12 -DTYPE16 = np.int16  
13 -ctypedef np.int16_t DTYPE16_t  
14 -  
15 -DTYPEF32 = np.float32  
16 -ctypedef np.float32_t DTYPEF32_t  
17 -  
18 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
19 -def lmip(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t tmin,  
20 - DTYPE16_t tmax, np.ndarray[DTYPE16_t, ndim=2] out):  
21 - cdef DTYPE16_t max  
22 - cdef int start  
23 - cdef int sz = image.shape[0]  
24 - cdef int sy = image.shape[1]  
25 - cdef int sx = image.shape[2]  
26 -  
27 - # AXIAL  
28 - if axis == 0:  
29 - for x in xrange(sx):  
30 - for y in xrange(sy):  
31 - max = image[0, y, x]  
32 - if max >= tmin and max <= tmax:  
33 - start = 1  
34 - else:  
35 - start = 0  
36 - for z in xrange(sz):  
37 - if image[z, y, x] > max:  
38 - max = image[z, y, x]  
39 -  
40 - elif image[z, y, x] < max and start:  
41 - break  
42 -  
43 - if image[z, y, x] >= tmin and image[z, y, x] <= tmax:  
44 - start = 1  
45 -  
46 - out[y, x] = max  
47 -  
48 - #CORONAL  
49 - elif axis == 1:  
50 - for z in xrange(sz):  
51 - for x in xrange(sx):  
52 - max = image[z, 0, x]  
53 - if max >= tmin and max <= tmax:  
54 - start = 1  
55 - else:  
56 - start = 0  
57 - for y in xrange(sy):  
58 - if image[z, y, x] > max:  
59 - max = image[z, y, x]  
60 -  
61 - elif image[z, y, x] < max and start:  
62 - break  
63 -  
64 - if image[z, y, x] >= tmin and image[z, y, x] <= tmax:  
65 - start = 1  
66 -  
67 - out[z, x] = max  
68 -  
69 - #CORONAL  
70 - elif axis == 2:  
71 - for z in xrange(sz):  
72 - for y in xrange(sy):  
73 - max = image[z, y, 0]  
74 - if max >= tmin and max <= tmax:  
75 - start = 1  
76 - else:  
77 - start = 0  
78 - for x in xrange(sx):  
79 - if image[z, y, x] > max:  
80 - max = image[z, y, x]  
81 -  
82 - elif image[z, y, x] < max and start:  
83 - break  
84 -  
85 - if image[z, y, x] >= tmin and image[z, y, x] <= tmax:  
86 - start = 1  
87 -  
88 - out[z, y] = max  
89 -  
90 -  
91 -cdef DTYPE16_t get_colour(DTYPE16_t vl, DTYPE16_t wl, DTYPE16_t ww):  
92 - cdef DTYPE16_t out_colour  
93 - cdef DTYPE16_t min_value = wl - (ww / 2)  
94 - cdef DTYPE16_t max_value = wl + (ww / 2)  
95 - if vl < min_value:  
96 - out_colour = min_value  
97 - elif vl > max_value:  
98 - out_colour = max_value  
99 - else:  
100 - out_colour = vl  
101 -  
102 - return out_colour  
103 -  
104 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
105 -@cython.cdivision(True)  
106 -cdef float get_opacity(DTYPE16_t vl, DTYPE16_t wl, DTYPE16_t ww) nogil:  
107 - cdef float out_opacity  
108 - cdef DTYPE16_t min_value = wl - (ww / 2)  
109 - cdef DTYPE16_t max_value = wl + (ww / 2)  
110 - if vl < min_value:  
111 - out_opacity = 0.0  
112 - elif vl > max_value:  
113 - out_opacity = 1.0  
114 - else:  
115 - out_opacity = 1.0/(max_value - min_value) * (vl - min_value)  
116 -  
117 - return out_opacity  
118 -  
119 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
120 -@cython.cdivision(True)  
121 -cdef float get_opacity_f32(DTYPEF32_t vl, DTYPE16_t wl, DTYPE16_t ww) nogil:  
122 - cdef float out_opacity  
123 - cdef DTYPE16_t min_value = wl - (ww / 2)  
124 - cdef DTYPE16_t max_value = wl + (ww / 2)  
125 - if vl < min_value:  
126 - out_opacity = 0.0  
127 - elif vl > max_value:  
128 - out_opacity = 1.0  
129 - else:  
130 - out_opacity = 1.0/(max_value - min_value) * (vl - min_value)  
131 -  
132 - return out_opacity  
133 -  
134 -  
135 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
136 -@cython.cdivision(True)  
137 -def mida(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t wl,  
138 - DTYPE16_t ww, np.ndarray[DTYPE16_t, ndim=2] out):  
139 - cdef int sz = image.shape[0]  
140 - cdef int sy = image.shape[1]  
141 - cdef int sx = image.shape[2]  
142 -  
143 - cdef DTYPE16_t min = image.min()  
144 - cdef DTYPE16_t max = image.max()  
145 - cdef DTYPE16_t vl  
146 -  
147 - cdef DTYPE16_t min_value = wl - (ww / 2)  
148 - cdef DTYPE16_t max_value = wl + (ww / 2)  
149 -  
150 - cdef float fmax=0.0  
151 - cdef float fpi  
152 - cdef float dl  
153 - cdef float bt  
154 -  
155 - cdef float alpha  
156 - cdef float alpha_p = 0.0  
157 - cdef float colour  
158 - cdef float colour_p = 0  
159 -  
160 - cdef int x, y, z  
161 -  
162 - # AXIAL  
163 - if axis == 0:  
164 - for x in prange(sx, nogil=True):  
165 - for y in xrange(sy):  
166 - fmax = 0.0  
167 - alpha_p = 0.0  
168 - colour_p = 0.0  
169 - for z in xrange(sz):  
170 - vl = image[z, y, x]  
171 - fpi = 1.0/(max - min) * (vl - min)  
172 - if fpi > fmax:  
173 - dl = fpi - fmax  
174 - fmax = fpi  
175 - else:  
176 - dl = 0.0  
177 -  
178 - bt = 1.0 - dl  
179 -  
180 - colour = fpi  
181 - alpha = get_opacity(vl, wl, ww)  
182 - colour = (bt * colour_p) + (1 - bt * alpha_p) * colour * alpha  
183 - alpha = (bt * alpha_p) + (1 - bt * alpha_p) * alpha  
184 -  
185 - colour_p = colour  
186 - alpha_p = alpha  
187 -  
188 - if alpha >= 1.0:  
189 - break  
190 -  
191 -  
192 - #out[y, x] = <DTYPE16_t>((max_value - min_value) * colour + min_value)  
193 - out[y, x] = <DTYPE16_t>((max - min) * colour + min)  
194 -  
195 -  
196 - #CORONAL  
197 - elif axis == 1:  
198 - for z in prange(sz, nogil=True):  
199 - for x in xrange(sx):  
200 - fmax = 0.0  
201 - alpha_p = 0.0  
202 - colour_p = 0.0  
203 - for y in xrange(sy):  
204 - vl = image[z, y, x]  
205 - fpi = 1.0/(max - min) * (vl - min)  
206 - if fpi > fmax:  
207 - dl = fpi - fmax  
208 - fmax = fpi  
209 - else:  
210 - dl = 0.0  
211 -  
212 - bt = 1.0 - dl  
213 -  
214 - colour = fpi  
215 - alpha = get_opacity(vl, wl, ww)  
216 - colour = (bt * colour_p) + (1 - bt * alpha_p) * colour * alpha  
217 - alpha = (bt * alpha_p) + (1 - bt * alpha_p) * alpha  
218 -  
219 - colour_p = colour  
220 - alpha_p = alpha  
221 -  
222 - if alpha >= 1.0:  
223 - break  
224 -  
225 - out[z, x] = <DTYPE16_t>((max - min) * colour + min)  
226 -  
227 - #AXIAL  
228 - elif axis == 2:  
229 - for z in prange(sz, nogil=True):  
230 - for y in xrange(sy):  
231 - fmax = 0.0  
232 - alpha_p = 0.0  
233 - colour_p = 0.0  
234 - for x in xrange(sx):  
235 - vl = image[z, y, x]  
236 - fpi = 1.0/(max - min) * (vl - min)  
237 - if fpi > fmax:  
238 - dl = fpi - fmax  
239 - fmax = fpi  
240 - else:  
241 - dl = 0.0  
242 -  
243 - bt = 1.0 - dl  
244 -  
245 - colour = fpi  
246 - alpha = get_opacity(vl, wl, ww)  
247 - colour = (bt * colour_p) + (1 - bt * alpha_p) * colour * alpha  
248 - alpha = (bt * alpha_p) + (1 - bt * alpha_p) * alpha  
249 -  
250 - colour_p = colour  
251 - alpha_p = alpha  
252 -  
253 - if alpha >= 1.0:  
254 - break  
255 -  
256 - out[z, y] = <DTYPE16_t>((max - min) * colour + min)  
257 -  
258 -  
259 -  
260 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
261 -@cython.cdivision(True)  
262 -cdef inline void finite_difference(DTYPE16_t[:, :, :] image,  
263 - int x, int y, int z, float h, float *g) nogil:  
264 - cdef int px, py, pz, fx, fy, fz  
265 -  
266 - cdef int sz = image.shape[0]  
267 - cdef int sy = image.shape[1]  
268 - cdef int sx = image.shape[2]  
269 -  
270 - cdef float gx, gy, gz  
271 -  
272 - if x == 0:  
273 - px = 0  
274 - fx = 1  
275 - elif x == sx - 1:  
276 - px = x - 1  
277 - fx = x  
278 - else:  
279 - px = x - 1  
280 - fx = x + 1  
281 -  
282 - if y == 0:  
283 - py = 0  
284 - fy = 1  
285 - elif y == sy - 1:  
286 - py = y - 1  
287 - fy = y  
288 - else:  
289 - py = y - 1  
290 - fy = y + 1  
291 -  
292 - if z == 0:  
293 - pz = 0  
294 - fz = 1  
295 - elif z == sz - 1:  
296 - pz = z - 1  
297 - fz = z  
298 - else:  
299 - pz = z - 1  
300 - fz = z + 1  
301 -  
302 - gx = (image[z, y, fx] - image[z, y, px]) / (2*h)  
303 - gy = (image[z, fy, x] - image[z, py, x]) / (2*h)  
304 - gz = (image[fz, y, x] - image[pz, y, x]) / (2*h)  
305 -  
306 - g[0] = gx  
307 - g[1] = gy  
308 - g[2] = gz  
309 -  
310 -  
311 -  
312 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
313 -@cython.cdivision(True)  
314 -cdef inline float calc_fcm_itensity(DTYPE16_t[:, :, :] image,  
315 - int x, int y, int z, float n, float* dir) nogil:  
316 - cdef float g[3]  
317 - finite_difference(image, x, y, z, 1.0, g)  
318 - cdef float gm = sqrt(g[0]*g[0] + g[1]*g[1] + g[2]*g[2])  
319 - cdef float d = g[0]*dir[0] + g[1]*dir[1] + g[2]*dir[2]  
320 - cdef float sf = (1.0 - fabs(d/gm))**n  
321 - #alpha = get_opacity_f32(gm, wl, ww)  
322 - cdef float vl = gm * sf  
323 - return vl  
324 -  
325 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
326 -@cython.cdivision(True)  
327 -def fast_countour_mip(np.ndarray[DTYPE16_t, ndim=3] image,  
328 - float n,  
329 - int axis,  
330 - DTYPE16_t wl, DTYPE16_t ww,  
331 - int tmip,  
332 - np.ndarray[DTYPE16_t, ndim=2] out):  
333 - cdef int sz = image.shape[0]  
334 - cdef int sy = image.shape[1]  
335 - cdef int sx = image.shape[2]  
336 - cdef float gm  
337 - cdef float alpha  
338 - cdef float sf  
339 - cdef float d  
340 -  
341 - cdef float* g  
342 - cdef float* dir = [ 0, 0, 0 ]  
343 -  
344 - cdef DTYPE16_t[:, :, :] vimage = image  
345 - cdef np.ndarray[DTYPE16_t, ndim=3] tmp = np.empty_like(image)  
346 -  
347 - cdef DTYPE16_t min = image.min()  
348 - cdef DTYPE16_t max = image.max()  
349 - cdef float fmin = <float>min  
350 - cdef float fmax = <float>max  
351 - cdef float vl  
352 - cdef DTYPE16_t V  
353 -  
354 - cdef int x, y, z  
355 -  
356 - if axis == 0:  
357 - dir[2] = 1.0  
358 - elif axis == 1:  
359 - dir[1] = 1.0  
360 - elif axis == 2:  
361 - dir[0] = 1.0  
362 -  
363 - for z in prange(sz, nogil=True):  
364 - for y in range(sy):  
365 - for x in range(sx):  
366 - vl = calc_fcm_itensity(vimage, x, y, z, n, dir)  
367 - tmp[z, y, x] = <DTYPE16_t>vl  
368 -  
369 - cdef DTYPE16_t tmin = tmp.min()  
370 - cdef DTYPE16_t tmax = tmp.max()  
371 -  
372 - #tmp = ((max - min)/<float>(tmax - tmin)) * (tmp - tmin) + min  
373 -  
374 - if tmip == 0:  
375 - out[:] = tmp.max(axis)  
376 - elif tmip == 1:  
377 - lmip(tmp, axis, 700, 3033, out)  
378 - elif tmip == 2:  
379 - mida(tmp, axis, wl, ww, out)  
invesalius/data/slice_.py
  Diff comments   View file @18c2bcd
@@ -28,13 +28,13 @@ from wx.lib.pubsub import pub as Publisher @@ -28,13 +28,13 @@ from wx.lib.pubsub import pub as Publisher
28 import invesalius.constants as const 28 import invesalius.constants as const
29 import invesalius.data.converters as converters 29 import invesalius.data.converters as converters
30 import invesalius.data.imagedata_utils as iu 30 import invesalius.data.imagedata_utils as iu
31 -import invesalius.data.transformations as transformations  
32 import invesalius.session as ses 31 import invesalius.session as ses
33 import invesalius.style as st 32 import invesalius.style as st
34 import invesalius.utils as utils 33 import invesalius.utils as utils
35 -from invesalius.data import mips, transforms 34 +from invesalius.data import transformations
36 from invesalius.data.mask import Mask 35 from invesalius.data.mask import Mask
37 from invesalius.project import Project 36 from invesalius.project import Project
  37 +from invesalius_cy import mips, transforms
38 38
39 OTHER = 0 39 OTHER = 0
40 PLIST = 1 40 PLIST = 1
invesalius/data/styles.py
  Diff comments   View file @18c2bcd
@@ -46,7 +46,7 @@ import invesalius.utils as utils @@ -46,7 +46,7 @@ import invesalius.utils as utils
46 from invesalius.data.measures import (CircleDensityMeasure, MeasureData, 46 from invesalius.data.measures import (CircleDensityMeasure, MeasureData,
47 PolygonDensityMeasure) 47 PolygonDensityMeasure)
48 48
49 -from . import floodfill 49 +from invesalius_cy import floodfill
50 50
51 ORIENTATIONS = { 51 ORIENTATIONS = {
52 "AXIAL": const.AXIAL, 52 "AXIAL": const.AXIAL,
invesalius/data/surface.py
  Diff comments   View file @18c2bcd
@@ -59,9 +59,9 @@ import invesalius.utils as utl @@ -59,9 +59,9 @@ import invesalius.utils as utl
59 import invesalius.data.vtk_utils as vu 59 import invesalius.data.vtk_utils as vu
60 60
61 from invesalius.gui import dialogs 61 from invesalius.gui import dialogs
  62 +from invesalius_cy import cy_mesh
62 63
63 -from invesalius.data import cy_mesh  
64 -# TODO: Verificar ReleaseDataFlagOn and SetSource 64 +# TODO: Verificar ReleaseDataFlagOn and SetSource
65 65
66 66
67 class Surface(): 67 class Surface():
invesalius/data/surface_process.py
  Diff comments   View file @18c2bcd
@@ -13,7 +13,7 @@ import vtk @@ -13,7 +13,7 @@ import vtk
13 13
14 import invesalius.i18n as i18n 14 import invesalius.i18n as i18n
15 import invesalius.data.converters as converters 15 import invesalius.data.converters as converters
16 -from invesalius.data import cy_mesh 16 +from invesalius_cy import cy_mesh
17 17
18 import weakref 18 import weakref
19 from scipy import ndimage 19 from scipy import ndimage
invesalius/data/transforms.pyx
View file @73ddbb0
@@ -1,174 +0,0 @@ @@ -1,174 +0,0 @@
1 -import numpy as np  
2 -cimport numpy as np  
3 -cimport cython  
4 -  
5 -from .cy_my_types cimport image_t  
6 -from .interpolation cimport interpolate, tricub_interpolate, tricubicInterpolate, lanczos3, nearest_neighbour_interp  
7 -  
8 -from libc.math cimport floor, ceil, sqrt, fabs, round  
9 -from cython.parallel import prange  
10 -  
11 -ctypedef double (*interp_function)(image_t[:, :, :], double, double, double) nogil  
12 -  
13 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
14 -@cython.cdivision(True)  
15 -@cython.wraparound(False)  
16 -cdef void mul_mat4_vec4(double[:, :] M,  
17 - double* coord,  
18 - double* out) nogil:  
19 -  
20 - out[0] = coord[0] * M[0, 0] + coord[1] * M[0, 1] + coord[2] * M[0, 2] + coord[3] * M[0, 3]  
21 - out[1] = coord[0] * M[1, 0] + coord[1] * M[1, 1] + coord[2] * M[1, 2] + coord[3] * M[1, 3]  
22 - out[2] = coord[0] * M[2, 0] + coord[1] * M[2, 1] + coord[2] * M[2, 2] + coord[3] * M[2, 3]  
23 - out[3] = coord[0] * M[3, 0] + coord[1] * M[3, 1] + coord[2] * M[3, 2] + coord[3] * M[3, 3]  
24 -  
25 -  
26 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
27 -@cython.cdivision(True)  
28 -@cython.wraparound(False)  
29 -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:  
30 -  
31 - cdef double coord[4]  
32 - coord[0] = z*sz  
33 - coord[1] = y*sy  
34 - coord[2] = x*sx  
35 - coord[3] = 1.0  
36 -  
37 - cdef double _ncoord[4]  
38 - _ncoord[3] = 1  
39 - # _ncoord[:] = [0.0, 0.0, 0.0, 1.0]  
40 -  
41 - cdef unsigned int dz, dy, dx  
42 - dz = volume.shape[0]  
43 - dy = volume.shape[1]  
44 - dx = volume.shape[2]  
45 -  
46 -  
47 - mul_mat4_vec4(M, coord, _ncoord)  
48 -  
49 - cdef double nz, ny, nx  
50 - nz = (_ncoord[0]/_ncoord[3])/sz  
51 - ny = (_ncoord[1]/_ncoord[3])/sy  
52 - nx = (_ncoord[2]/_ncoord[3])/sx  
53 -  
54 - cdef double v  
55 -  
56 - if 0 <= nz <= (dz-1) and 0 <= ny <= (dy-1) and 0 <= nx <= (dx-1):  
57 - return <image_t>f_interp(volume, nx, ny, nz)  
58 - # if minterpol == 0:  
59 - # return <image_t>nearest_neighbour_interp(volume, nx, ny, nz)  
60 - # elif minterpol == 1:  
61 - # return <image_t>interpolate(volume, nx, ny, nz)  
62 - # elif minterpol == 2:  
63 - # v = tricubicInterpolate(volume, nx, ny, nz)  
64 - # if (v < cval):  
65 - # v = cval  
66 - # return <image_t>v  
67 - # else:  
68 - # v = lanczos3(volume, nx, ny, nz)  
69 - # if (v < cval):  
70 - # v = cval  
71 - # return <image_t>v  
72 - else:  
73 - return cval  
74 -  
75 -  
76 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
77 -@cython.cdivision(True)  
78 -@cython.wraparound(False)  
79 -def apply_view_matrix_transform(image_t[:, :, :] volume,  
80 - spacing,  
81 - double[:, :] M,  
82 - unsigned int n, str orientation,  
83 - int minterpol,  
84 - image_t cval,  
85 - image_t[:, :, :] out):  
86 -  
87 - cdef int dz, dy, dx  
88 - cdef int z, y, x  
89 - dz = volume.shape[0]  
90 - dy = volume.shape[1]  
91 - dx = volume.shape[2]  
92 -  
93 - cdef unsigned int odz, ody, odx  
94 - odz = out.shape[0]  
95 - ody = out.shape[1]  
96 - odx = out.shape[2]  
97 -  
98 - cdef unsigned int count = 0  
99 -  
100 - cdef double sx, sy, sz  
101 - sx = spacing[0]  
102 - sy = spacing[1]  
103 - sz = spacing[2]  
104 -  
105 - cdef interp_function f_interp  
106 -  
107 - if minterpol == 0:  
108 - f_interp = nearest_neighbour_interp  
109 - elif minterpol == 1:  
110 - f_interp = interpolate  
111 - elif minterpol == 2:  
112 - f_interp = tricubicInterpolate  
113 - else:  
114 - f_interp = lanczos3  
115 -  
116 - if orientation == 'AXIAL':  
117 - for z in xrange(n, n+odz):  
118 - for y in prange(dy, nogil=True):  
119 - for x in xrange(dx):  
120 - out[count, y, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval)  
121 - count += 1  
122 -  
123 - elif orientation == 'CORONAL':  
124 - for y in xrange(n, n+ody):  
125 - for z in prange(dz, nogil=True):  
126 - for x in xrange(dx):  
127 - out[z, count, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval)  
128 - count += 1  
129 -  
130 - elif orientation == 'SAGITAL':  
131 - for x in xrange(n, n+odx):  
132 - for z in prange(dz, nogil=True):  
133 - for y in xrange(dy):  
134 - out[z, y, count] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval)  
135 - count += 1  
136 -  
137 -  
138 -@cython.boundscheck(False) # turn of bounds-checking for entire function  
139 -@cython.cdivision(True)  
140 -@cython.wraparound(False)  
141 -def convolve_non_zero(image_t[:, :, :] volume,  
142 - image_t[:, :, :] kernel,  
143 - image_t cval):  
144 - cdef Py_ssize_t x, y, z, sx, sy, sz, kx, ky, kz, skx, sky, skz, i, j, k  
145 - cdef image_t v  
146 -  
147 - cdef image_t[:, :, :] out = np.zeros_like(volume)  
148 -  
149 - sz = volume.shape[0]  
150 - sy = volume.shape[1]  
151 - sx = volume.shape[2]  
152 -  
153 - skz = kernel.shape[0]  
154 - sky = kernel.shape[1]  
155 - skx = kernel.shape[2]  
156 -  
157 - for z in prange(sz, nogil=True):  
158 - for y in xrange(sy):  
159 - for x in xrange(sx):  
160 - if volume[z, y, x] != 0:  
161 - for k in xrange(skz):  
162 - kz = z - skz // 2 + k  
163 - for j in xrange(sky):  
164 - ky = y - sky // 2 + j  
165 - for i in xrange(skx):  
166 - kx = x - skx // 2 + i  
167 -  
168 - if 0 <= kz < sz and 0 <= ky < sy and 0 <= kx < sx:  
169 - v = volume[kz, ky, kx]  
170 - else:  
171 - v = cval  
172 -  
173 - out[z, y, x] += (v * kernel[k, j, i])  
174 - return np.asarray(out)  
invesalius_cy/cy_mesh.pyx 0 → 100644
  Diff comments   View file @18c2bcd
@@ -0,0 +1,466 @@ @@ -0,0 +1,466 @@
  1 +# distutils: language = c++
  2 +#cython: boundscheck=False
  3 +#cython: wraparound=False
  4 +#cython: initializedcheck=False
  5 +#cython: cdivision=True
  6 +#cython: nonecheck=False
  7 +
  8 +import sys
  9 +import time
  10 +cimport numpy as np
  11 +
  12 +from libc.math cimport sin, cos, acos, exp, sqrt, fabs, M_PI
  13 +from libc.stdlib cimport abs as cabs
  14 +from cython.operator cimport dereference as deref, preincrement as inc
  15 +from libcpp.map cimport map
  16 +from libcpp.unordered_map cimport unordered_map
  17 +from libcpp.set cimport set
  18 +from libcpp.vector cimport vector
  19 +from libcpp.pair cimport pair
  20 +from libcpp cimport bool
  21 +from libcpp.deque cimport deque as cdeque
  22 +from cython.parallel import prange
  23 +
  24 +from cy_my_types cimport vertex_t, normal_t, vertex_id_t
  25 +
  26 +import numpy as np
  27 +import vtk
  28 +
  29 +from vtk.util import numpy_support
  30 +
  31 +ctypedef float weight_t
  32 +
  33 +cdef struct Point:
  34 + vertex_t x
  35 + vertex_t y
  36 + vertex_t z
  37 +
  38 +ctypedef pair[vertex_id_t, vertex_id_t] key
  39 +
  40 +
  41 +cdef class Mesh:
  42 + cdef vertex_t[:, :] vertices
  43 + cdef vertex_id_t[:, :] faces
  44 + cdef normal_t[:, :] normals
  45 +
  46 + cdef unordered_map[int, vector[vertex_id_t]] map_vface
  47 + cdef unordered_map[vertex_id_t, int] border_vertices
  48 +
  49 + cdef bool _initialized
  50 +
  51 + def __cinit__(self, pd=None, other=None):
  52 + cdef int i
  53 + cdef map[key, int] edge_nfaces
  54 + cdef map[key, int].iterator it
  55 + if pd:
  56 + self._initialized = True
  57 + _vertices = numpy_support.vtk_to_numpy(pd.GetPoints().GetData())
  58 + _vertices.shape = -1, 3
  59 +
  60 + _faces = numpy_support.vtk_to_numpy(pd.GetPolys().GetData())
  61 + _faces.shape = -1, 4
  62 +
  63 + _normals = numpy_support.vtk_to_numpy(pd.GetCellData().GetArray("Normals"))
  64 + _normals.shape = -1, 3
  65 +
  66 + self.vertices = _vertices
  67 + self.faces = _faces
  68 + self.normals = _normals
  69 +
  70 + for i in xrange(_faces.shape[0]):
  71 + self.map_vface[self.faces[i, 1]].push_back(i)
  72 + self.map_vface[self.faces[i, 2]].push_back(i)
  73 + self.map_vface[self.faces[i, 3]].push_back(i)
  74 +
  75 + edge_nfaces[key(min(self.faces[i, 1], self.faces[i, 2]), max(self.faces[i, 1], self.faces[i, 2]))] += 1
  76 + edge_nfaces[key(min(self.faces[i, 2], self.faces[i, 3]), max(self.faces[i, 2], self.faces[i, 3]))] += 1
  77 + edge_nfaces[key(min(self.faces[i, 1], self.faces[i, 3]), max(self.faces[i, 1], self.faces[i, 3]))] += 1
  78 +
  79 + it = edge_nfaces.begin()
  80 +
  81 + while it != edge_nfaces.end():
  82 + if deref(it).second == 1:
  83 + self.border_vertices[deref(it).first.first] = 1
  84 + self.border_vertices[deref(it).first.second] = 1
  85 +
  86 + inc(it)
  87 +
  88 + elif other:
  89 + _other = <Mesh>other
  90 + self._initialized = True
  91 + self.vertices = _other.vertices.copy()
  92 + self.faces = _other.faces.copy()
  93 + self.normals = _other.normals.copy()
  94 + self.map_vface = unordered_map[int, vector[vertex_id_t]](_other.map_vface)
  95 + self.border_vertices = unordered_map[vertex_id_t, int](_other.border_vertices)
  96 + else:
  97 + self._initialized = False
  98 +
  99 + cdef void copy_to(self, Mesh other):
  100 + """
  101 + Copies self content to other.
  102 + """
  103 + if self._initialized:
  104 + other.vertices[:] = self.vertices
  105 + other.faces[:] = self.faces
  106 + other.normals[:] = self.normals
  107 + other.map_vface = unordered_map[int, vector[vertex_id_t]](self.map_vface)
  108 + other.border_vertices = unordered_map[vertex_id_t, int](self.border_vertices)
  109 + else:
  110 + other.vertices = self.vertices.copy()
  111 + other.faces = self.faces.copy()
  112 + other.normals = self.normals.copy()
  113 +
  114 + other.map_vface = self.map_vface
  115 + other.border_vertices = self.border_vertices
  116 +
  117 + def to_vtk(self):
  118 + """
  119 + Converts Mesh to vtkPolyData.
  120 + """
  121 + vertices = np.asarray(self.vertices)
  122 + faces = np.asarray(self.faces)
  123 + normals = np.asarray(self.normals)
  124 +
  125 + points = vtk.vtkPoints()
  126 + points.SetData(numpy_support.numpy_to_vtk(vertices))
  127 +
  128 + id_triangles = numpy_support.numpy_to_vtkIdTypeArray(faces)
  129 + triangles = vtk.vtkCellArray()
  130 + triangles.SetCells(faces.shape[0], id_triangles)
  131 +
  132 + pd = vtk.vtkPolyData()
  133 + pd.SetPoints(points)
  134 + pd.SetPolys(triangles)
  135 +
  136 + return pd
  137 +
  138 + cdef vector[vertex_id_t]* get_faces_by_vertex(self, int v_id) nogil:
  139 + """
  140 + Returns the faces whose vertex `v_id' is part.
  141 + """
  142 + return &self.map_vface[v_id]
  143 +
  144 + cdef set[vertex_id_t]* get_ring1(self, vertex_id_t v_id) nogil:
  145 + """
  146 + Returns the ring1 of vertex `v_id'
  147 + """
  148 + cdef vertex_id_t f_id
  149 + cdef set[vertex_id_t]* ring1 = new set[vertex_id_t]()
  150 + cdef vector[vertex_id_t].iterator it = self.map_vface[v_id].begin()
  151 +
  152 + while it != self.map_vface[v_id].end():
  153 + f_id = deref(it)
  154 + inc(it)
  155 + if self.faces[f_id, 1] != v_id:
  156 + ring1.insert(self.faces[f_id, 1])
  157 + if self.faces[f_id, 2] != v_id:
  158 + ring1.insert(self.faces[f_id, 2])
  159 + if self.faces[f_id, 3] != v_id:
  160 + ring1.insert(self.faces[f_id, 3])
  161 +
  162 + return ring1
  163 +
  164 + cdef bool is_border(self, vertex_id_t v_id) nogil:
  165 + """
  166 + Check if vertex `v_id' is a vertex border.
  167 + """
  168 + return self.border_vertices.find(v_id) != self.border_vertices.end()
  169 +
  170 + cdef vector[vertex_id_t]* get_near_vertices_to_v(self, vertex_id_t v_id, float dmax) nogil:
  171 + """
  172 + Returns all vertices with distance at most `d' to the vertex `v_id'
  173 +
  174 + Params:
  175 + v_id: id of the vertex
  176 + dmax: the maximum distance.
  177 + """
  178 + cdef vector[vertex_id_t]* idfaces
  179 + cdef vector[vertex_id_t]* near_vertices = new vector[vertex_id_t]()
  180 +
  181 + cdef cdeque[vertex_id_t] to_visit
  182 + cdef unordered_map[vertex_id_t, bool] status_v
  183 + cdef unordered_map[vertex_id_t, bool] status_f
  184 +
  185 + cdef vertex_t *vip
  186 + cdef vertex_t *vjp
  187 +
  188 + cdef float distance
  189 + cdef int nf, nid, j
  190 + cdef vertex_id_t f_id, vj
  191 +
  192 + vip = &self.vertices[v_id, 0]
  193 + to_visit.push_back(v_id)
  194 + while(not to_visit.empty()):
  195 + v_id = to_visit.front()
  196 + to_visit.pop_front()
  197 +
  198 + status_v[v_id] = True
  199 +
  200 + idfaces = self.get_faces_by_vertex(v_id)
  201 + nf = idfaces.size()
  202 +
  203 + for nid in xrange(nf):
  204 + f_id = deref(idfaces)[nid]
  205 + if status_f.find(f_id) == status_f.end():
  206 + status_f[f_id] = True
  207 +
  208 + for j in xrange(3):
  209 + vj = self.faces[f_id, j+1]
  210 + if status_v.find(vj) == status_v.end():
  211 + status_v[vj] = True
  212 + vjp = &self.vertices[vj, 0]
  213 + distance = sqrt((vip[0] - vjp[0]) * (vip[0] - vjp[0]) \
  214 + + (vip[1] - vjp[1]) * (vip[1] - vjp[1]) \
  215 + + (vip[2] - vjp[2]) * (vip[2] - vjp[2]))
  216 + if distance <= dmax:
  217 + near_vertices.push_back(vj)
  218 + to_visit.push_back(vj)
  219 +
  220 + return near_vertices
  221 +
  222 +
  223 +cdef vector[weight_t]* calc_artifacts_weight(Mesh mesh, vector[vertex_id_t]& vertices_staircase, float tmax, float bmin) nogil:
  224 + """
  225 + Calculate the artifact weight based on distance of each vertex to its
  226 + nearest staircase artifact vertex.
  227 +
  228 + Params:
  229 + mesh: Mesh
  230 + vertices_staircase: the identified staircase artifact vertices
  231 + tmax: max distance the vertex must be to its nearest artifact vertex
  232 + to considered to calculate the weight
  233 + bmin: The minimum weight.
  234 + """
  235 + cdef int vi_id, vj_id, nnv, n_ids, i, j
  236 + cdef vector[vertex_id_t]* near_vertices
  237 + cdef weight_t value
  238 + cdef float d
  239 + n_ids = vertices_staircase.size()
  240 +
  241 + cdef vertex_t* vi
  242 + cdef vertex_t* vj
  243 + cdef size_t msize
  244 +
  245 + msize = mesh.vertices.shape[0]
  246 + cdef vector[weight_t]* weights = new vector[weight_t](msize)
  247 + weights.assign(msize, bmin)
  248 +
  249 + for i in prange(n_ids, nogil=True):
  250 + vi_id = vertices_staircase[i]
  251 + deref(weights)[vi_id] = 1.0
  252 +
  253 + vi = &mesh.vertices[vi_id, 0]
  254 + near_vertices = mesh.get_near_vertices_to_v(vi_id, tmax)
  255 + nnv = near_vertices.size()
  256 +
  257 + for j in xrange(nnv):
  258 + vj_id = deref(near_vertices)[j]
  259 + vj = &mesh.vertices[vj_id, 0]
  260 +
  261 + d = sqrt((vi[0] - vj[0]) * (vi[0] - vj[0])\
  262 + + (vi[1] - vj[1]) * (vi[1] - vj[1])\
  263 + + (vi[2] - vj[2]) * (vi[2] - vj[2]))
  264 + value = (1.0 - d/tmax) * (1.0 - bmin) + bmin
  265 +
  266 + if value > deref(weights)[vj_id]:
  267 + deref(weights)[vj_id] = value
  268 +
  269 + del near_vertices
  270 +
  271 + # for i in xrange(msize):
  272 + # if mesh.is_border(i):
  273 + # deref(weights)[i] = 0.0
  274 +
  275 + # cdef vertex_id_t v0, v1, v2
  276 + # for i in xrange(mesh.faces.shape[0]):
  277 + # for j in xrange(1, 4):
  278 + # v0 = mesh.faces[i, j]
  279 + # vi = &mesh.vertices[v0, 0]
  280 + # if mesh.is_border(v0):
  281 + # deref(weights)[v0] = 0.0
  282 + # v1 = mesh.faces[i, (j + 1) % 3 + 1]
  283 + # if mesh.is_border(v1):
  284 + # vi = &mesh.vertices[v1, 0]
  285 + # deref(weights)[v0] = 0.0
  286 +
  287 + return weights
  288 +
  289 +
  290 +cdef inline Point calc_d(Mesh mesh, vertex_id_t v_id) nogil:
  291 + cdef Point D
  292 + cdef int nf, f_id, nid
  293 + cdef float n=0
  294 + cdef int i
  295 + cdef vertex_t* vi
  296 + cdef vertex_t* vj
  297 + cdef set[vertex_id_t]* vertices
  298 + cdef set[vertex_id_t].iterator it
  299 + cdef vertex_id_t vj_id
  300 +
  301 + D.x = 0.0
  302 + D.y = 0.0
  303 + D.z = 0.0
  304 +
  305 + vertices = mesh.get_ring1(v_id)
  306 + vi = &mesh.vertices[v_id, 0]
  307 +
  308 + if mesh.is_border(v_id):
  309 + it = vertices.begin()
  310 + while it != vertices.end():
  311 + vj_id = deref(it)
  312 + if mesh.is_border(vj_id):
  313 + vj = &mesh.vertices[vj_id, 0]
  314 +
  315 + D.x = D.x + (vi[0] - vj[0])
  316 + D.y = D.y + (vi[1] - vj[1])
  317 + D.z = D.z + (vi[2] - vj[2])
  318 + n += 1.0
  319 +
  320 + inc(it)
  321 + else:
  322 + it = vertices.begin()
  323 + while it != vertices.end():
  324 + vj_id = deref(it)
  325 + vj = &mesh.vertices[vj_id, 0]
  326 +
  327 + D.x = D.x + (vi[0] - vj[0])
  328 + D.y = D.y + (vi[1] - vj[1])
  329 + D.z = D.z + (vi[2] - vj[2])
  330 + n += 1.0
  331 +
  332 + inc(it)
  333 +
  334 + del vertices
  335 +
  336 + D.x = D.x / n
  337 + D.y = D.y / n
  338 + D.z = D.z / n
  339 + return D
  340 +
  341 +cdef vector[vertex_id_t]* find_staircase_artifacts(Mesh mesh, double[3] stack_orientation, double T) nogil:
  342 + """
  343 + This function is used to find vertices at staircase artifacts, which are
  344 + those vertices whose incident faces' orientation differences are
  345 + greater than T.
  346 +
  347 + Params:
  348 + mesh: Mesh
  349 + stack_orientation: orientation of slice stacking
  350 + T: Min angle (between vertex faces and stack_orientation) to consider a
  351 + vertex a staircase artifact.
  352 + """
  353 + cdef int nv, nf, f_id, v_id
  354 + cdef double of_z, of_y, of_x, min_z, max_z, min_y, max_y, min_x, max_x;
  355 + cdef vector[vertex_id_t]* f_ids
  356 + cdef normal_t* normal
  357 +
  358 + cdef vector[vertex_id_t]* output = new vector[vertex_id_t]()
  359 + cdef int i
  360 +
  361 + nv = mesh.vertices.shape[0]
  362 +
  363 + for v_id in xrange(nv):
  364 + max_z = -10000
  365 + min_z = 10000
  366 + max_y = -10000
  367 + min_y = 10000
  368 + max_x = -10000
  369 + min_x = 10000
  370 +
  371 + f_ids = mesh.get_faces_by_vertex(v_id)
  372 + nf = deref(f_ids).size()
  373 +
  374 + for i in xrange(nf):
  375 + f_id = deref(f_ids)[i]
  376 + normal = &mesh.normals[f_id, 0]
  377 +
  378 + of_z = 1 - fabs(normal[0]*stack_orientation[0] + normal[1]*stack_orientation[1] + normal[2]*stack_orientation[2]);
  379 + of_y = 1 - fabs(normal[0]*0 + normal[1]*1 + normal[2]*0);
  380 + of_x = 1 - fabs(normal[0]*1 + normal[1]*0 + normal[2]*0);
  381 +
  382 + if (of_z > max_z):
  383 + max_z = of_z
  384 +
  385 + if (of_z < min_z):
  386 + min_z = of_z
  387 +
  388 + if (of_y > max_y):
  389 + max_y = of_y
  390 +
  391 + if (of_y < min_y):
  392 + min_y = of_y
  393 +
  394 + if (of_x > max_x):
  395 + max_x = of_x
  396 +
  397 + if (of_x < min_x):
  398 + min_x = of_x
  399 +
  400 +
  401 + if ((fabs(max_z - min_z) >= T) or (fabs(max_y - min_y) >= T) or (fabs(max_x - min_x) >= T)):
  402 + output.push_back(v_id)
  403 + break
  404 + return output
  405 +
  406 +
  407 +cdef void taubin_smooth(Mesh mesh, vector[weight_t]& weights, float l, float m, int steps):
  408 + """
  409 + Implementation of Taubin's smooth algorithm described in the paper "A
  410 + Signal Processing Approach To Fair Surface Design". His benefeat is it
  411 + avoids surface shrinking.
  412 + """
  413 + cdef int s, i, nvertices
  414 + nvertices = mesh.vertices.shape[0]
  415 + cdef vector[Point] D = vector[Point](nvertices)
  416 + cdef vertex_t* vi
  417 + for s in xrange(steps):
  418 + for i in prange(nvertices, nogil=True):
  419 + D[i] = calc_d(mesh, i)
  420 +
  421 + for i in prange(nvertices, nogil=True):
  422 + mesh.vertices[i, 0] += weights[i]*l*D[i].x;
  423 + mesh.vertices[i, 1] += weights[i]*l*D[i].y;
  424 + mesh.vertices[i, 2] += weights[i]*l*D[i].z;
  425 +
  426 + for i in prange(nvertices, nogil=True):
  427 + D[i] = calc_d(mesh, i)
  428 +
  429 + for i in prange(nvertices, nogil=True):
  430 + mesh.vertices[i, 0] += weights[i]*m*D[i].x;
  431 + mesh.vertices[i, 1] += weights[i]*m*D[i].y;
  432 + mesh.vertices[i, 2] += weights[i]*m*D[i].z;
  433 +
  434 +
  435 +def ca_smoothing(Mesh mesh, double T, double tmax, double bmin, int n_iters):
  436 + """
  437 + This is a implementation of the paper "Context-aware mesh smoothing for
  438 + biomedical applications". It can be used to smooth meshes generated by
  439 + binary images to remove its staircase artifacts and keep the fine features.
  440 +
  441 + Params:
  442 + mesh: Mesh
  443 + T: Min angle (between vertex faces and stack_orientation) to consider a
  444 + vertex a staircase artifact
  445 + tmax: max distance the vertex must be to its nearest artifact vertex
  446 + to considered to calculate the weight
  447 + bmin: The minimum weight
  448 + n_iters: Number of iterations.
  449 + """
  450 + cdef double[3] stack_orientation = [0.0, 0.0, 1.0]
  451 +
  452 + t0 = time.time()
  453 + cdef vector[vertex_id_t]* vertices_staircase = find_staircase_artifacts(mesh, stack_orientation, T)
  454 + print "vertices staircase", time.time() - t0
  455 +
  456 + t0 = time.time()
  457 + cdef vector[weight_t]* weights = calc_artifacts_weight(mesh, deref(vertices_staircase), tmax, bmin)
  458 + print "Weights", time.time() - t0
  459 +
  460 + del vertices_staircase
  461 +
  462 + t0 = time.time()
  463 + taubin_smooth(mesh, deref(weights), 0.5, -0.53, n_iters)
  464 + print "taubin", time.time() - t0
  465 +
  466 + del weights
invesalius_cy/cy_my_types.pxd 0 → 100644
  Diff comments   View file @18c2bcd
@@ -0,0 +1,21 @@ @@ -0,0 +1,21 @@
  1 +import numpy as np
  2 +cimport numpy as np
  3 +cimport cython
  4 +
  5 +# ctypedef np.uint16_t image_t
  6 +
  7 +ctypedef fused image_t:
  8 + np.float64_t
  9 + np.int16_t
  10 + np.uint8_t
  11 +
  12 +ctypedef np.uint8_t mask_t
  13 +
  14 +ctypedef np.float32_t vertex_t
  15 +ctypedef np.float32_t normal_t
  16 +
  17 +# To compile in windows 64
  18 +IF UNAME_MACHINE == 'AMD64':
  19 + ctypedef np.int64_t vertex_id_t
  20 +ELSE:
  21 + ctypedef np.int_t vertex_id_t
invesalius_cy/floodfill.pyx 0 → 100644
  Diff comments   View file @18c2bcd
@@ -0,0 +1,274 @@ @@ -0,0 +1,274 @@
  1 +import numpy as np
  2 +cimport numpy as np
  3 +cimport cython
  4 +
  5 +from collections import deque
  6 +
  7 +from cython.parallel import prange
  8 +from libc.math cimport floor, ceil
  9 +from libcpp cimport bool
  10 +from libcpp.deque cimport deque as cdeque
  11 +from libcpp.vector cimport vector
  12 +
  13 +from cy_my_types cimport image_t, mask_t
  14 +
  15 +cdef struct s_coord:
  16 + int x
  17 + int y
  18 + int z
  19 +
  20 +ctypedef s_coord coord
  21 +
  22 +
  23 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  24 +@cython.wraparound(False)
  25 +@cython.nonecheck(False)
  26 +@cython.cdivision(True)
  27 +cdef inline void append_queue(cdeque[int]& stack, int x, int y, int z, int d, int h, int w) nogil:
  28 + stack.push_back(z*h*w + y*w + x)
  29 +
  30 +
  31 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  32 +@cython.wraparound(False)
  33 +@cython.nonecheck(False)
  34 +@cython.cdivision(True)
  35 +cdef inline void pop_queue(cdeque[int]& stack, int* x, int* y, int* z, int d, int h, int w) nogil:
  36 + cdef int i = stack.front()
  37 + stack.pop_front()
  38 + x[0] = i % w
  39 + y[0] = (i / w) % h
  40 + z[0] = i / (h * w)
  41 +
  42 +
  43 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  44 +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):
  45 +
  46 + cdef int to_return = 0
  47 + if out is None:
  48 + out = np.zeros_like(data)
  49 + to_return = 1
  50 +
  51 + cdef int x, y, z
  52 + cdef int w, h, d
  53 +
  54 + d = data.shape[0]
  55 + h = data.shape[1]
  56 + w = data.shape[2]
  57 +
  58 + stack = [(i, j, k), ]
  59 + out[k, j, i] = fill
  60 +
  61 + while stack:
  62 + x, y, z = stack.pop()
  63 +
  64 + if z + 1 < d and data[z + 1, y, x] == v and out[z + 1, y, x] != fill:
  65 + out[z + 1, y, x] = fill
  66 + stack.append((x, y, z + 1))
  67 +
  68 + if z - 1 >= 0 and data[z - 1, y, x] == v and out[z - 1, y, x] != fill:
  69 + out[z - 1, y, x] = fill
  70 + stack.append((x, y, z - 1))
  71 +
  72 + if y + 1 < h and data[z, y + 1, x] == v and out[z, y + 1, x] != fill:
  73 + out[z, y + 1, x] = fill
  74 + stack.append((x, y + 1, z))
  75 +
  76 + if y - 1 >= 0 and data[z, y - 1, x] == v and out[z, y - 1, x] != fill:
  77 + out[z, y - 1, x] = fill
  78 + stack.append((x, y - 1, z))
  79 +
  80 + if x + 1 < w and data[z, y, x + 1] == v and out[z, y, x + 1] != fill:
  81 + out[z, y, x + 1] = fill
  82 + stack.append((x + 1, y, z))
  83 +
  84 + if x - 1 >= 0 and data[z, y, x - 1] == v and out[z, y, x - 1] != fill:
  85 + out[z, y, x - 1] = fill
  86 + stack.append((x - 1, y, z))
  87 +
  88 + if to_return:
  89 + return out
  90 +
  91 +
  92 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  93 +@cython.wraparound(False)
  94 +@cython.nonecheck(False)
  95 +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):
  96 +
  97 + cdef int to_return = 0
  98 + if out is None:
  99 + out = np.zeros_like(data)
  100 + to_return = 1
  101 +
  102 + cdef int x, y, z
  103 + cdef int dx, dy, dz
  104 + cdef int odx, ody, odz
  105 + cdef int xo, yo, zo
  106 + cdef int i, j, k
  107 + cdef int offset_x, offset_y, offset_z
  108 +
  109 + dz = data.shape[0]
  110 + dy = data.shape[1]
  111 + dx = data.shape[2]
  112 +
  113 + odz = strct.shape[0]
  114 + ody = strct.shape[1]
  115 + odx = strct.shape[2]
  116 +
  117 + cdef cdeque[coord] stack
  118 + cdef coord c
  119 +
  120 + offset_z = odz / 2
  121 + offset_y = ody / 2
  122 + offset_x = odx / 2
  123 +
  124 + for i, j, k in seeds:
  125 + if data[k, j, i] >= t0 and data[k, j, i] <= t1:
  126 + c.x = i
  127 + c.y = j
  128 + c.z = k
  129 + stack.push_back(c)
  130 + out[k, j, i] = fill
  131 +
  132 + with nogil:
  133 + while stack.size():
  134 + c = stack.back()
  135 + stack.pop_back()
  136 +
  137 + x = c.x
  138 + y = c.y
  139 + z = c.z
  140 +
  141 + out[z, y, x] = fill
  142 +
  143 + for k in xrange(odz):
  144 + zo = z + k - offset_z
  145 + for j in xrange(ody):
  146 + yo = y + j - offset_y
  147 + for i in xrange(odx):
  148 + if strct[k, j, i]:
  149 + xo = x + i - offset_x
  150 + 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:
  151 + out[zo, yo, xo] = fill
  152 + c.x = xo
  153 + c.y = yo
  154 + c.z = zo
  155 + stack.push_back(c)
  156 +
  157 + if to_return:
  158 + return out
  159 +
  160 +
  161 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  162 +@cython.wraparound(False)
  163 +@cython.nonecheck(False)
  164 +def floodfill_auto_threshold(np.ndarray[image_t, ndim=3] data, list seeds, float p, int fill, np.ndarray[mask_t, ndim=3] out):
  165 +
  166 + cdef int to_return = 0
  167 + if out is None:
  168 + out = np.zeros_like(data)
  169 + to_return = 1
  170 +
  171 + cdef cdeque[int] stack
  172 + cdef int x, y, z
  173 + cdef int w, h, d
  174 + cdef int xo, yo, zo
  175 + cdef int t0, t1
  176 +
  177 + cdef int i, j, k
  178 +
  179 + d = data.shape[0]
  180 + h = data.shape[1]
  181 + w = data.shape[2]
  182 +
  183 +
  184 + # stack = deque()
  185 +
  186 + x = 0
  187 + y = 0
  188 + z = 0
  189 +
  190 +
  191 + for i, j, k in seeds:
  192 + append_queue(stack, i, j, k, d, h, w)
  193 + out[k, j, i] = fill
  194 + print i, j, k, d, h, w
  195 +
  196 + with nogil:
  197 + while stack.size():
  198 + pop_queue(stack, &x, &y, &z, d, h, w)
  199 +
  200 + # print x, y, z, d, h, w
  201 +
  202 + xo = x
  203 + yo = y
  204 + zo = z
  205 +
  206 + t0 = <int>ceil(data[z, y, x] * (1 - p))
  207 + t1 = <int>floor(data[z, y, x] * (1 + p))
  208 +
  209 + 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:
  210 + out[zo + 1, yo, xo] = fill
  211 + append_queue(stack, x, y, z+1, d, h, w)
  212 +
  213 + 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:
  214 + out[zo - 1, yo, xo] = fill
  215 + append_queue(stack, x, y, z-1, d, h, w)
  216 +
  217 + 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:
  218 + out[zo, yo + 1, xo] = fill
  219 + append_queue(stack, x, y+1, z, d, h, w)
  220 +
  221 + 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:
  222 + out[zo, yo - 1, xo] = fill
  223 + append_queue(stack, x, y-1, z, d, h, w)
  224 +
  225 + 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:
  226 + out[zo, yo, xo + 1] = fill
  227 + append_queue(stack, x+1, y, z, d, h, w)
  228 +
  229 + 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:
  230 + out[zo, yo, xo - 1] = fill
  231 + append_queue(stack, x-1, y, z, d, h, w)
  232 +
  233 + if to_return:
  234 + return out
  235 +
  236 +
  237 +@cython.boundscheck(False)
  238 +@cython.wraparound(False)
  239 +@cython.nonecheck(False)
  240 +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):
  241 + """
  242 + Fill mask holes automatically. The hole must <= max_size. Return True if any hole were filled.
  243 + """
  244 + cdef np.ndarray[np.uint32_t, ndim=1] sizes = np.zeros(shape=(nlabels + 1), dtype=np.uint32)
  245 + cdef int x, y, z
  246 + cdef int dx, dy, dz
  247 + cdef int i
  248 +
  249 + cdef bool modified = False
  250 +
  251 + dz = mask.shape[0]
  252 + dy = mask.shape[1]
  253 + dx = mask.shape[2]
  254 +
  255 + for z in xrange(dz):
  256 + for y in xrange(dy):
  257 + for x in xrange(dx):
  258 + sizes[labels[z, y, x]] += 1
  259 +
  260 + #Checking if any hole will be filled
  261 + for i in xrange(nlabels + 1):
  262 + if sizes[i] <= max_size:
  263 + modified = True
  264 +
  265 + if not modified:
  266 + return 0
  267 +
  268 + for z in prange(dz, nogil=True):
  269 + for y in xrange(dy):
  270 + for x in xrange(dx):
  271 + if sizes[labels[z, y, x]] <= max_size:
  272 + mask[z, y, x] = 254
  273 +
  274 + return modified
invesalius_cy/interpolation.pxd 0 → 100644
  Diff comments   View file @18c2bcd
@@ -0,0 +1,8 @@ @@ -0,0 +1,8 @@
  1 +from .cy_my_types cimport image_t
  2 +
  3 +cdef double interpolate(image_t[:, :, :], double, double, double) nogil
  4 +cdef double tricub_interpolate(image_t[:, :, :], double, double, double) nogil
  5 +cdef double tricubicInterpolate (image_t[:, :, :], double, double, double) nogil
  6 +cdef double lanczos3 (image_t[:, :, :], double, double, double) nogil
  7 +
  8 +cdef double nearest_neighbour_interp(image_t[:, :, :], double, double, double) nogil
invesalius_cy/interpolation.pyx 0 → 100644
  Diff comments   View file @18c2bcd
@@ -0,0 +1,392 @@ @@ -0,0 +1,392 @@
  1 +# from interpolation cimport interpolate
  2 +
  3 +import numpy as np
  4 +cimport numpy as np
  5 +cimport cython
  6 +
  7 +from libc.math cimport floor, ceil, sqrt, fabs, sin, M_PI
  8 +from cython.parallel import prange
  9 +
  10 +DEF LANCZOS_A = 4
  11 +DEF SIZE_LANCZOS_TMP = LANCZOS_A * 2 - 1
  12 +
  13 +cdef double[64][64] temp = [
  14 + [ 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],
  15 + [ 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],
  16 + [-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],
  17 + [ 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],
  18 + [ 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],
  19 + [ 0, 0, 0, 0, 0, 0, 0, 0, 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],
  20 + [ 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],
  21 + [ 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],
  22 + [-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],
  23 + [ 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],
  24 + [ 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],
  25 + [-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],
  26 + [ 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],
  27 + [ 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],
  28 + [-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],
  29 + [ 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],
  30 + [ 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],
  31 + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],
  32 + [ 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],
  33 + [ 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],
  34 + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],
  35 + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],
  36 + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],
  37 + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],
  38 + [ 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],
  39 + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],
  40 + [ 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],
  41 + [ 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],
  42 + [ 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],
  43 + [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 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],
  44 + [ 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],
  45 + [ 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],
  46 + [-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],
  47 + [ 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],
  48 + [ 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],
  49 + [-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],
  50 + [ 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],
  51 + [ 0, 0, 0, 0, 0, 0, 0, 0, 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],
  52 + [ 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],
  53 + [ 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],
  54 + [ 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],
  55 + [ 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],
  56 + [-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],
  57 + [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],
  58 + [-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],
  59 + [ 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],
  60 + [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],
  61 + [-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],
  62 + [ 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],
  63 + [ 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],
  64 + [-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],
  65 + [ 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],
  66 + [ 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],
  67 + [ 0, 0, 0, 0, 0, 0, 0, 0, 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],
  68 + [ 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],
  69 + [ 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],
  70 + [-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],
  71 + [ 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],
  72 + [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],
  73 + [-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],
  74 + [ 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],
  75 + [ 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],
  76 + [-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],
  77 + [ 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]
  78 +]
  79 +
  80 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  81 +@cython.cdivision(True)
  82 +@cython.wraparound(False)
  83 +cdef double nearest_neighbour_interp(image_t[:, :, :] V, double x, double y, double z) nogil:
  84 + return V[<int>(z), <int>(y), <int>(x)]
  85 +
  86 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  87 +@cython.cdivision(True)
  88 +@cython.wraparound(False)
  89 +cdef double interpolate(image_t[:, :, :] V, double x, double y, double z) nogil:
  90 + cdef double xd, yd, zd
  91 + cdef double c00, c10, c01, c11
  92 + cdef double c0, c1
  93 + cdef double c
  94 +
  95 + cdef int x0 = <int>floor(x)
  96 + cdef int x1 = x0 + 1
  97 +
  98 + cdef int y0 = <int>floor(y)
  99 + cdef int y1 = y0 + 1
  100 +
  101 + cdef int z0 = <int>floor(z)
  102 + cdef int z1 = z0 + 1
  103 +
  104 + if x0 == x1:
  105 + xd = 1.0
  106 + else:
  107 + xd = (x - x0) / (x1 - x0)
  108 +
  109 + if y0 == y1:
  110 + yd = 1.0
  111 + else:
  112 + yd = (y - y0) / (y1 - y0)
  113 +
  114 + if z0 == z1:
  115 + zd = 1.0
  116 + else:
  117 + zd = (z - z0) / (z1 - z0)
  118 +
  119 + c00 = _G(V, x0, y0, z0)*(1 - xd) + _G(V, x1, y0, z0)*xd
  120 + c10 = _G(V, x0, y1, z0)*(1 - xd) + _G(V, x1, y1, z0)*xd
  121 + c01 = _G(V, x0, y0, z1)*(1 - xd) + _G(V, x1, y0, z1)*xd
  122 + c11 = _G(V, x0, y1, z1)*(1 - xd) + _G(V, x1, y1, z1)*xd
  123 +
  124 + c0 = c00*(1 - yd) + c10*yd
  125 + c1 = c01*(1 - yd) + c11*yd
  126 +
  127 + c = c0*(1 - zd) + c1*zd
  128 +
  129 + return c
  130 +
  131 +
  132 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  133 +@cython.cdivision(True)
  134 +@cython.wraparound(False)
  135 +cdef inline double lanczos3_L(double x, int a) nogil:
  136 + if x == 0:
  137 + return 1.0
  138 + elif -a <= x < a:
  139 + return (a * sin(M_PI * x) * sin(M_PI * (x / a)))/(M_PI**2 * x**2)
  140 + else:
  141 + return 0.0
  142 +
  143 +
  144 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  145 +@cython.cdivision(True)
  146 +@cython.wraparound(False)
  147 +cdef double lanczos3(image_t[:, :, :] V, double x, double y, double z) nogil:
  148 + cdef int a = LANCZOS_A
  149 +
  150 + cdef int xd = <int>floor(x)
  151 + cdef int yd = <int>floor(y)
  152 + cdef int zd = <int>floor(z)
  153 +
  154 + cdef int xi = xd - a + 1
  155 + cdef int xf = xd + a
  156 +
  157 + cdef int yi = yd - a + 1
  158 + cdef int yf = yd + a
  159 +
  160 + cdef int zi = zd - a + 1
  161 + cdef int zf = zd + a
  162 +
  163 + cdef double lx = 0.0
  164 + cdef double ly = 0.0
  165 + cdef double lz = 0.0
  166 +
  167 + cdef double[SIZE_LANCZOS_TMP][SIZE_LANCZOS_TMP] temp_x
  168 + cdef double[SIZE_LANCZOS_TMP] temp_y
  169 +
  170 + cdef int i, j, k
  171 + cdef int m, n, o
  172 +
  173 + m = 0
  174 + for k in xrange(zi, zf):
  175 + n = 0
  176 + for j in xrange(yi, yf):
  177 + lx = 0
  178 + for i in xrange(xi, xf):
  179 + lx += _G(V, i, j, k) * lanczos3_L(x - i, a)
  180 + temp_x[m][n] = lx
  181 + n += 1
  182 + m += 1
  183 +
  184 + m = 0
  185 + for k in xrange(zi, zf):
  186 + n = 0
  187 + ly = 0
  188 + for j in xrange(yi, yf):
  189 + ly += temp_x[m][n] * lanczos3_L(y - j, a)
  190 + n += 1
  191 + temp_y[m] = ly
  192 + m += 1
  193 +
  194 + m = 0
  195 + for k in xrange(zi, zf):
  196 + lz += temp_y[m] * lanczos3_L(z - k, a)
  197 + m += 1
  198 +
  199 + return lz
  200 +
  201 +
  202 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  203 +@cython.cdivision(True)
  204 +@cython.wraparound(False)
  205 +cdef image_t _G(image_t[:, :, :] V, int x, int y, int z) nogil:
  206 + cdef int dz, dy, dx
  207 + dz = V.shape[0] - 1
  208 + dy = V.shape[1] - 1
  209 + dx = V.shape[2] - 1
  210 +
  211 + if x < 0:
  212 + x = dx + x + 1
  213 + elif x > dx:
  214 + x = x - dx - 1
  215 +
  216 + if y < 0:
  217 + y = dy + y + 1
  218 + elif y > dy:
  219 + y = y - dy - 1
  220 +
  221 + if z < 0:
  222 + z = dz + z + 1
  223 + elif z > dz:
  224 + z = z - dz - 1
  225 +
  226 + return V[z, y, x]
  227 +
  228 +
  229 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  230 +@cython.cdivision(True)
  231 +@cython.wraparound(False)
  232 +cdef void calc_coef_tricub(image_t[:, :, :] V, double x, double y, double z, double [64] coef) nogil:
  233 + cdef int xi = <int>floor(x)
  234 + cdef int yi = <int>floor(y)
  235 + cdef int zi = <int>floor(z)
  236 +
  237 + cdef double[64] _x
  238 +
  239 + cdef int i, j
  240 +
  241 + _x[0] = _G(V, xi, yi, zi)
  242 + _x[1] = _G(V, xi + 1, yi, zi)
  243 + _x[2] = _G(V, xi, yi + 1, zi)
  244 + _x[3] = _G(V, xi + 1, yi + 1, zi)
  245 + _x[4] = _G(V, xi, yi, zi + 1)
  246 + _x[5] = _G(V, xi + 1, yi, zi + 1)
  247 + _x[6] = _G(V, xi, yi + 1, zi + 1)
  248 + _x[7] = _G(V, xi + 1, yi + 1, zi + 1)
  249 +
  250 + _x[8] = 0.5*(_G(V, xi+1,yi,zi) - _G(V, xi-1, yi, zi))
  251 + _x[9] = 0.5*(_G(V, xi+2,yi,zi) - _G(V, xi, yi, zi))
  252 + _x[10] = 0.5*(_G(V, xi+1, yi+1,zi) - _G(V, xi-1, yi+1, zi))
  253 + _x[11] = 0.5*(_G(V, xi+2, yi+1,zi) - _G(V, xi, yi+1, zi))
  254 + _x[12] = 0.5*(_G(V, xi+1, yi,zi+1) - _G(V, xi-1, yi, zi+1))
  255 + _x[13] = 0.5*(_G(V, xi+2, yi,zi+1) - _G(V, xi, yi, zi+1))
  256 + _x[14] = 0.5*(_G(V, xi+1, yi+1,zi+1) - _G(V, xi-1, yi+1, zi+1))
  257 + _x[15] = 0.5*(_G(V, xi+2, yi+1,zi+1) - _G(V, xi, yi+1, zi+1))
  258 + _x[16] = 0.5*(_G(V, xi, yi+1,zi) - _G(V, xi, yi-1, zi))
  259 + _x[17] = 0.5*(_G(V, xi+1, yi+1,zi) - _G(V, xi+1, yi-1, zi))
  260 + _x[18] = 0.5*(_G(V, xi, yi+2,zi) - _G(V, xi, yi, zi))
  261 + _x[19] = 0.5*(_G(V, xi+1, yi+2,zi) - _G(V, xi+1, yi, zi))
  262 + _x[20] = 0.5*(_G(V, xi, yi+1,zi+1) - _G(V, xi, yi-1, zi+1))
  263 + _x[21] = 0.5*(_G(V, xi+1, yi+1,zi+1) - _G(V, xi+1, yi-1, zi+1))
  264 + _x[22] = 0.5*(_G(V, xi, yi+2,zi+1) - _G(V, xi, yi, zi+1))
  265 + _x[23] = 0.5*(_G(V, xi+1, yi+2,zi+1) - _G(V, xi+1, yi, zi+1))
  266 + _x[24] = 0.5*(_G(V, xi, yi,zi+1) - _G(V, xi, yi, zi-1))
  267 + _x[25] = 0.5*(_G(V, xi+1, yi,zi+1) - _G(V, xi+1, yi, zi-1))
  268 + _x[26] = 0.5*(_G(V, xi, yi+1,zi+1) - _G(V, xi, yi+1, zi-1))
  269 + _x[27] = 0.5*(_G(V, xi+1, yi+1,zi+1) - _G(V, xi+1, yi+1, zi-1))
  270 + _x[28] = 0.5*(_G(V, xi, yi,zi+2) - _G(V, xi, yi, zi))
  271 + _x[29] = 0.5*(_G(V, xi+1, yi,zi+2) - _G(V, xi+1, yi, zi))
  272 + _x[30] = 0.5*(_G(V, xi, yi+1,zi+2) - _G(V, xi, yi+1, zi))
  273 + _x[31] = 0.5*(_G(V, xi+1, yi+1,zi+2) - _G(V, xi+1, yi+1, zi))
  274 +
  275 + _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))
  276 + _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))
  277 + _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))
  278 + _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))
  279 + _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))
  280 + _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))
  281 + _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))
  282 + _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))
  283 + _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))
  284 + _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))
  285 + _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))
  286 + _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))
  287 + _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))
  288 + _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))
  289 + _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))
  290 + _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))
  291 + _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))
  292 + _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))
  293 + _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))
  294 + _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))
  295 + _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))
  296 + _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))
  297 + _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))
  298 + _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))
  299 +
  300 + _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))
  301 + _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))
  302 + _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))
  303 + _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))
  304 + _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))
  305 + _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))
  306 + _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))
  307 + _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))
  308 +
  309 + for j in prange(64):
  310 + coef[j] = 0.0
  311 + for i in xrange(64):
  312 + coef[j] += (temp[j][i] * _x[i])
  313 +
  314 +
  315 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  316 +@cython.cdivision(True)
  317 +@cython.wraparound(False)
  318 +cdef double tricub_interpolate(image_t[:, :, :] V, double x, double y, double z) nogil:
  319 + # From: Tricubic interpolation in three dimensions. Lekien and Marsden
  320 + cdef double[64] coef
  321 + cdef double result = 0.0
  322 + calc_coef_tricub(V, x, y, z, coef)
  323 +
  324 + cdef int i, j, k
  325 +
  326 + cdef int xi = <int>floor(x)
  327 + cdef int yi = <int>floor(y)
  328 + cdef int zi = <int>floor(z)
  329 +
  330 + for i in xrange(4):
  331 + for j in xrange(4):
  332 + for k in xrange(4):
  333 + result += (coef[i+4*j+16*k] * ((x-xi)**i) * ((y-yi)**j) * ((z-zi)**k))
  334 + # return V[<int>z, <int>y, <int>x]
  335 + # with gil:
  336 + # print result
  337 + return result
  338 +
  339 +
  340 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  341 +@cython.cdivision(True)
  342 +@cython.wraparound(False)
  343 +cdef double cubicInterpolate(double p[4], double x) nogil:
  344 + 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])))
  345 +
  346 +
  347 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  348 +@cython.cdivision(True)
  349 +@cython.wraparound(False)
  350 +cdef double bicubicInterpolate (double p[4][4], double x, double y) nogil:
  351 + cdef double arr[4]
  352 + arr[0] = cubicInterpolate(p[0], y)
  353 + arr[1] = cubicInterpolate(p[1], y)
  354 + arr[2] = cubicInterpolate(p[2], y)
  355 + arr[3] = cubicInterpolate(p[3], y)
  356 + return cubicInterpolate(arr, x)
  357 +
  358 +
  359 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  360 +@cython.cdivision(True)
  361 +@cython.wraparound(False)
  362 +cdef double tricubicInterpolate(image_t[:, :, :] V, double x, double y, double z) nogil:
  363 + # From http://www.paulinternet.nl/?page=bicubic
  364 + cdef double p[4][4][4]
  365 +
  366 + cdef int xi = <int>floor(x)
  367 + cdef int yi = <int>floor(y)
  368 + cdef int zi = <int>floor(z)
  369 +
  370 + cdef int i, j, k
  371 +
  372 + for i in xrange(4):
  373 + for j in xrange(4):
  374 + for k in xrange(4):
  375 + p[i][j][k] = _G(V, xi + i -1, yi + j -1, zi + k - 1)
  376 +
  377 + cdef double arr[4]
  378 + arr[0] = bicubicInterpolate(p[0], y-yi, z-zi)
  379 + arr[1] = bicubicInterpolate(p[1], y-yi, z-zi)
  380 + arr[2] = bicubicInterpolate(p[2], y-yi, z-zi)
  381 + arr[3] = bicubicInterpolate(p[3], y-yi, z-zi)
  382 + return cubicInterpolate(arr, x-xi)
  383 +
  384 +
  385 +def tricub_interpolate_py(image_t[:, :, :] V, double x, double y, double z):
  386 + return tricub_interpolate(V, x, y, z)
  387 +
  388 +def tricub_interpolate2_py(image_t[:, :, :] V, double x, double y, double z):
  389 + return tricubicInterpolate(V, x, y, z)
  390 +
  391 +def trilin_interpolate_py(image_t[:, :, :] V, double x, double y, double z):
  392 + return interpolate(V, x, y, z)
invesalius_cy/mips.pyx 0 → 100644
  Diff comments   View file @18c2bcd
@@ -0,0 +1,379 @@ @@ -0,0 +1,379 @@
  1 +#http://en.wikipedia.org/wiki/Local_maximum_intensity_projection
  2 +import numpy as np
  3 +cimport numpy as np
  4 +cimport cython
  5 +
  6 +from libc.math cimport floor, ceil, sqrt, fabs
  7 +from cython.parallel import prange
  8 +
  9 +DTYPE = np.uint8
  10 +ctypedef np.uint8_t DTYPE_t
  11 +
  12 +DTYPE16 = np.int16
  13 +ctypedef np.int16_t DTYPE16_t
  14 +
  15 +DTYPEF32 = np.float32
  16 +ctypedef np.float32_t DTYPEF32_t
  17 +
  18 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  19 +def lmip(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t tmin,
  20 + DTYPE16_t tmax, np.ndarray[DTYPE16_t, ndim=2] out):
  21 + cdef DTYPE16_t max
  22 + cdef int start
  23 + cdef int sz = image.shape[0]
  24 + cdef int sy = image.shape[1]
  25 + cdef int sx = image.shape[2]
  26 +
  27 + # AXIAL
  28 + if axis == 0:
  29 + for x in xrange(sx):
  30 + for y in xrange(sy):
  31 + max = image[0, y, x]
  32 + if max >= tmin and max <= tmax:
  33 + start = 1
  34 + else:
  35 + start = 0
  36 + for z in xrange(sz):
  37 + if image[z, y, x] > max:
  38 + max = image[z, y, x]
  39 +
  40 + elif image[z, y, x] < max and start:
  41 + break
  42 +
  43 + if image[z, y, x] >= tmin and image[z, y, x] <= tmax:
  44 + start = 1
  45 +
  46 + out[y, x] = max
  47 +
  48 + #CORONAL
  49 + elif axis == 1:
  50 + for z in xrange(sz):
  51 + for x in xrange(sx):
  52 + max = image[z, 0, x]
  53 + if max >= tmin and max <= tmax:
  54 + start = 1
  55 + else:
  56 + start = 0
  57 + for y in xrange(sy):
  58 + if image[z, y, x] > max:
  59 + max = image[z, y, x]
  60 +
  61 + elif image[z, y, x] < max and start:
  62 + break
  63 +
  64 + if image[z, y, x] >= tmin and image[z, y, x] <= tmax:
  65 + start = 1
  66 +
  67 + out[z, x] = max
  68 +
  69 + #CORONAL
  70 + elif axis == 2:
  71 + for z in xrange(sz):
  72 + for y in xrange(sy):
  73 + max = image[z, y, 0]
  74 + if max >= tmin and max <= tmax:
  75 + start = 1
  76 + else:
  77 + start = 0
  78 + for x in xrange(sx):
  79 + if image[z, y, x] > max:
  80 + max = image[z, y, x]
  81 +
  82 + elif image[z, y, x] < max and start:
  83 + break
  84 +
  85 + if image[z, y, x] >= tmin and image[z, y, x] <= tmax:
  86 + start = 1
  87 +
  88 + out[z, y] = max
  89 +
  90 +
  91 +cdef DTYPE16_t get_colour(DTYPE16_t vl, DTYPE16_t wl, DTYPE16_t ww):
  92 + cdef DTYPE16_t out_colour
  93 + cdef DTYPE16_t min_value = wl - (ww / 2)
  94 + cdef DTYPE16_t max_value = wl + (ww / 2)
  95 + if vl < min_value:
  96 + out_colour = min_value
  97 + elif vl > max_value:
  98 + out_colour = max_value
  99 + else:
  100 + out_colour = vl
  101 +
  102 + return out_colour
  103 +
  104 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  105 +@cython.cdivision(True)
  106 +cdef float get_opacity(DTYPE16_t vl, DTYPE16_t wl, DTYPE16_t ww) nogil:
  107 + cdef float out_opacity
  108 + cdef DTYPE16_t min_value = wl - (ww / 2)
  109 + cdef DTYPE16_t max_value = wl + (ww / 2)
  110 + if vl < min_value:
  111 + out_opacity = 0.0
  112 + elif vl > max_value:
  113 + out_opacity = 1.0
  114 + else:
  115 + out_opacity = 1.0/(max_value - min_value) * (vl - min_value)
  116 +
  117 + return out_opacity
  118 +
  119 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  120 +@cython.cdivision(True)
  121 +cdef float get_opacity_f32(DTYPEF32_t vl, DTYPE16_t wl, DTYPE16_t ww) nogil:
  122 + cdef float out_opacity
  123 + cdef DTYPE16_t min_value = wl - (ww / 2)
  124 + cdef DTYPE16_t max_value = wl + (ww / 2)
  125 + if vl < min_value:
  126 + out_opacity = 0.0
  127 + elif vl > max_value:
  128 + out_opacity = 1.0
  129 + else:
  130 + out_opacity = 1.0/(max_value - min_value) * (vl - min_value)
  131 +
  132 + return out_opacity
  133 +
  134 +
  135 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  136 +@cython.cdivision(True)
  137 +def mida(np.ndarray[DTYPE16_t, ndim=3] image, int axis, DTYPE16_t wl,
  138 + DTYPE16_t ww, np.ndarray[DTYPE16_t, ndim=2] out):
  139 + cdef int sz = image.shape[0]
  140 + cdef int sy = image.shape[1]
  141 + cdef int sx = image.shape[2]
  142 +
  143 + cdef DTYPE16_t min = image.min()
  144 + cdef DTYPE16_t max = image.max()
  145 + cdef DTYPE16_t vl
  146 +
  147 + cdef DTYPE16_t min_value = wl - (ww / 2)
  148 + cdef DTYPE16_t max_value = wl + (ww / 2)
  149 +
  150 + cdef float fmax=0.0
  151 + cdef float fpi
  152 + cdef float dl
  153 + cdef float bt
  154 +
  155 + cdef float alpha
  156 + cdef float alpha_p = 0.0
  157 + cdef float colour
  158 + cdef float colour_p = 0
  159 +
  160 + cdef int x, y, z
  161 +
  162 + # AXIAL
  163 + if axis == 0:
  164 + for x in prange(sx, nogil=True):
  165 + for y in xrange(sy):
  166 + fmax = 0.0
  167 + alpha_p = 0.0
  168 + colour_p = 0.0
  169 + for z in xrange(sz):
  170 + vl = image[z, y, x]
  171 + fpi = 1.0/(max - min) * (vl - min)
  172 + if fpi > fmax:
  173 + dl = fpi - fmax
  174 + fmax = fpi
  175 + else:
  176 + dl = 0.0
  177 +
  178 + bt = 1.0 - dl
  179 +
  180 + colour = fpi
  181 + alpha = get_opacity(vl, wl, ww)
  182 + colour = (bt * colour_p) + (1 - bt * alpha_p) * colour * alpha
  183 + alpha = (bt * alpha_p) + (1 - bt * alpha_p) * alpha
  184 +
  185 + colour_p = colour
  186 + alpha_p = alpha
  187 +
  188 + if alpha >= 1.0:
  189 + break
  190 +
  191 +
  192 + #out[y, x] = <DTYPE16_t>((max_value - min_value) * colour + min_value)
  193 + out[y, x] = <DTYPE16_t>((max - min) * colour + min)
  194 +
  195 +
  196 + #CORONAL
  197 + elif axis == 1:
  198 + for z in prange(sz, nogil=True):
  199 + for x in xrange(sx):
  200 + fmax = 0.0
  201 + alpha_p = 0.0
  202 + colour_p = 0.0
  203 + for y in xrange(sy):
  204 + vl = image[z, y, x]
  205 + fpi = 1.0/(max - min) * (vl - min)
  206 + if fpi > fmax:
  207 + dl = fpi - fmax
  208 + fmax = fpi
  209 + else:
  210 + dl = 0.0
  211 +
  212 + bt = 1.0 - dl
  213 +
  214 + colour = fpi
  215 + alpha = get_opacity(vl, wl, ww)
  216 + colour = (bt * colour_p) + (1 - bt * alpha_p) * colour * alpha
  217 + alpha = (bt * alpha_p) + (1 - bt * alpha_p) * alpha
  218 +
  219 + colour_p = colour
  220 + alpha_p = alpha
  221 +
  222 + if alpha >= 1.0:
  223 + break
  224 +
  225 + out[z, x] = <DTYPE16_t>((max - min) * colour + min)
  226 +
  227 + #AXIAL
  228 + elif axis == 2:
  229 + for z in prange(sz, nogil=True):
  230 + for y in xrange(sy):
  231 + fmax = 0.0
  232 + alpha_p = 0.0
  233 + colour_p = 0.0
  234 + for x in xrange(sx):
  235 + vl = image[z, y, x]
  236 + fpi = 1.0/(max - min) * (vl - min)
  237 + if fpi > fmax:
  238 + dl = fpi - fmax
  239 + fmax = fpi
  240 + else:
  241 + dl = 0.0
  242 +
  243 + bt = 1.0 - dl
  244 +
  245 + colour = fpi
  246 + alpha = get_opacity(vl, wl, ww)
  247 + colour = (bt * colour_p) + (1 - bt * alpha_p) * colour * alpha
  248 + alpha = (bt * alpha_p) + (1 - bt * alpha_p) * alpha
  249 +
  250 + colour_p = colour
  251 + alpha_p = alpha
  252 +
  253 + if alpha >= 1.0:
  254 + break
  255 +
  256 + out[z, y] = <DTYPE16_t>((max - min) * colour + min)
  257 +
  258 +
  259 +
  260 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  261 +@cython.cdivision(True)
  262 +cdef inline void finite_difference(DTYPE16_t[:, :, :] image,
  263 + int x, int y, int z, float h, float *g) nogil:
  264 + cdef int px, py, pz, fx, fy, fz
  265 +
  266 + cdef int sz = image.shape[0]
  267 + cdef int sy = image.shape[1]
  268 + cdef int sx = image.shape[2]
  269 +
  270 + cdef float gx, gy, gz
  271 +
  272 + if x == 0:
  273 + px = 0
  274 + fx = 1
  275 + elif x == sx - 1:
  276 + px = x - 1
  277 + fx = x
  278 + else:
  279 + px = x - 1
  280 + fx = x + 1
  281 +
  282 + if y == 0:
  283 + py = 0
  284 + fy = 1
  285 + elif y == sy - 1:
  286 + py = y - 1
  287 + fy = y
  288 + else:
  289 + py = y - 1
  290 + fy = y + 1
  291 +
  292 + if z == 0:
  293 + pz = 0
  294 + fz = 1
  295 + elif z == sz - 1:
  296 + pz = z - 1
  297 + fz = z
  298 + else:
  299 + pz = z - 1
  300 + fz = z + 1
  301 +
  302 + gx = (image[z, y, fx] - image[z, y, px]) / (2*h)
  303 + gy = (image[z, fy, x] - image[z, py, x]) / (2*h)
  304 + gz = (image[fz, y, x] - image[pz, y, x]) / (2*h)
  305 +
  306 + g[0] = gx
  307 + g[1] = gy
  308 + g[2] = gz
  309 +
  310 +
  311 +
  312 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  313 +@cython.cdivision(True)
  314 +cdef inline float calc_fcm_itensity(DTYPE16_t[:, :, :] image,
  315 + int x, int y, int z, float n, float* dir) nogil:
  316 + cdef float g[3]
  317 + finite_difference(image, x, y, z, 1.0, g)
  318 + cdef float gm = sqrt(g[0]*g[0] + g[1]*g[1] + g[2]*g[2])
  319 + cdef float d = g[0]*dir[0] + g[1]*dir[1] + g[2]*dir[2]
  320 + cdef float sf = (1.0 - fabs(d/gm))**n
  321 + #alpha = get_opacity_f32(gm, wl, ww)
  322 + cdef float vl = gm * sf
  323 + return vl
  324 +
  325 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  326 +@cython.cdivision(True)
  327 +def fast_countour_mip(np.ndarray[DTYPE16_t, ndim=3] image,
  328 + float n,
  329 + int axis,
  330 + DTYPE16_t wl, DTYPE16_t ww,
  331 + int tmip,
  332 + np.ndarray[DTYPE16_t, ndim=2] out):
  333 + cdef int sz = image.shape[0]
  334 + cdef int sy = image.shape[1]
  335 + cdef int sx = image.shape[2]
  336 + cdef float gm
  337 + cdef float alpha
  338 + cdef float sf
  339 + cdef float d
  340 +
  341 + cdef float* g
  342 + cdef float* dir = [ 0, 0, 0 ]
  343 +
  344 + cdef DTYPE16_t[:, :, :] vimage = image
  345 + cdef np.ndarray[DTYPE16_t, ndim=3] tmp = np.empty_like(image)
  346 +
  347 + cdef DTYPE16_t min = image.min()
  348 + cdef DTYPE16_t max = image.max()
  349 + cdef float fmin = <float>min
  350 + cdef float fmax = <float>max
  351 + cdef float vl
  352 + cdef DTYPE16_t V
  353 +
  354 + cdef int x, y, z
  355 +
  356 + if axis == 0:
  357 + dir[2] = 1.0
  358 + elif axis == 1:
  359 + dir[1] = 1.0
  360 + elif axis == 2:
  361 + dir[0] = 1.0
  362 +
  363 + for z in prange(sz, nogil=True):
  364 + for y in range(sy):
  365 + for x in range(sx):
  366 + vl = calc_fcm_itensity(vimage, x, y, z, n, dir)
  367 + tmp[z, y, x] = <DTYPE16_t>vl
  368 +
  369 + cdef DTYPE16_t tmin = tmp.min()
  370 + cdef DTYPE16_t tmax = tmp.max()
  371 +
  372 + #tmp = ((max - min)/<float>(tmax - tmin)) * (tmp - tmin) + min
  373 +
  374 + if tmip == 0:
  375 + out[:] = tmp.max(axis)
  376 + elif tmip == 1:
  377 + lmip(tmp, axis, 700, 3033, out)
  378 + elif tmip == 2:
  379 + mida(tmp, axis, wl, ww, out)
invesalius_cy/transforms.pyx 0 → 100644
  Diff comments   View file @18c2bcd
@@ -0,0 +1,174 @@ @@ -0,0 +1,174 @@
  1 +import numpy as np
  2 +cimport numpy as np
  3 +cimport cython
  4 +
  5 +from .cy_my_types cimport image_t
  6 +from .interpolation cimport interpolate, tricub_interpolate, tricubicInterpolate, lanczos3, nearest_neighbour_interp
  7 +
  8 +from libc.math cimport floor, ceil, sqrt, fabs, round
  9 +from cython.parallel import prange
  10 +
  11 +ctypedef double (*interp_function)(image_t[:, :, :], double, double, double) nogil
  12 +
  13 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  14 +@cython.cdivision(True)
  15 +@cython.wraparound(False)
  16 +cdef void mul_mat4_vec4(double[:, :] M,
  17 + double* coord,
  18 + double* out) nogil:
  19 +
  20 + out[0] = coord[0] * M[0, 0] + coord[1] * M[0, 1] + coord[2] * M[0, 2] + coord[3] * M[0, 3]
  21 + out[1] = coord[0] * M[1, 0] + coord[1] * M[1, 1] + coord[2] * M[1, 2] + coord[3] * M[1, 3]
  22 + out[2] = coord[0] * M[2, 0] + coord[1] * M[2, 1] + coord[2] * M[2, 2] + coord[3] * M[2, 3]
  23 + out[3] = coord[0] * M[3, 0] + coord[1] * M[3, 1] + coord[2] * M[3, 2] + coord[3] * M[3, 3]
  24 +
  25 +
  26 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  27 +@cython.cdivision(True)
  28 +@cython.wraparound(False)
  29 +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:
  30 +
  31 + cdef double coord[4]
  32 + coord[0] = z*sz
  33 + coord[1] = y*sy
  34 + coord[2] = x*sx
  35 + coord[3] = 1.0
  36 +
  37 + cdef double _ncoord[4]
  38 + _ncoord[3] = 1
  39 + # _ncoord[:] = [0.0, 0.0, 0.0, 1.0]
  40 +
  41 + cdef unsigned int dz, dy, dx
  42 + dz = volume.shape[0]
  43 + dy = volume.shape[1]
  44 + dx = volume.shape[2]
  45 +
  46 +
  47 + mul_mat4_vec4(M, coord, _ncoord)
  48 +
  49 + cdef double nz, ny, nx
  50 + nz = (_ncoord[0]/_ncoord[3])/sz
  51 + ny = (_ncoord[1]/_ncoord[3])/sy
  52 + nx = (_ncoord[2]/_ncoord[3])/sx
  53 +
  54 + cdef double v
  55 +
  56 + if 0 <= nz <= (dz-1) and 0 <= ny <= (dy-1) and 0 <= nx <= (dx-1):
  57 + return <image_t>f_interp(volume, nx, ny, nz)
  58 + # if minterpol == 0:
  59 + # return <image_t>nearest_neighbour_interp(volume, nx, ny, nz)
  60 + # elif minterpol == 1:
  61 + # return <image_t>interpolate(volume, nx, ny, nz)
  62 + # elif minterpol == 2:
  63 + # v = tricubicInterpolate(volume, nx, ny, nz)
  64 + # if (v < cval):
  65 + # v = cval
  66 + # return <image_t>v
  67 + # else:
  68 + # v = lanczos3(volume, nx, ny, nz)
  69 + # if (v < cval):
  70 + # v = cval
  71 + # return <image_t>v
  72 + else:
  73 + return cval
  74 +
  75 +
  76 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  77 +@cython.cdivision(True)
  78 +@cython.wraparound(False)
  79 +def apply_view_matrix_transform(image_t[:, :, :] volume,
  80 + spacing,
  81 + double[:, :] M,
  82 + unsigned int n, str orientation,
  83 + int minterpol,
  84 + image_t cval,
  85 + image_t[:, :, :] out):
  86 +
  87 + cdef int dz, dy, dx
  88 + cdef int z, y, x
  89 + dz = volume.shape[0]
  90 + dy = volume.shape[1]
  91 + dx = volume.shape[2]
  92 +
  93 + cdef unsigned int odz, ody, odx
  94 + odz = out.shape[0]
  95 + ody = out.shape[1]
  96 + odx = out.shape[2]
  97 +
  98 + cdef unsigned int count = 0
  99 +
  100 + cdef double sx, sy, sz
  101 + sx = spacing[0]
  102 + sy = spacing[1]
  103 + sz = spacing[2]
  104 +
  105 + cdef interp_function f_interp
  106 +
  107 + if minterpol == 0:
  108 + f_interp = nearest_neighbour_interp
  109 + elif minterpol == 1:
  110 + f_interp = interpolate
  111 + elif minterpol == 2:
  112 + f_interp = tricubicInterpolate
  113 + else:
  114 + f_interp = lanczos3
  115 +
  116 + if orientation == 'AXIAL':
  117 + for z in xrange(n, n+odz):
  118 + for y in prange(dy, nogil=True):
  119 + for x in xrange(dx):
  120 + out[count, y, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval)
  121 + count += 1
  122 +
  123 + elif orientation == 'CORONAL':
  124 + for y in xrange(n, n+ody):
  125 + for z in prange(dz, nogil=True):
  126 + for x in xrange(dx):
  127 + out[z, count, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval)
  128 + count += 1
  129 +
  130 + elif orientation == 'SAGITAL':
  131 + for x in xrange(n, n+odx):
  132 + for z in prange(dz, nogil=True):
  133 + for y in xrange(dy):
  134 + out[z, y, count] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval)
  135 + count += 1
  136 +
  137 +
  138 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  139 +@cython.cdivision(True)
  140 +@cython.wraparound(False)
  141 +def convolve_non_zero(image_t[:, :, :] volume,
  142 + image_t[:, :, :] kernel,
  143 + image_t cval):
  144 + cdef Py_ssize_t x, y, z, sx, sy, sz, kx, ky, kz, skx, sky, skz, i, j, k
  145 + cdef image_t v
  146 +
  147 + cdef image_t[:, :, :] out = np.zeros_like(volume)
  148 +
  149 + sz = volume.shape[0]
  150 + sy = volume.shape[1]
  151 + sx = volume.shape[2]
  152 +
  153 + skz = kernel.shape[0]
  154 + sky = kernel.shape[1]
  155 + skx = kernel.shape[2]
  156 +
  157 + for z in prange(sz, nogil=True):
  158 + for y in xrange(sy):
  159 + for x in xrange(sx):
  160 + if volume[z, y, x] != 0:
  161 + for k in xrange(skz):
  162 + kz = z - skz // 2 + k
  163 + for j in xrange(sky):
  164 + ky = y - sky // 2 + j
  165 + for i in xrange(skx):
  166 + kx = x - skx // 2 + i
  167 +
  168 + if 0 <= kz < sz and 0 <= ky < sy and 0 <= kx < sx:
  169 + v = volume[kz, ky, kx]
  170 + else:
  171 + v = cval
  172 +
  173 + out[z, y, x] += (v * kernel[k, j, i])
  174 + return np.asarray(out)
@@ -40,25 +40,25 @@ setup( @@ -40,25 +40,25 @@ setup(
40 ext_modules=cythonize( 40 ext_modules=cythonize(
41 [ 41 [
42 Extension( 42 Extension(
43 - "invesalius.data.mips",  
44 - ["invesalius/data/mips.pyx"], 43 + "invesalius_cy.mips",
  44 + ["invesalius_cy/mips.pyx"],
45 ), 45 ),
46 Extension( 46 Extension(
47 - "invesalius.data.interpolation",  
48 - ["invesalius/data/interpolation.pyx"], 47 + "invesalius_cy.interpolation",
  48 + ["invesalius_cy/interpolation.pyx"],
49 ), 49 ),
50 Extension( 50 Extension(
51 - "invesalius.data.transforms",  
52 - ["invesalius/data/transforms.pyx"], 51 + "invesalius_cy.transforms",
  52 + ["invesalius_cy/transforms.pyx"],
53 ), 53 ),
54 Extension( 54 Extension(
55 - "invesalius.data.floodfill",  
56 - ["invesalius/data/floodfill.pyx"], 55 + "invesalius_cy.floodfill",
  56 + ["invesalius_cy/floodfill.pyx"],
57 language="c++", 57 language="c++",
58 ), 58 ),
59 Extension( 59 Extension(
60 - "invesalius.data.cy_mesh",  
61 - ["invesalius/data/cy_mesh.pyx"], 60 + "invesalius_cy.cy_mesh",
  61 + ["invesalius_cy/cy_mesh.pyx"],
62 language="c++", 62 language="c++",
63 ), 63 ),
64 ] 64 ]