Download as
Browse Code »

Commit 1f66dad197c97f63b733001eacd6ef789ddba1c4

Authored by Thiago Franco de Moraes
1 parent ae0a6c50
Exists in smoothing_bin_pll

Using memcpy and openmp

Showing 1 changed file with 83 additions and 32 deletions   Show diff stats
invesalius/data/smooth_cy.pyx
  Diff comments   View file @1f66dad
@@ -3,6 +3,7 @@ cimport numpy as np @@ -3,6 +3,7 @@ cimport numpy as np
3 cimport cython 3 cimport cython
4 4
5 from libc.math cimport floor, ceil, sqrt, fabs, round 5 from libc.math cimport floor, ceil, sqrt, fabs, round
  6 +from libc.string cimport memcpy
6 from cython.parallel import prange 7 from cython.parallel import prange
7 8
8 DTYPE8 = np.uint8 9 DTYPE8 = np.uint8
@@ -15,6 +16,7 @@ ctypedef np.float64_t DTYPEF64_t @@ -15,6 +16,7 @@ ctypedef np.float64_t DTYPEF64_t
15 16
16 @cython.boundscheck(False) # turn of bounds-checking for entire function 17 @cython.boundscheck(False) # turn of bounds-checking for entire function
17 @cython.cdivision(True) 18 @cython.cdivision(True)
  19 +@cython.wraparound(False)
18 cdef inline DTYPEF64_t GS(DTYPEF64_t[:, :, :] I, int z, int y, int x) nogil: 20 cdef inline DTYPEF64_t GS(DTYPEF64_t[:, :, :] I, int z, int y, int x) nogil:
19 cdef int dz = I.shape[0] 21 cdef int dz = I.shape[0]
20 cdef int dy = I.shape[1] 22 cdef int dy = I.shape[1]
@@ -31,6 +33,7 @@ cdef inline DTYPEF64_t GS(DTYPEF64_t[:, :, :] I, int z, int y, int x) nogil: @@ -31,6 +33,7 @@ cdef inline DTYPEF64_t GS(DTYPEF64_t[:, :, :] I, int z, int y, int x) nogil:
31 33
32 @cython.boundscheck(False) # turn of bounds-checking for entire function 34 @cython.boundscheck(False) # turn of bounds-checking for entire function
33 @cython.cdivision(True) 35 @cython.cdivision(True)
  36 +@cython.wraparound(False)
34 cdef void perim(DTYPE8_t[:, :, :] image, 37 cdef void perim(DTYPE8_t[:, :, :] image,
35 DTYPE8_t[:, :, :] out) nogil: 38 DTYPE8_t[:, :, :] out) nogil:
36 39
@@ -58,6 +61,7 @@ cdef void perim(DTYPE8_t[:, :, :] image, @@ -58,6 +61,7 @@ cdef void perim(DTYPE8_t[:, :, :] image,
58 61
59 @cython.boundscheck(False) # turn of bounds-checking for entire function 62 @cython.boundscheck(False) # turn of bounds-checking for entire function
60 @cython.cdivision(True) 63 @cython.cdivision(True)
  64 +@cython.wraparound(False)
61 cdef DTYPEF64_t calculate_H(DTYPEF64_t[:, :, :] I, int z, int y, int x) nogil: 65 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 66 # 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 67 cdef DTYPEF64_t fx, fy, fz, fxx, fyy, fzz, fxy, fxz, fyz, H
@@ -101,6 +105,7 @@ cdef DTYPEF64_t calculate_H(DTYPEF64_t[:, :, :] I, int z, int y, int x) nogil: @@ -101,6 +105,7 @@ cdef DTYPEF64_t calculate_H(DTYPEF64_t[:, :, :] I, int z, int y, int x) nogil:
101 105
102 @cython.boundscheck(False) # turn of bounds-checking for entire function 106 @cython.boundscheck(False) # turn of bounds-checking for entire function
103 @cython.cdivision(True) 107 @cython.cdivision(True)
  108 +@cython.wraparound(False)
104 cdef void replicate(DTYPEF64_t[:, :, :] source, DTYPEF64_t[:, :, :] dest) nogil: 109 cdef void replicate(DTYPEF64_t[:, :, :] source, DTYPEF64_t[:, :, :] dest) nogil:
105 cdef int dz = source.shape[0] 110 cdef int dz = source.shape[0]
106 cdef int dy = source.shape[1] 111 cdef int dy = source.shape[1]
@@ -111,28 +116,71 @@ cdef void replicate(DTYPEF64_t[:, :, :] source, DTYPEF64_t[:, :, :] dest) nogil: @@ -111,28 +116,71 @@ cdef void replicate(DTYPEF64_t[:, :, :] source, DTYPEF64_t[:, :, :] dest) nogil:
111 for x in xrange(dx): 116 for x in xrange(dx):
112 dest[z, y, x] = source[z, y, x] 117 dest[z, y, x] = source[z, y, x]
113 118
  119 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  120 +@cython.cdivision(True)
  121 +@cython.wraparound(False)
  122 +cdef void replicate8(DTYPE8_t[:, :, :] source, DTYPE8_t[:, :, :] dest) nogil:
  123 + cdef int dz = source.shape[0]
  124 + cdef int dy = source.shape[1]
  125 + cdef int dx = source.shape[2]
  126 + cdef int x, y, z
  127 + for z in prange(dz, nogil=True):
  128 + for y in xrange(dy):
  129 + for x in xrange(dx):
  130 + dest[z, y, x] = source[z, y, x]
  131 +
  132 +
  133 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  134 +@cython.cdivision(True)
  135 +@cython.wraparound(False)
  136 +cdef void _smooth(DTYPE8_t[:, :, :] image, DTYPEF64_t[:, :, :] aux, DTYPE8_t[:, :, :] mask, int x, int y, int z, DTYPEF64_t[:, :, :] out) nogil:
  137 +
  138 + cdef DTYPEF64_t H, v, cn
  139 + cdef DTYPEF64_t dt=1/6.0
  140 + H = calculate_H(aux, z, y, x)
  141 + v = aux[z, y, x] + dt*H
  142 +
  143 + if image[z, y, x]:
  144 + if v < 0:
  145 + out[z, y, x] = 0.00001
  146 + else:
  147 + out[z, y, x] = v
  148 + else:
  149 + if v > 0:
  150 + out[z, y, x] = -0.00001
  151 + else:
  152 + out[z, y, x] = v
  153 +
114 154
115 @cython.boundscheck(False) # turn of bounds-checking for entire function 155 @cython.boundscheck(False) # turn of bounds-checking for entire function
116 @cython.cdivision(True) 156 @cython.cdivision(True)
117 -def smooth(np.ndarray[DTYPE8_t, ndim=3] image, 157 +@cython.wraparound(False)
  158 +def smooth(DTYPE8_t[:, :, :] image,
118 int n, int bsize, 159 int n, int bsize,
119 - np.ndarray[DTYPEF64_t, ndim=3] out): 160 + DTYPEF64_t[:, :, :] out):
120 161
121 - cdef np.ndarray[DTYPE8_t, ndim=3] mask = np.zeros_like(image)  
122 - cdef np.ndarray[DTYPE8_t, ndim=3] _mask = np.zeros_like(image)  
123 - cdef np.ndarray[DTYPEF64_t, ndim=3] aux = np.zeros_like(out) 162 + cdef DTYPE8_t[:, :, :] mask = np.zeros_like(image)
  163 + cdef DTYPE8_t[:, :, :] _mask = np.zeros_like(image)
  164 + cdef DTYPEF64_t[:, :, :] aux = np.zeros_like(out)
124 165
125 cdef int i, x, y, z, S 166 cdef int i, x, y, z, S
126 cdef DTYPEF64_t H, v, cn 167 cdef DTYPEF64_t H, v, cn
127 cdef DTYPEF64_t diff=0.0 168 cdef DTYPEF64_t diff=0.0
128 cdef DTYPEF64_t dt=1/6.0 169 cdef DTYPEF64_t dt=1/6.0
129 170
  171 +
130 cdef DTYPEF64_t E = 0.001 172 cdef DTYPEF64_t E = 0.001
131 173
132 - _mask[:] = image 174 + print ">>>>>>>>>", image.size
  175 +
  176 + # _mask[:] = image
  177 + # replicate8(image, _mask)
  178 + memcpy(&_mask[0, 0, 0], &image[0, 0, 0], image.nbytes)
133 for i in xrange(bsize): 179 for i in xrange(bsize):
134 perim(_mask, mask) 180 perim(_mask, mask)
135 - _mask[:] = mask 181 + # _mask[:] = mask
  182 + # replicate8(mask, _mask)
  183 + memcpy(&_mask[0, 0, 0], &mask[0, 0, 0], mask.nbytes)
136 print i 184 print i
137 185
138 # out[:] = mask 186 # out[:] = mask
@@ -158,33 +206,36 @@ def smooth(np.ndarray[DTYPE8_t, ndim=3] image, @@ -158,33 +206,36 @@ def smooth(np.ndarray[DTYPE8_t, ndim=3] image,
158 S += 1 206 S += 1
159 207
160 for i in xrange(n): 208 for i in xrange(n):
161 - replicate(out, aux) 209 + # replicate(out, aux)
  210 + memcpy(&aux[0, 0, 0], &out[0, 0, 0], out.nbytes)
162 diff = 0.0 211 diff = 0.0
163 212
164 - for z in xrange(dz): 213 + for z in prange(dz, nogil=True):
165 for y in xrange(dy): 214 for y in xrange(dy):
166 for x in xrange(dx): 215 for x in xrange(dx):
167 if mask[z, y, x]: 216 if mask[z, y, x]:
168 - H = calculate_H(aux, z, y, x)  
169 - v = aux[z, y, x] + dt*H  
170 -  
171 - if image[z, y, x]:  
172 - if v < 0:  
173 - out[z, y, x] = 0.00001  
174 - else:  
175 - out[z, y, x] = v  
176 - else:  
177 - if v > 0:  
178 - out[z, y, x] = -0.00001  
179 - else:  
180 - out[z, y, x] = v  
181 -  
182 - diff += (out[z, y, x] - aux[z, y, x])*(out[z, y, x] - aux[z, y, x])  
183 -  
184 - cn = sqrt((1.0/S) * diff)  
185 - print "%d - CN: %.28f - diff: %.28f\n" % (i, cn, diff)  
186 -  
187 - if cn <= E:  
188 - break  
189 -  
190 - return mask 217 + _smooth(image, aux, mask, x, y, z, out)
  218 + # H = calculate_H(aux, z, y, x)
  219 + # v = aux[z, y, x] + dt*H
  220 +
  221 + # if image[z, y, x]:
  222 + # # if v < 0:
  223 + # # out[z, y, x] = 0.00001
  224 + # # else:
  225 + # out[z, y, x] = v
  226 + # else:
  227 + # # if v > 0:
  228 + # # out[z, y, x] = -0.00001
  229 + # # else:
  230 + # out[z, y, x] = v
  231 +
  232 + # diff += (out[z, y, x] - aux[z, y, x])*(out[z, y, x] - aux[z, y, x])
  233 +
  234 + # cn = sqrt((1.0/S) * diff)
  235 + # print "%d - CN: %.28f - diff: %.28f\n" % (i, cn, diff)
  236 + print "Step %d" % i
  237 +
  238 + # if cn <= E:
  239 + # break
  240 +
  241 + return np.asarray(mask)