Download as
Browse Code »

Commit 8bdc2a60fedf36c4fa80f5cf4f49dd0d827c5b62

Authored by Thiago Franco de Moraes
1 parent d0d68211
Exists in master and in 3 other branches inv3.0.1, usp-navegador, wxgtk3_b5

Adding mips.pyx to invesalius project

Showing 2 changed files with 392 additions and 0 deletions   Show diff stats
invesalius/data/mips.pyx 0 → 100644
  Diff comments   View file @8bdc2a6
@@ -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)
setup.py 0 → 100644
  Diff comments   View file @8bdc2a6
@@ -0,0 +1,13 @@ @@ -0,0 +1,13 @@
  1 +from distutils.core import setup
  2 +from distutils.extension import Extension
  3 +from Cython.Distutils import build_ext
  4 +
  5 +import numpy
  6 +
  7 +setup(
  8 + cmdclass = {'build_ext': build_ext},
  9 + ext_modules = [ Extension("invesalius/data/mips", ["invesalius/data/mips.pyx"],
  10 + include_dirs = [numpy.get_include()],
  11 + extra_compile_args=['-fopenmp'],
  12 + extra_link_args=['-fopenmp'],)]
  13 + )