diff --git a/invesalius/data/interpolation.pxd b/invesalius/data/interpolation.pxd index 6dde24a..b63cc41 100644 --- a/invesalius/data/interpolation.pxd +++ b/invesalius/data/interpolation.pxd @@ -3,3 +3,6 @@ from .cy_my_types cimport image_t cdef double interpolate(image_t[:, :, :], double, double, double) nogil cdef double tricub_interpolate(image_t[:, :, :], double, double, double) nogil cdef double tricubicInterpolate (image_t[:, :, :], double, double, double) nogil +cdef double lanczos3 (image_t[:, :, :], double, double, double) nogil + +cdef double nearest_neighbour_interp(image_t[:, :, :], double, double, double) nogil diff --git a/invesalius/data/interpolation.pyx b/invesalius/data/interpolation.pyx index 1f7a99e..8aaea8d 100644 --- a/invesalius/data/interpolation.pyx +++ b/invesalius/data/interpolation.pyx @@ -4,9 +4,12 @@ import numpy as np cimport numpy as np cimport cython -from libc.math cimport floor, ceil, sqrt, fabs, round +from libc.math cimport floor, ceil, sqrt, fabs, round, sin, M_PI from cython.parallel import prange +DEF LANCZOS_A = 4 +DEF SIZE_LANCZOS_TMP = LANCZOS_A * 2 - 1 + cdef double[64][64] temp = [ [ 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], [ 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 = [ @cython.boundscheck(False) # turn of bounds-checking for entire function @cython.cdivision(True) @cython.wraparound(False) +cdef double nearest_neighbour_interp(image_t[:, :, :] V, double x, double y, double z) nogil: + return V[round(z), round(y), round(x)] + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) cdef double interpolate(image_t[:, :, :] V, double x, double y, double z) nogil: cdef double xd, yd, zd cdef double c00, c10, c01, c11 @@ -123,6 +132,76 @@ cdef double interpolate(image_t[:, :, :] V, double x, double y, double z) nogil: @cython.boundscheck(False) # turn of bounds-checking for entire function @cython.cdivision(True) @cython.wraparound(False) +cdef inline double lanczos3_L(double x, int a) nogil: + if x == 0: + return 1.0 + elif -a <= x < a: + return (a * sin(M_PI * x) * sin(M_PI * (x / a)))/(M_PI**2 * x**2) + else: + return 0.0 + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) +cdef double lanczos3(image_t[:, :, :] V, double x, double y, double z) nogil: + cdef int a = LANCZOS_A + + cdef int xd = floor(x) + cdef int yd = floor(y) + cdef int zd = floor(z) + + cdef int xi = xd - a + 1 + cdef int xf = xd + a + + cdef int yi = yd - a + 1 + cdef int yf = yd + a + + cdef int zi = zd - a + 1 + cdef int zf = zd + a + + cdef double lx = 0.0 + cdef double ly = 0.0 + cdef double lz = 0.0 + + cdef double[SIZE_LANCZOS_TMP][SIZE_LANCZOS_TMP] temp_x + cdef double[SIZE_LANCZOS_TMP] temp_y + + cdef int i, j, k + cdef int m, n, o + + m = 0 + for k in xrange(zi, zf): + n = 0 + for j in xrange(yi, yf): + lx = 0 + for i in xrange(xi, xf): + lx += _G(V, i, j, k) * lanczos3_L(x - i, a) + temp_x[m][n] = lx + n += 1 + m += 1 + + m = 0 + for k in xrange(zi, zf): + n = 0 + ly = 0 + for j in xrange(yi, yf): + ly += temp_x[m][n] * lanczos3_L(y - j, a) + n += 1 + temp_y[m] = ly + m += 1 + + m = 0 + for k in xrange(zi, zf): + lz += temp_y[m] * lanczos3_L(z - k, a) + m += 1 + + return lz + + +@cython.boundscheck(False) # turn of bounds-checking for entire function +@cython.cdivision(True) +@cython.wraparound(False) cdef image_t _G(image_t[:, :, :] V, int x, int y, int z) nogil: cdef int dz, dy, dx dz = V.shape[0] - 1 @@ -311,4 +390,3 @@ def tricub_interpolate2_py(image_t[:, :, :] V, double x, double y, double z): def trilin_interpolate_py(image_t[:, :, :] V, double x, double y, double z): return interpolate(V, x, y, z) - diff --git a/invesalius/data/slice_.py b/invesalius/data/slice_.py index 656ed42..0a6183d 100644 --- a/invesalius/data/slice_.py +++ b/invesalius/data/slice_.py @@ -91,6 +91,8 @@ class Slice(object): self._type_projection = const.PROJECTION_NORMAL self.n_border = const.PROJECTION_BORDER_SIZE + self.interp_method = 2 + self._spacing = (1.0, 1.0, 1.0) self.center = [0, 0, 0] @@ -198,6 +200,8 @@ class Slice(object): Publisher.subscribe(self._fill_holes_auto, 'Fill holes automatically') + Publisher.subscribe(self._set_interpolation_method, 'Set interpolation method') + def GetMaxSliceNumber(self, orientation): shape = self.matrix.shape @@ -624,7 +628,8 @@ class Slice(object): if orientation == 'AXIAL': tmp_array = np.array(self.matrix[slice_number:slice_number + number_slices]) if np.any(self.q_orientation[1::]): - transforms.apply_view_matrix_transform(self.matrix, self.spacing, M, slice_number, orientation, 2, self.matrix.min(), tmp_array) + transforms.apply_view_matrix_transform(self.matrix, self.spacing, M, slice_number, orientation, self.interp_method, self.matrix.min(), tmp_array) + print ">>>", tmp_array.min(), tmp_array.max() if self._type_projection == const.PROJECTION_NORMAL: n_image = tmp_array.squeeze() else: @@ -671,7 +676,7 @@ class Slice(object): elif orientation == 'CORONAL': tmp_array = np.array(self.matrix[:, slice_number: slice_number + number_slices, :]) if np.any(self.q_orientation[1::]): - transforms.apply_view_matrix_transform(self.matrix, self.spacing, M, slice_number, orientation, 2, self.matrix.min(), tmp_array) + transforms.apply_view_matrix_transform(self.matrix, self.spacing, M, slice_number, orientation, self.interp_method, self.matrix.min(), tmp_array) if self._type_projection == const.PROJECTION_NORMAL: n_image = tmp_array.squeeze() @@ -721,7 +726,7 @@ class Slice(object): elif orientation == 'SAGITAL': tmp_array = np.array(self.matrix[:, :, slice_number: slice_number + number_slices]) if np.any(self.q_orientation[1::]): - transforms.apply_view_matrix_transform(self.matrix, self.spacing, M, slice_number, orientation, 2, self.matrix.min(), tmp_array) + transforms.apply_view_matrix_transform(self.matrix, self.spacing, M, slice_number, orientation, self.interp_method, self.matrix.min(), tmp_array) if self._type_projection == const.PROJECTION_NORMAL: n_image = tmp_array.squeeze() @@ -953,6 +958,10 @@ class Slice(object): tprojection = pubsub_evt.data self.SetTypeProjection(tprojection) + def _set_interpolation_method(self, pubsub_evt): + interp_method = pubsub_evt.data + self.SetInterpolationMethod(interp_method) + def SetTypeProjection(self, tprojection): if self._type_projection != tprojection: if self._type_projection == const.PROJECTION_NORMAL: @@ -969,6 +978,13 @@ class Slice(object): Publisher.sendMessage('Check projection menu', tprojection) + def SetInterpolationMethod(self, interp_method): + if self.interp_method != interp_method: + self.interp_method = interp_method + for buffer_ in self.buffer_slices.values(): + buffer_.discard_buffer() + Publisher.sendMessage('Reload actual slice') + def UpdateWindowLevelBackground(self, pubsub_evt): window, level = pubsub_evt.data self.window_width = window @@ -1423,7 +1439,7 @@ class Slice(object): T1 = transformations.translation_matrix((cz, cy, cx)) M = transformations.concatenate_matrices(T1, R.T, T0) - transforms.apply_view_matrix_transform(mcopy, self.spacing, M, 0, 'AXIAL', 2, mcopy.min(), self.matrix) + transforms.apply_view_matrix_transform(mcopy, self.spacing, M, 0, 'AXIAL', self.interp_method, mcopy.min(), self.matrix) del mcopy os.remove(temp_file) diff --git a/invesalius/data/transforms.pyx b/invesalius/data/transforms.pyx index c8a46f8..a815f23 100644 --- a/invesalius/data/transforms.pyx +++ b/invesalius/data/transforms.pyx @@ -3,11 +3,12 @@ cimport numpy as np cimport cython from .cy_my_types cimport image_t -from .interpolation cimport interpolate, tricub_interpolate, tricubicInterpolate +from .interpolation cimport interpolate, tricub_interpolate, tricubicInterpolate, lanczos3, nearest_neighbour_interp from libc.math cimport floor, ceil, sqrt, fabs, round from cython.parallel import prange +ctypedef double (*interp_function)(image_t[:, :, :], double, double, double) nogil @cython.boundscheck(False) # turn of bounds-checking for entire function @cython.cdivision(True) @@ -25,7 +26,7 @@ cdef void mul_mat4_vec4(double[:, :] M, @cython.boundscheck(False) # turn of bounds-checking for entire function @cython.cdivision(True) @cython.wraparound(False) -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: +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: cdef double coord[4] coord[0] = z*sz @@ -53,15 +54,21 @@ cdef image_t coord_transform(image_t[:, :, :] volume, double[:, :] M, int x, int cdef double v if 0 <= nz <= (dz-1) and 0 <= ny <= (dy-1) and 0 <= nx <= (dx-1): - if minterpol == 0: - return volume[round(nz), round(ny), round(nx)] - elif minterpol == 1: - return interpolate(volume, nx, ny, nz) - else: - v = tricubicInterpolate(volume, nx, ny, nz) - if (v < cval): - v = cval - return v + return f_interp(volume, nx, ny, nz) + # if minterpol == 0: + # return nearest_neighbour_interp(volume, nx, ny, nz) + # elif minterpol == 1: + # return interpolate(volume, nx, ny, nz) + # elif minterpol == 2: + # v = tricubicInterpolate(volume, nx, ny, nz) + # if (v < cval): + # v = cval + # return v + # else: + # v = lanczos3(volume, nx, ny, nz) + # if (v < cval): + # v = cval + # return v else: return cval @@ -95,23 +102,34 @@ def apply_view_matrix_transform(image_t[:, :, :] volume, sy = spacing[1] sz = spacing[2] + cdef interp_function f_interp + + if minterpol == 0: + f_interp = nearest_neighbour_interp + elif minterpol == 1: + f_interp = interpolate + elif minterpol == 2: + f_interp = tricubicInterpolate + else: + f_interp = lanczos3 + if orientation == 'AXIAL': for z in xrange(n, n+odz): for y in prange(dy, nogil=True): for x in xrange(dx): - out[count, y, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, minterpol, cval) + out[count, y, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval) count += 1 elif orientation == 'CORONAL': for y in xrange(n, n+ody): for z in prange(dz, nogil=True): for x in xrange(dx): - out[z, count, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, minterpol, cval) + out[z, count, x] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval) count += 1 elif orientation == 'SAGITAL': for x in xrange(n, n+odx): for z in prange(dz, nogil=True): for y in xrange(dy): - out[z, y, count] = coord_transform(volume, M, x, y, z, sx, sy, sz, minterpol,cval) + out[z, y, count] = coord_transform(volume, M, x, y, z, sx, sy, sz, f_interp, cval) count += 1 diff --git a/invesalius/gui/dialogs.py b/invesalius/gui/dialogs.py index d69f0fc..337a739 100644 --- a/invesalius/gui/dialogs.py +++ b/invesalius/gui/dialogs.py @@ -1883,26 +1883,38 @@ class ReorientImageDialog(wx.Dialog): self._bind_events_wx() def _init_gui(self): + interp_methods_choices = ((_("Nearest Neighbour"), 0), + (_("Trilinear"), 1), + (_("Tricubic"), 2), + (_("Lanczos"), 3)) + self.interp_method = wx.ComboBox(self, -1, choices=[], style=wx.CB_READONLY) + for txt, im_code in interp_methods_choices: + self.interp_method.Append(txt, im_code) + self.interp_method.SetValue(interp_methods_choices[2][0]) + self.anglex = wx.TextCtrl(self, -1, "0.0", style=wx.TE_READONLY) self.angley = wx.TextCtrl(self, -1, "0.0", style=wx.TE_READONLY) self.anglez = wx.TextCtrl(self, -1, "0.0", style=wx.TE_READONLY) self.btnapply = wx.Button(self, -1, _("Apply")) - sizer = wx.BoxSizer(wx.HORIZONTAL) + sizer = wx.BoxSizer(wx.VERTICAL) - sizer.Add(wx.StaticText(self, -1, _("Angle X")), 0, wx.EXPAND | wx.ALIGN_CENTER_VERTICAL | wx.ALL, 5) - sizer.Add(self.anglex, 0, wx.EXPAND | wx.ALL, 5) - sizer.AddSpacer(5) + angles_sizer = wx.FlexGridSizer(3, 2, 5, 5) + angles_sizer.AddMany([ + (wx.StaticText(self, -1, _("Angle X")), 1, wx.ALIGN_CENTER_VERTICAL | wx.ALL, 5), + (self.anglex, 0, wx.EXPAND | wx.ALL, 5), - sizer.Add(wx.StaticText(self, -1, _("Angle Y")), 0, wx.EXPAND | wx.ALIGN_CENTER_VERTICAL | wx.ALL, 5) - sizer.Add(self.angley, 0, wx.EXPAND | wx.ALL, 5) - sizer.AddSpacer(5) + (wx.StaticText(self, -1, _("Angle Y")), 1, wx.ALIGN_CENTER_VERTICAL | wx.ALL, 5), + (self.angley, 0, wx.EXPAND | wx.ALL, 5), - sizer.Add(wx.StaticText(self, -1, _("Angle Z")), 0, wx.EXPAND | wx.ALIGN_CENTER_VERTICAL | wx.ALL, 5) - sizer.Add(self.anglez, 0, wx.EXPAND | wx.ALL, 5) - sizer.AddSpacer(5) + (wx.StaticText(self, -1, _("Angle Z")), 1, wx.ALIGN_CENTER_VERTICAL | wx.ALL, 5), + (self.anglez, 0, wx.EXPAND | wx.ALL, 5), + ]) + sizer.Add(wx.StaticText(self, -1, _("Interpolation method:")), 0, wx.EXPAND | wx.TOP | wx.LEFT | wx.RIGHT, 5) + sizer.Add(self.interp_method, 0, wx.EXPAND | wx.ALL, 5) + sizer.Add(angles_sizer, 0, wx.EXPAND | wx.ALL, 5) sizer.Add(self.btnapply, 0, wx.EXPAND | wx.ALL, 5) sizer.AddSpacer(5) @@ -1914,6 +1926,7 @@ class ReorientImageDialog(wx.Dialog): Publisher.subscribe(self._close_dialog, 'Close reorient dialog') def _bind_events_wx(self): + self.interp_method.Bind(wx.EVT_COMBOBOX, self.OnSelect) self.btnapply.Bind(wx.EVT_BUTTON, self.apply_reorientation) self.Bind(wx.EVT_CLOSE, self.OnClose) @@ -1935,6 +1948,9 @@ class ReorientImageDialog(wx.Dialog): Publisher.sendMessage('Enable style', const.STATE_DEFAULT) self.Destroy() + def OnSelect(self, evt): + im_code = self.interp_method.GetClientData(self.interp_method.GetSelection()) + Publisher.sendMessage('Set interpolation method', im_code) class ImportBitmapParameters(wx.Dialog): -- libgit2 0.21.2