Download as
Browse Code »

Commit c6a22019532ade5b966532f4ad235d7c816d09bc

Authored by Thiago Franco de Moraes
Committed by GitHub
1 parent f4c771e9
Exists in master

Lanczos 3d (#108) 

* Testing with lanczos 3D

* Using lanczos to interpolate after applying a reorientation

* More options and pointer to function

* lanczos a=4
Showing 5 changed files with 161 additions and 30 deletions   Show diff stats
invesalius/data/interpolation.pxd
  Diff comments   View file @c6a2201
... ... @@ -3,3 +3,6 @@ from .cy_my_types cimport image_t
3 3 cdef double interpolate(image_t[:, :, :], double, double, double) nogil
4 4 cdef double tricub_interpolate(image_t[:, :, :], double, double, double) nogil
5 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
  Diff comments   View file @c6a2201
... ... @@ -4,9 +4,12 @@ import numpy as np
4 4 cimport numpy as np
5 5 cimport cython
6 6  
7   -from libc.math cimport floor, ceil, sqrt, fabs, round
  7 +from libc.math cimport floor, ceil, sqrt, fabs, round, sin, M_PI
8 8 from cython.parallel import prange
9 9  
  10 +DEF LANCZOS_A = 4
  11 +DEF SIZE_LANCZOS_TMP = LANCZOS_A * 2 - 1
  12 +
10 13 cdef double[64][64] temp = [
11 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],
12 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],
... ... @@ -77,6 +80,12 @@ cdef double[64][64] temp = [
77 80 @cython.boundscheck(False) # turn of bounds-checking for entire function
78 81 @cython.cdivision(True)
79 82 @cython.wraparound(False)
  83 +cdef double nearest_neighbour_interp(image_t[:, :, :] V, double x, double y, double z) nogil:
  84 + return V[<int>round(z), <int>round(y), <int>round(x)]
  85 +
  86 +@cython.boundscheck(False) # turn of bounds-checking for entire function
  87 +@cython.cdivision(True)
  88 +@cython.wraparound(False)
80 89 cdef double interpolate(image_t[:, :, :] V, double x, double y, double z) nogil:
81 90 cdef double xd, yd, zd
82 91 cdef double c00, c10, c01, c11
... ... @@ -123,6 +132,76 @@ cdef double interpolate(image_t[:, :, :] V, double x, double y, double z) nogil:
123 132 @cython.boundscheck(False) # turn of bounds-checking for entire function
124 133 @cython.cdivision(True)
125 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)
126 205 cdef image_t _G(image_t[:, :, :] V, int x, int y, int z) nogil:
127 206 cdef int dz, dy, dx
128 207 dz = V.shape[0] - 1
... ... @@ -311,4 +390,3 @@ def tricub_interpolate2_py(image_t[:, :, :] V, double x, double y, double z):
311 390  
312 391 def trilin_interpolate_py(image_t[:, :, :] V, double x, double y, double z):
313 392 return interpolate(V, x, y, z)
314   -
... ...
invesalius/data/slice_.py
  Diff comments   View file @c6a2201
... ... @@ -91,6 +91,8 @@ class Slice(object):
91 91 self._type_projection = const.PROJECTION_NORMAL
92 92 self.n_border = const.PROJECTION_BORDER_SIZE
93 93  
  94 + self.interp_method = 2
  95 +
94 96 self._spacing = (1.0, 1.0, 1.0)
95 97 self.center = [0, 0, 0]
96 98  
... ... @@ -198,6 +200,8 @@ class Slice(object):
198 200  
199 201 Publisher.subscribe(self._fill_holes_auto, 'Fill holes automatically')
200 202  
  203 + Publisher.subscribe(self._set_interpolation_method, 'Set interpolation method')
  204 +
201 205 def GetMaxSliceNumber(self, orientation):
202 206 shape = self.matrix.shape
203 207  
... ... @@ -624,7 +628,8 @@ class Slice(object):
624 628 if orientation == 'AXIAL':
625 629 tmp_array = np.array(self.matrix[slice_number:slice_number + number_slices])
626 630 if np.any(self.q_orientation[1::]):
627   - transforms.apply_view_matrix_transform(self.matrix, self.spacing, M, slice_number, orientation, 2, self.matrix.min(), tmp_array)
  631 + transforms.apply_view_matrix_transform(self.matrix, self.spacing, M, slice_number, orientation, self.interp_method, self.matrix.min(), tmp_array)
  632 + print ">>>", tmp_array.min(), tmp_array.max()
628 633 if self._type_projection == const.PROJECTION_NORMAL:
629 634 n_image = tmp_array.squeeze()
630 635 else:
... ... @@ -671,7 +676,7 @@ class Slice(object):
671 676 elif orientation == 'CORONAL':
672 677 tmp_array = np.array(self.matrix[:, slice_number: slice_number + number_slices, :])
673 678 if np.any(self.q_orientation[1::]):
674   - transforms.apply_view_matrix_transform(self.matrix, self.spacing, M, slice_number, orientation, 2, self.matrix.min(), tmp_array)
  679 + transforms.apply_view_matrix_transform(self.matrix, self.spacing, M, slice_number, orientation, self.interp_method, self.matrix.min(), tmp_array)
675 680  
676 681 if self._type_projection == const.PROJECTION_NORMAL:
677 682 n_image = tmp_array.squeeze()
... ... @@ -721,7 +726,7 @@ class Slice(object):
721 726 elif orientation == 'SAGITAL':
722 727 tmp_array = np.array(self.matrix[:, :, slice_number: slice_number + number_slices])
723 728 if np.any(self.q_orientation[1::]):
724   - transforms.apply_view_matrix_transform(self.matrix, self.spacing, M, slice_number, orientation, 2, self.matrix.min(), tmp_array)
  729 + transforms.apply_view_matrix_transform(self.matrix, self.spacing, M, slice_number, orientation, self.interp_method, self.matrix.min(), tmp_array)
725 730  
726 731 if self._type_projection == const.PROJECTION_NORMAL:
727 732 n_image = tmp_array.squeeze()
... ... @@ -953,6 +958,10 @@ class Slice(object):
953 958 tprojection = pubsub_evt.data
954 959 self.SetTypeProjection(tprojection)
955 960  
  961 + def _set_interpolation_method(self, pubsub_evt):
  962 + interp_method = pubsub_evt.data
  963 + self.SetInterpolationMethod(interp_method)
  964 +
956 965 def SetTypeProjection(self, tprojection):
957 966 if self._type_projection != tprojection:
958 967 if self._type_projection == const.PROJECTION_NORMAL:
... ... @@ -969,6 +978,13 @@ class Slice(object):
969 978  
970 979 Publisher.sendMessage('Check projection menu', tprojection)
971 980  
  981 + def SetInterpolationMethod(self, interp_method):
  982 + if self.interp_method != interp_method:
  983 + self.interp_method = interp_method
  984 + for buffer_ in self.buffer_slices.values():
  985 + buffer_.discard_buffer()
  986 + Publisher.sendMessage('Reload actual slice')
  987 +
972 988 def UpdateWindowLevelBackground(self, pubsub_evt):
973 989 window, level = pubsub_evt.data
974 990 self.window_width = window
... ... @@ -1423,7 +1439,7 @@ class Slice(object):
1423 1439 T1 = transformations.translation_matrix((cz, cy, cx))
1424 1440 M = transformations.concatenate_matrices(T1, R.T, T0)
1425 1441  
1426   - transforms.apply_view_matrix_transform(mcopy, self.spacing, M, 0, 'AXIAL', 2, mcopy.min(), self.matrix)
  1442 + transforms.apply_view_matrix_transform(mcopy, self.spacing, M, 0, 'AXIAL', self.interp_method, mcopy.min(), self.matrix)
1427 1443  
1428 1444 del mcopy
1429 1445 os.remove(temp_file)
... ...
invesalius/data/transforms.pyx
  Diff comments   View file @c6a2201
... ... @@ -3,11 +3,12 @@ cimport numpy as np
3 3 cimport cython
4 4  
5 5 from .cy_my_types cimport image_t
6   -from .interpolation cimport interpolate, tricub_interpolate, tricubicInterpolate
  6 +from .interpolation cimport interpolate, tricub_interpolate, tricubicInterpolate, lanczos3, nearest_neighbour_interp
7 7  
8 8 from libc.math cimport floor, ceil, sqrt, fabs, round
9 9 from cython.parallel import prange
10 10  
  11 +ctypedef double (*interp_function)(image_t[:, :, :], double, double, double) nogil
11 12  
12 13 @cython.boundscheck(False) # turn of bounds-checking for entire function
13 14 @cython.cdivision(True)
... ... @@ -25,7 +26,7 @@ cdef void mul_mat4_vec4(double[:, :] M,
25 26 @cython.boundscheck(False) # turn of bounds-checking for entire function
26 27 @cython.cdivision(True)
27 28 @cython.wraparound(False)
28   -cdef image_t coord_transform(image_t[:, :, :] volume, double[:, :] M, int x, int y, int z, double sx, double sy, double sz, short minterpol, image_t cval) nogil:
  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:
29 30  
30 31 cdef double coord[4]
31 32 coord[0] = z*sz
... ... @@ -53,15 +54,21 @@ cdef image_t coord_transform(image_t[:, :, :] volume, double[:, :] M, int x, int
53 54 cdef double v
54 55  
55 56 if 0 <= nz <= (dz-1) and 0 <= ny <= (dy-1) and 0 <= nx <= (dx-1):
56   - if minterpol == 0:
57   - return volume[<int>round(nz), <int>round(ny), <int>round(nx)]
58   - elif minterpol == 1:
59   - return <image_t>interpolate(volume, nx, ny, nz)
60   - else:
61   - v = tricubicInterpolate(volume, nx, ny, nz)
62   - if (v < cval):
63   - v = cval
64   - return <image_t>v
  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
65 72 else:
66 73 return cval
67 74  
... ... @@ -95,23 +102,34 @@ def apply_view_matrix_transform(image_t[:, :, :] volume,
95 102 sy = spacing[1]
96 103 sz = spacing[2]
97 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 +
98 116 if orientation == 'AXIAL':
99 117 for z in xrange(n, n+odz):
100 118 for y in prange(dy, nogil=True):
101 119 for x in xrange(dx):
102   - out[count, y, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, minterpol, cval)
  120 + out[count, y, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval)
103 121 count += 1
104 122  
105 123 elif orientation == 'CORONAL':
106 124 for y in xrange(n, n+ody):
107 125 for z in prange(dz, nogil=True):
108 126 for x in xrange(dx):
109   - out[z, count, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, minterpol, cval)
  127 + out[z, count, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval)
110 128 count += 1
111 129  
112 130 elif orientation == 'SAGITAL':
113 131 for x in xrange(n, n+odx):
114 132 for z in prange(dz, nogil=True):
115 133 for y in xrange(dy):
116   - out[z, y, count] = coord_transform(volume, M, x, y, z, sx, sy, sz, minterpol,cval)
  134 + out[z, y, count] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval)
117 135 count += 1
... ...
invesalius/gui/dialogs.py
  Diff comments   View file @c6a2201
... ... @@ -1883,26 +1883,38 @@ class ReorientImageDialog(wx.Dialog):
1883 1883 self._bind_events_wx()
1884 1884  
1885 1885 def _init_gui(self):
  1886 + interp_methods_choices = ((_("Nearest Neighbour"), 0),
  1887 + (_("Trilinear"), 1),
  1888 + (_("Tricubic"), 2),
  1889 + (_("Lanczos"), 3))
  1890 + self.interp_method = wx.ComboBox(self, -1, choices=[], style=wx.CB_READONLY)
  1891 + for txt, im_code in interp_methods_choices:
  1892 + self.interp_method.Append(txt, im_code)
  1893 + self.interp_method.SetValue(interp_methods_choices[2][0])
  1894 +
1886 1895 self.anglex = wx.TextCtrl(self, -1, "0.0", style=wx.TE_READONLY)
1887 1896 self.angley = wx.TextCtrl(self, -1, "0.0", style=wx.TE_READONLY)
1888 1897 self.anglez = wx.TextCtrl(self, -1, "0.0", style=wx.TE_READONLY)
1889 1898  
1890 1899 self.btnapply = wx.Button(self, -1, _("Apply"))
1891 1900  
1892   - sizer = wx.BoxSizer(wx.HORIZONTAL)
  1901 + sizer = wx.BoxSizer(wx.VERTICAL)
1893 1902  
1894   - sizer.Add(wx.StaticText(self, -1, _("Angle X")), 0, wx.EXPAND | wx.ALIGN_CENTER_VERTICAL | wx.ALL, 5)
1895   - sizer.Add(self.anglex, 0, wx.EXPAND | wx.ALL, 5)
1896   - sizer.AddSpacer(5)
  1903 + angles_sizer = wx.FlexGridSizer(3, 2, 5, 5)
  1904 + angles_sizer.AddMany([
  1905 + (wx.StaticText(self, -1, _("Angle X")), 1, wx.ALIGN_CENTER_VERTICAL | wx.ALL, 5),
  1906 + (self.anglex, 0, wx.EXPAND | wx.ALL, 5),
1897 1907  
1898   - sizer.Add(wx.StaticText(self, -1, _("Angle Y")), 0, wx.EXPAND | wx.ALIGN_CENTER_VERTICAL | wx.ALL, 5)
1899   - sizer.Add(self.angley, 0, wx.EXPAND | wx.ALL, 5)
1900   - sizer.AddSpacer(5)
  1908 + (wx.StaticText(self, -1, _("Angle Y")), 1, wx.ALIGN_CENTER_VERTICAL | wx.ALL, 5),
  1909 + (self.angley, 0, wx.EXPAND | wx.ALL, 5),
1901 1910  
1902   - sizer.Add(wx.StaticText(self, -1, _("Angle Z")), 0, wx.EXPAND | wx.ALIGN_CENTER_VERTICAL | wx.ALL, 5)
1903   - sizer.Add(self.anglez, 0, wx.EXPAND | wx.ALL, 5)
1904   - sizer.AddSpacer(5)
  1911 + (wx.StaticText(self, -1, _("Angle Z")), 1, wx.ALIGN_CENTER_VERTICAL | wx.ALL, 5),
  1912 + (self.anglez, 0, wx.EXPAND | wx.ALL, 5),
  1913 + ])
1905 1914  
  1915 + sizer.Add(wx.StaticText(self, -1, _("Interpolation method:")), 0, wx.EXPAND | wx.TOP | wx.LEFT | wx.RIGHT, 5)
  1916 + sizer.Add(self.interp_method, 0, wx.EXPAND | wx.ALL, 5)
  1917 + sizer.Add(angles_sizer, 0, wx.EXPAND | wx.ALL, 5)
1906 1918 sizer.Add(self.btnapply, 0, wx.EXPAND | wx.ALL, 5)
1907 1919 sizer.AddSpacer(5)
1908 1920  
... ... @@ -1914,6 +1926,7 @@ class ReorientImageDialog(wx.Dialog):
1914 1926 Publisher.subscribe(self._close_dialog, 'Close reorient dialog')
1915 1927  
1916 1928 def _bind_events_wx(self):
  1929 + self.interp_method.Bind(wx.EVT_COMBOBOX, self.OnSelect)
1917 1930 self.btnapply.Bind(wx.EVT_BUTTON, self.apply_reorientation)
1918 1931 self.Bind(wx.EVT_CLOSE, self.OnClose)
1919 1932  
... ... @@ -1935,6 +1948,9 @@ class ReorientImageDialog(wx.Dialog):
1935 1948 Publisher.sendMessage('Enable style', const.STATE_DEFAULT)
1936 1949 self.Destroy()
1937 1950  
  1951 + def OnSelect(self, evt):
  1952 + im_code = self.interp_method.GetClientData(self.interp_method.GetSelection())
  1953 + Publisher.sendMessage('Set interpolation method', im_code)
1938 1954  
1939 1955  
1940 1956 class ImportBitmapParameters(wx.Dialog):
... ...