Download as
Browse Code »

Commit 6716954c39efab9846c28e8562f28e5b6275e315

Authored by Thiago Franco de Moraes
1 parent 37109ae4
Exists in smoothing_bin_pll

Added smooth_cy (only the algorithm) to InVesalius

Showing 2 changed files with 203 additions and 9 deletions   Show diff stats
invesalius/data/smooth_cy.pyx 0 → 100644
  Diff comments   View file @6716954
@@ -0,0 +1,182 @@ @@ -0,0 +1,182 @@
  1 +import numpy as np
  2 +cimport numpy as np
  3 +cimport cython
  4 +
  5 +from libc.math cimport floor, ceil, sqrt, fabs, round
  6 +from cython.parallel import prange
  7 +
  8 +DTYPE8 = np.uint8
  9 +ctypedef np.uint8_t DTYPE8_t
  10 +
  11 +DTYPEF64 = np.float64
  12 +ctypedef np.float64_t DTYPEF64_t
  13 +
  14 +
  15 +
  16 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  17 +@cython.cdivision(True)
  18 +cdef inline DTYPEF64_t GS(DTYPEF64_t[:, :, :] I, int z, int y, int x) nogil:
  19 + cdef int dz = I.shape[0]
  20 + cdef int dy = I.shape[1]
  21 + cdef int dx = I.shape[2]
  22 +
  23 + if 0 <= x < dx \
  24 + and 0 <= y < dy \
  25 + and 0 <= z < dz:
  26 + return I[z, y, x]
  27 + else:
  28 + return 0
  29 +
  30 +
  31 +
  32 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  33 +@cython.cdivision(True)
  34 +cdef void perim(DTYPE8_t[:, :, :] image,
  35 + DTYPE8_t[:, :, :] out) nogil:
  36 +
  37 + cdef int dz = image.shape[0]
  38 + cdef int dy = image.shape[1]
  39 + cdef int dx = image.shape[2]
  40 +
  41 + cdef int z, y, x
  42 + cdef int z_, y_, x_
  43 +
  44 + for z in prange(dz, nogil=True):
  45 + for y in xrange(dy):
  46 + for x in xrange(dx):
  47 + for z_ in xrange(z-1, z+2, 2):
  48 + for y_ in xrange(y-1, y+2, 2):
  49 + for x_ in xrange(x-1, x+2, 2):
  50 + if 0 <= x_ < dx \
  51 + and 0 <= y_ < dy \
  52 + and 0 <= z_ < dz \
  53 + and image[z, y, x] != image[z_, y_, x_]:
  54 + out[z, y, x] = 1
  55 + break
  56 +
  57 +
  58 +
  59 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  60 +@cython.cdivision(True)
  61 +cdef DTYPEF64_t calculate_H(DTYPEF64_t[:, :, :] I, int z, int y, int x) nogil:
  62 + # double fx, fy, fz, fxx, fyy, fzz, fxy, fxz, fyz, H
  63 + cdef DTYPEF64_t fx, fy, fz, fxx, fyy, fzz, fxy, fxz, fyz, H
  64 + # int h, k, l
  65 + cdef int h = 1
  66 + cdef int k = 1
  67 + cdef int l = 1
  68 +
  69 + fx = (GS(I, z, y, x + h) - GS(I, z, y, x - h)) / (2.0*h)
  70 + fy = (GS(I, z, y + k, x) - GS(I, z, y - k, x)) / (2.0*k)
  71 + fz = (GS(I, z + l, y, x) - GS(I, z - l, y, x)) / (2.0*l)
  72 +
  73 + fxx = (GS(I, z, y, x + h) - 2*GS(I, z, y, x) + GS(I, z, y, x - h)) / (h*h)
  74 + fyy = (GS(I, z, y + k, x) - 2*GS(I, z, y, x) + GS(I, z, y - k, x)) / (k*k)
  75 + fzz = (GS(I, z + l, y, x) - 2*GS(I, z, y, x) + GS(I, z - l, y, x)) / (l*l)
  76 +
  77 + fxy = (GS(I, z, y + k, x + h) - GS(I, z, y - k, x + h) \
  78 + - GS(I, z, y + k, x - h) + GS(I, z, y - k, x - h)) \
  79 + / (4.0*h*k)
  80 + fxz = (GS(I, z + l, y, x + h) - GS(I, z + l, y, x - h) \
  81 + - GS(I, z - l, y, x + h) + GS(I, z - l, y, x - h)) \
  82 + / (4.0*h*l)
  83 + fyz = (GS(I, z + l, y + k, x) - GS(I, z + l, y - k, x) \
  84 + - GS(I, z - l, y + k, x) + GS(I, z - l, y - k, x)) \
  85 + / (4.0*k*l)
  86 +
  87 + if fx*fx + fy*fy + fz*fz > 0:
  88 + H = ((fy*fy + fz*fz)*fxx + (fx*fx + fz*fz)*fyy \
  89 + + (fx*fx + fy*fy)*fzz - 2*(fx*fy*fxy \
  90 + + fx*fz*fxz + fy*fz*fyz)) \
  91 + / (fx*fx + fy*fy + fz*fz)
  92 + else:
  93 + H = 0.0
  94 +
  95 + return H
  96 +
  97 +
  98 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  99 +@cython.cdivision(True)
  100 +cdef void replicate(DTYPEF64_t[:, :, :] source, DTYPEF64_t[:, :, :] dest) nogil:
  101 + cdef int dz = source.shape[0]
  102 + cdef int dy = source.shape[1]
  103 + cdef int dx = source.shape[2]
  104 + cdef int x, y, z
  105 + for z in prange(dz, nogil=True):
  106 + for y in xrange(dy):
  107 + for x in xrange(dx):
  108 + dest[z, y, x] = source[z, y, x]
  109 +
  110 +
  111 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  112 +@cython.cdivision(True)
  113 +def smooth(np.ndarray[DTYPE8_t, ndim=3] image,
  114 + int n, int bsize,
  115 + np.ndarray[DTYPEF64_t, ndim=3] out):
  116 +
  117 + cdef np.ndarray[DTYPE8_t, ndim=3] mask = np.zeros_like(image)
  118 + cdef np.ndarray[DTYPE8_t, ndim=3] _mask = np.zeros_like(image)
  119 + cdef np.ndarray[DTYPEF64_t, ndim=3] aux = np.zeros_like(out)
  120 +
  121 + cdef int i, x, y, z, S;
  122 + cdef DTYPEF64_t H, v, cn
  123 + cdef DTYPEF64_t diff=0.0
  124 + cdef DTYPEF64_t dt=1/6.0
  125 +
  126 + cdef DTYPEF64_t E = 0.001
  127 +
  128 + _mask[:] = image
  129 + for i in xrange(bsize):
  130 + perim(_mask, mask)
  131 + _mask[:] = mask
  132 + print i
  133 +
  134 + # out[:] = mask
  135 +
  136 + del _mask
  137 +
  138 + cdef int dz = image.shape[0]
  139 + cdef int dy = image.shape[1]
  140 + cdef int dx = image.shape[2]
  141 +
  142 + S = 0
  143 + for z in prange(dz, nogil=True):
  144 + for y in xrange(dy):
  145 + for x in xrange(dx):
  146 + if image[z, y, x]:
  147 + out[z, y, x] = 1.0
  148 + else:
  149 + out[z, y, x] = -1.0
  150 +
  151 + if mask[z, y, x]:
  152 + S += 1
  153 +
  154 + for i in xrange(n):
  155 + replicate(out, aux);
  156 + diff = 0.0;
  157 +
  158 + for z in xrange(dz):
  159 + for y in xrange(dy):
  160 + for x in xrange(dx):
  161 + if mask[z, y, x]:
  162 + H = calculate_H(aux, z, y, x);
  163 + v = aux[z, y, x] + dt*H;
  164 +
  165 + if image[z, y, x]:
  166 + if v > 0:
  167 + out[z, y, x] = v
  168 + else:
  169 + out[z, y, x] = 0.0001
  170 + else:
  171 + if v < 0:
  172 + out[z, y, x] = v
  173 + else:
  174 + out[z, y, x] = -0.0001
  175 +
  176 + diff += (out[z, y, x] - aux[z, y, x])*(out[z, y, x] - aux[z, y, x])
  177 +
  178 + cn = sqrt((1.0/S) * diff);
  179 + print "%d - CN: %.28f - diff: %.28f\n" % (i, cn, diff)
  180 +
  181 + if cn <= E:
  182 + break;
@@ -12,14 +12,14 @@ if sys.platform == &#39;linux2&#39;: @@ -12,14 +12,14 @@ if sys.platform == &#39;linux2&#39;:
12 setup( 12 setup(
13 cmdclass = {'build_ext': build_ext}, 13 cmdclass = {'build_ext': build_ext},
14 ext_modules = cythonize([ Extension("invesalius.data.mips", ["invesalius/data/mips.pyx"], 14 ext_modules = cythonize([ Extension("invesalius.data.mips", ["invesalius/data/mips.pyx"],
15 - include_dirs = [numpy.get_include()],  
16 - extra_compile_args=['-fopenmp'],  
17 - extra_link_args=['-fopenmp']), 15 + include_dirs = [numpy.get_include()],
  16 + extra_compile_args=['-fopenmp'],
  17 + extra_link_args=['-fopenmp']),
18 18
19 - Extension("invesalius.data.interpolation", ["invesalius/data/interpolation.pyx"],  
20 - include_dirs=[numpy.get_include()],  
21 - extra_compile_args=['-fopenmp',],  
22 - extra_link_args=['-fopenmp',]), 19 + Extension("invesalius.data.interpolation", ["invesalius/data/interpolation.pyx"],
  20 + include_dirs=[numpy.get_include()],
  21 + extra_compile_args=['-fopenmp',],
  22 + extra_link_args=['-fopenmp',]),
23 23
24 Extension("invesalius.data.transforms", ["invesalius/data/transforms.pyx"], 24 Extension("invesalius.data.transforms", ["invesalius/data/transforms.pyx"],
25 include_dirs=[numpy.get_include()], 25 include_dirs=[numpy.get_include()],
@@ -29,8 +29,20 @@ if sys.platform == &#39;linux2&#39;: @@ -29,8 +29,20 @@ if sys.platform == &#39;linux2&#39;:
29 Extension("invesalius.data.floodfill", ["invesalius/data/floodfill.pyx"], 29 Extension("invesalius.data.floodfill", ["invesalius/data/floodfill.pyx"],
30 include_dirs=[numpy.get_include()], 30 include_dirs=[numpy.get_include()],
31 language='c++',), 31 language='c++',),
32 - ])  
33 - ) 32 + Extension("invesalius.data.transforms", ["invesalius/data/transforms.pyx"],
  33 + include_dirs=[numpy.get_include()],
  34 + extra_compile_args=['-fopenmp',],
  35 + extra_link_args=['-fopenmp',]),
  36 +
  37 + # "Reducing Aliasing Artifacts in Iso-Surfaces of Binary Volumes"
  38 + Extension("invesalius.data.smooth_cy",
  39 + ["invesalius/data/smooth_cy.pyx"],
  40 + include_dirs=[numpy.get_include()],
  41 + extra_compile_args=['-fopenmp'],
  42 + extra_link_args=['-fopenmp'],),
  43 +
  44 + ])
  45 + )
34 46
35 elif sys.platform == 'win32': 47 elif sys.platform == 'win32':
36 setup( 48 setup(