Added a cython version to threshold

Thiago Franco de Moraes
1 parent 3867177c
Showing 3 changed files with 48 additions and 7 deletions Show diff stats
invesalius/data/slice_.py
invesalius/data/threshold.pyx
setup.py
@@ -35,6 +35,8 @@ from mask import Mask
 from project import Project
 from data import mips
  
+from data import threshold
+
 OTHER=0
 PLIST=1
 WIDGET=2
@@ -1206,12 +1208,14 @@ class Slice(object):
         given slice_matrix.
         """
         thresh_min, thresh_max = self.current_mask.threshold_range
-        m = (((slice_matrix >= thresh_min) & (slice_matrix <= thresh_max)) * 255)
-        m[mask == 1] = 1
-        m[mask == 2] = 2
-        m[mask == 253] = 253
-        m[mask == 254] = 254
-        return m.astype('uint8')
+        #m = (((slice_matrix >= thresh_min) & (slice_matrix <= thresh_max)) * 255)
+        #m[mask == 1] = 1
+        #m[mask == 2] = 2
+        #m[mask == 253] = 253
+        #m[mask == 254] = 254
+        m = numpy.zeros_like(mask, dtype='uint8')
+        threshold.threshold(slice_matrix, m, thresh_min, thresh_max)
+        return m
  
     def do_threshold_to_all_slices(self):
         mask = self.current_mask
@@ -0,0 +1,32 @@
+import numpy as np
+cimport numpy as np
+cimport cython
+
+from libc.math cimport floor, ceil, sqrt, fabs
+from cython.parallel import prange
+
+DTYPE8 = np.uint8
+ctypedef np.uint8_t DTYPE8_t
+
+DTYPE16 = np.int16
+ctypedef np.int16_t DTYPE16_t
+
+DTYPEF32 = np.float32
+ctypedef np.float32_t DTYPEF32_t
+
+@cython.boundscheck(False) # turn of bounds-checking for entire function
+@cython.cdivision(True)
+@cython.wraparound(False)
+def threshold(DTYPE16_t[:, :] image, DTYPE8_t[:, :] mask, DTYPE16_t low, DTYPE16_t high):
+    cdef int sy = image.shape[0]
+    cdef int sx = image.shape[1]
+    cdef int x, y
+    cdef DTYPE16_t v
+    for y in prange(sy, nogil=True):
+        for x in xrange(sx):
+            v = mask[y, x]
+            if not v:
+                if v >= low and v <= high:
+                    mask[y, x] = 255
+                else:
+                    mask[y, x] = 0
@@ -12,7 +12,12 @@ if sys.platform == &#39;linux2&#39;:
         ext_modules = [ Extension("invesalius.data.mips", ["invesalius/data/mips.pyx"],
                                   include_dirs =  [numpy.get_include()],
                                   extra_compile_args=['-fopenmp'],
-                                  extra_link_args=['-fopenmp'],)]
+                                  extra_link_args=['-fopenmp'],),
+                       Extension("invesalius.data.threshold", ["invesalius/data/threshold.pyx"],
+                                  include_dirs =  [numpy.get_include()],
+                                  extra_compile_args=['-fopenmp'],
+                                  extra_link_args=['-fopenmp'],),
+                      ]
          )
  
 elif sys.platform == 'win32':
...	...	@@ -35,6 +35,8 @@ from mask import Mask
35	35	from project import Project
36	36	from data import mips
37	37
	38	+from data import threshold
	39	+
38	40	OTHER=0
39	41	PLIST=1
40	42	WIDGET=2
...	...	@@ -1206,12 +1208,14 @@ class Slice(object):
1206	1208	given slice_matrix.
1207	1209	"""
1208	1210	thresh_min, thresh_max = self.current_mask.threshold_range
1209		- m = (((slice_matrix >= thresh_min) & (slice_matrix <= thresh_max)) * 255)
1210		- m[mask == 1] = 1
1211		- m[mask == 2] = 2
1212		- m[mask == 253] = 253
1213		- m[mask == 254] = 254
1214		- return m.astype('uint8')
	1211	+ #m = (((slice_matrix >= thresh_min) & (slice_matrix <= thresh_max)) * 255)
	1212	+ #m[mask == 1] = 1
	1213	+ #m[mask == 2] = 2
	1214	+ #m[mask == 253] = 253
	1215	+ #m[mask == 254] = 254
	1216	+ m = numpy.zeros_like(mask, dtype='uint8')
	1217	+ threshold.threshold(slice_matrix, m, thresh_min, thresh_max)
	1218	+ return m
1215	1219
1216	1220	def do_threshold_to_all_slices(self):
1217	1221	mask = self.current_mask
...	...
...	...	@@ -0,0 +1,32 @@
	1	+import numpy as np
	2	+cimport numpy as np
	3	+cimport cython
	4	+
	5	+from libc.math cimport floor, ceil, sqrt, fabs
	6	+from cython.parallel import prange
	7	+
	8	+DTYPE8 = np.uint8
	9	+ctypedef np.uint8_t DTYPE8_t
	10	+
	11	+DTYPE16 = np.int16
	12	+ctypedef np.int16_t DTYPE16_t
	13	+
	14	+DTYPEF32 = np.float32
	15	+ctypedef np.float32_t DTYPEF32_t
	16	+
	17	+@cython.boundscheck(False) # turn of bounds-checking for entire function
	18	+@cython.cdivision(True)
	19	+@cython.wraparound(False)
	20	+def threshold(DTYPE16_t[:, :] image, DTYPE8_t[:, :] mask, DTYPE16_t low, DTYPE16_t high):
	21	+ cdef int sy = image.shape[0]
	22	+ cdef int sx = image.shape[1]
	23	+ cdef int x, y
	24	+ cdef DTYPE16_t v
	25	+ for y in prange(sy, nogil=True):
	26	+ for x in xrange(sx):
	27	+ v = mask[y, x]
	28	+ if not v:
	29	+ if v >= low and v <= high:
	30	+ mask[y, x] = 255
	31	+ else:
	32	+ mask[y, x] = 0
...	...
...	...	@@ -12,7 +12,12 @@ if sys.platform == 'linux2':
12	12	ext_modules = [ Extension("invesalius.data.mips", ["invesalius/data/mips.pyx"],
13	13	include_dirs = [numpy.get_include()],
14	14	extra_compile_args=['-fopenmp'],
15		- extra_link_args=['-fopenmp'],)]
	15	+ extra_link_args=['-fopenmp'],),
	16	+ Extension("invesalius.data.threshold", ["invesalius/data/threshold.pyx"],
	17	+ include_dirs = [numpy.get_include()],
	18	+ extra_compile_args=['-fopenmp'],
	19	+ extra_link_args=['-fopenmp'],),
	20	+ ]
16	21	)
17	22
18	23	elif sys.platform == 'win32':
...	...