From 8af3aa25a11d2c8e2ef92ce1286b03d1b3ea2d0f Mon Sep 17 00:00:00 2001 From: Yousef Moazzam Date: Thu, 19 Jan 2023 11:18:10 +0000 Subject: [PATCH 1/3] Port large stripe removal from TomoPy to CuPy Similar to other ports to CuPy, the residual is non-zero between the CuPy version and the CPU TomoPy version, but small. Notably, the result of both implementations when using default values for all parameters produces projection images with noisy artifacts. It is primarily the regions of noise which have the greatest difference in the residual. Currently, the reason for these artifacts is unknown, but one possibility is that this particular stripe removal method is not suited for the test data (ie, possibly no "large" stripes exist in the test data). --- src/httomolib/prep/stripe.py | 118 ++++++++++++++++++++++++++++++++++- 1 file changed, 115 insertions(+), 3 deletions(-) diff --git a/src/httomolib/prep/stripe.py b/src/httomolib/prep/stripe.py index 1a867524..83f4cafa 100644 --- a/src/httomolib/prep/stripe.py +++ b/src/httomolib/prep/stripe.py @@ -23,7 +23,7 @@ import cupy as cp import numpy as np from cupy import abs, mean, ndarray -from cupyx.scipy.ndimage import median_filter +from cupyx.scipy.ndimage import median_filter, binary_dilation __all__ = [ 'detect_stripes', @@ -32,8 +32,120 @@ 'remove_stripes_titarenko_cupy', ] -# TODO: port 'remove_all_stripe', 'remove_large_stripe' and 'remove_dead_stripe' -# from https://github.com/tomopy/tomopy/blob/master/source/tomopy/prep/stripe.py +# TODO: port 'remove_all_stripe' and 'remove_dead_stripe' from +# https://github.com/tomopy/tomopy/blob/master/source/tomopy/prep/stripe.py + + +## %%%%%%%%%%%%%%%%%%%%% remove_large_stripe_cupy %%%%%%%%%%%%%%%%%%%%%%%%% ## +## Naive CuPy port of the NumPy implementation in TomoPy +def remove_large_stripe_cupy(tomo: ndarray, snr: float=3, size: int=51, + drop_ratio: float=0.1, norm: bool=True) -> ndarray: + """ + Remove large stripe artifacts from sinogram using Nghia Vo's + approach :cite:`Vo:18` (algorithm 5). + + Parameters + ---------- + tomo : ndarray + 3D tomographic data. + + snr : float, optional + Ratio used to locate of large stripes. + Greater is less sensitive. + + size : int, optional + Window size of the median filter. + + drop_ratio : float, optional + Ratio of pixels to be dropped, which is used to reduce the false + detection of stripes. + + norm : bool, optional + Apply normalization if True. + + Returns + ------- + ndarray + Corrected 3D tomographic data. + """ + matindex = _create_matindex(tomo.shape[2], tomo.shape[0]) + for m in range(tomo.shape[1]): + sino = tomo[:, m, :] + tomo[:, m, :] = _rs_large(sino, snr, size, matindex, drop_ratio, norm) + + return tomo + + +def _rs_large(sinogram: ndarray, snr: float, size: int, matindex: ndarray, + drop_ratio: float=0.1, norm: bool=True) -> ndarray: + drop_ratio = cp.clip(cp.asarray(drop_ratio, dtype=cp.float32), 0.0, 0.8) + (nrow, _) = sinogram.shape + ndrop = int(0.5 * drop_ratio * nrow) + # Note: CuPy's docs + # https://docs.cupy.dev/en/stable/reference/generated/cupy.sort.html refer + # to the default option of `kind=None` being a stable algorithm. NumPy docs + # https://numpy.org/doc/stable/reference/generated/numpy.sort.html#numpy.sort + # on the other hand when given the default value `kind=None` uses quicksort. + sinosort = cp.sort(sinogram, axis=0) + sinosmooth = median_filter(sinosort, (1, size)) + list1 = mean(sinosort[ndrop:nrow - ndrop], axis=0) + list2 = mean(sinosmooth[ndrop:nrow - ndrop], axis=0) + # TODO: Using the `out` parameter in conjunction with the `where` parameter + # could decrease memory usage via avoiding creating new array objects; + # however, something isn't quite working with the value being passed for the + # `where` parameter, requires a bit of investigation. + listfact = cp.divide(list1, list2, + #out=cp.ones_like(list1), + #where=list2 != 0 + ) + + # Locate stripes + listmask = _detect_stripe(listfact, snr) + listmask = binary_dilation(listmask, iterations=1).astype(listmask.dtype) + matfact = cp.tile(listfact, (nrow, 1)) + + # Normalize + if norm is True: + sinogram = sinogram / matfact + sinogram1 = cp.transpose(sinogram) + matcombine = cp.asarray(cp.dstack((matindex, sinogram1))) + matsort = cp.asarray( + [row[row[:, 1].argsort()] for row in matcombine]) + matsort[:, :, 1] = cp.transpose(sinosmooth) + matsortback = cp.asarray( + [row[row[:, 0].argsort()] for row in matsort]) + sino_corrected = cp.transpose(matsortback[:, :, 1]) + # TODO: Taking into account NumPy docs for `np.where()` + # https://numpy.org/doc/stable/reference/generated/numpy.where.html#numpy.where + # and CuPy docs for `cp.where()` + # https://docs.cupy.dev/en/stable/reference/generated/cupy.where.html?highlight=where, + # a possibly better alternative is: + # listxmiss = cp.asarray(listmask > 0.0).nonzero()[0] + listxmiss = cp.where(listmask > 0.0)[0] + sinogram[:, listxmiss] = sino_corrected[:, listxmiss] + return sinogram + + +def _detect_stripe(listdata: ndarray, snr: float) -> ndarray: + numdata = len(listdata) + listsorted = cp.sort(listdata)[::-1] + xlist = cp.arange(0, numdata, 1.0) + ndrop = cp.int16(0.25 * numdata) + (_slope, _intercept) = cp.polyfit( + xlist[ndrop:-ndrop - 1], listsorted[ndrop:-ndrop - 1], 1) + numt1 = _intercept + _slope * xlist[-1] + noiselevel = abs(numt1 - _intercept) + noiselevel = cp.clip(cp.asarray(noiselevel, dtype=cp.float32), 1e-6, None) + val1 = abs(listsorted[0] - _intercept) / noiselevel + val2 = abs(listsorted[-1] - numt1) / noiselevel + listmask = cp.zeros_like(listdata) + if (val1 >= snr): + upper_thresh = _intercept + noiselevel * snr * 0.5 + listmask[listdata > upper_thresh] = 1.0 + if (val2 >= snr): + lower_thresh = numt1 - noiselevel * snr * 0.5 + listmask[listdata <= lower_thresh] = 1.0 + return listmask ## %%%%%%%%%%%%%%%%% remove_stripe_based_sorting_cupy %%%%%%%%%%%%%%%%%%%%% ## From e18feafedb3c0f3aacda8d93dc34d9a6c4aee227 Mon Sep 17 00:00:00 2001 From: Yousef Moazzam Date: Fri, 20 Jan 2023 10:09:08 +0000 Subject: [PATCH 2/3] Port dead stripe removal from TomoPy to CuPy Similar comments here as the large stripe removal port in 8af3aa2 regarding the residual between the TomoPy/CPU implementation and the CuPy/GPU implementation. Something else worth noting is that the TomoPy version of this method uses a 2D interpolation function from scipy https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp2d.html. The current latest version of CuPy (v11.5.0) does not have an analogue of this function, but the pre-release version v.12.0.0b3 does, it has an n-dimensional interpolation function that serves well as a GPU replacement of the scipy `interp2d()` function https://github.com/cupy/cupy/releases/tag/v12.0.0b3. The required function `interpn()` in the CuPy codebase and its dependent functions https://github.com/cupy/cupy/blob/1a9c91411fe6297ea81c6ad103cd8211134f6446/cupyx/scipy/interpolate/_rgi.py have been copied into the httomolib codebase for the moment, as it was the simplest way to incorporate the required interpolation function. --- src/httomolib/prep/_rgi.py | 536 +++++++++++++++++++++++++++++++++++ src/httomolib/prep/stripe.py | 79 +++++- 2 files changed, 613 insertions(+), 2 deletions(-) create mode 100644 src/httomolib/prep/_rgi.py diff --git a/src/httomolib/prep/_rgi.py b/src/httomolib/prep/_rgi.py new file mode 100644 index 00000000..ae289d32 --- /dev/null +++ b/src/httomolib/prep/_rgi.py @@ -0,0 +1,536 @@ +# Code taken from +# https://github.com/cupy/cupy/blob/1a9c91411fe6297ea81c6ad103cd8211134f6446/cupyx/scipy/interpolate/_rgi.py +# since the `interpn()` function is only available in CuPy v12.0.0b3 pre-release +# https://github.com/cupy/cupy/releases/tag/v12.0.0b3 +__all__ = ['RegularGridInterpolator', 'interpn'] + +import itertools +import cupy as cp + + +def _ndim_coords_from_arrays(points, ndim=None): + """ + Convert a tuple of coordinate arrays to a (..., ndim)-shaped array. + """ + if isinstance(points, tuple) and len(points) == 1: + # handle argument tuple + points = points[0] + if isinstance(points, tuple): + p = cp.broadcast_arrays(*points) + n = len(p) + for j in range(1, n): + if p[j].shape != p[0].shape: + raise ValueError( + "coordinate arrays do not have the same shape") + points = cp.empty(p[0].shape + (len(points),), dtype=float) + for j, item in enumerate(p): + points[..., j] = item + else: + points = cp.asanyarray(points) + if points.ndim == 1: + if ndim is None: + points = points.reshape(-1, 1) + else: + points = points.reshape(-1, ndim) + return points + + +def _check_points(points): + descending_dimensions = [] + grid = [] + for i, p in enumerate(points): + # early make points float + # see https://github.com/scipy/scipy/pull/17230 + p = cp.asarray(p, dtype=float) + if not cp.all(p[1:] > p[:-1]): + if cp.all(p[1:] < p[:-1]): + # input is descending, so make it ascending + descending_dimensions.append(i) + p = cp.flip(p) + p = cp.ascontiguousarray(p) + else: + raise ValueError( + "The points in dimension %d must be strictly " + "ascending or descending" % i) + grid.append(p) + return tuple(grid), tuple(descending_dimensions) + + +def _check_dimensionality(points, values): + if len(points) > values.ndim: + raise ValueError("There are %d point arrays, but values has %d " + "dimensions" % (len(points), values.ndim)) + for i, p in enumerate(points): + if not cp.asarray(p).ndim == 1: + raise ValueError("The points in dimension %d must be " + "1-dimensional" % i) + if not values.shape[i] == len(p): + raise ValueError("There are %d points and %d values in " + "dimension %d" % (len(p), values.shape[i], i)) + + +class RegularGridInterpolator: + """ + Interpolation on a regular or rectilinear grid in arbitrary dimensions. + + The data must be defined on a rectilinear grid; that is, a rectangular + grid with even or uneven spacing. Linear and nearest-neighbor + interpolations are supported. After setting up the interpolator object, + the interpolation method may be chosen at each evaluation. + + Parameters + ---------- + points : tuple of ndarray of float, with shapes (m1, ), ..., (mn, ) + The points defining the regular grid in n dimensions. The points in + each dimension (i.e. every elements of the points tuple) must be + strictly ascending or descending. + + values : ndarray, shape (m1, ..., mn, ...) + The data on the regular grid in n dimensions. Complex data can be + acceptable. + + method : str, optional + The method of interpolation to perform. Supported are "linear" and + "nearest". This parameter will become the default for the object's + ``__call__`` method. Default is "linear". + + bounds_error : bool, optional + If True, when interpolated values are requested outside of the + domain of the input data, a ValueError is raised. + If False, then `fill_value` is used. + Default is True. + + fill_value : float or None, optional + The value to use for points outside of the interpolation domain. + If None, values outside the domain are extrapolated. + Default is ``cp.nan``. + + Notes + ----- + Contrary to scipy's `LinearNDInterpolator` and `NearestNDInterpolator`, + this class avoids expensive triangulation of the input data by taking + advantage of the regular grid structure. + + In other words, this class assumes that the data is defined on a + *rectilinear* grid. + + If the input data is such that dimensions have incommensurate + units and differ by many orders of magnitude, the interpolant may have + numerical artifacts. Consider rescaling the data before interpolating. + + Examples + -------- + **Evaluate a function on the points of a 3-D grid** + + As a first example, we evaluate a simple example function on the points of + a 3-D grid: + + >>> from cupyx.scipy.interpolate import RegularGridInterpolator + >>> import cupy as cp + >>> def f(x, y, z): + ... return 2 * x**3 + 3 * y**2 - z + >>> x = cp.linspace(1, 4, 11) + >>> y = cp.linspace(4, 7, 22) + >>> z = cp.linspace(7, 9, 33) + >>> xg, yg ,zg = cp.meshgrid(x, y, z, indexing='ij', sparse=True) + >>> data = f(xg, yg, zg) + + ``data`` is now a 3-D array with ``data[i, j, k] = f(x[i], y[j], z[k])``. + Next, define an interpolating function from this data: + + >>> interp = RegularGridInterpolator((x, y, z), data) + + Evaluate the interpolating function at the two points + ``(x,y,z) = (2.1, 6.2, 8.3)`` and ``(3.3, 5.2, 7.1)``: + + >>> pts = cp.array([[2.1, 6.2, 8.3], + ... [3.3, 5.2, 7.1]]) + >>> interp(pts) + array([ 125.80469388, 146.30069388]) + + which is indeed a close approximation to + + >>> f(2.1, 6.2, 8.3), f(3.3, 5.2, 7.1) + (125.54200000000002, 145.894) + + **Interpolate and extrapolate a 2D dataset** + + As a second example, we interpolate and extrapolate a 2D data set: + + >>> x, y = cp.array([-2, 0, 4]), cp.array([-2, 0, 2, 5]) + >>> def ff(x, y): + ... return x**2 + y**2 + + >>> xg, yg = cp.meshgrid(x, y, indexing='ij') + >>> data = ff(xg, yg) + >>> interp = RegularGridInterpolator((x, y), data, + ... bounds_error=False, fill_value=None) + + >>> import matplotlib.pyplot as plt + >>> fig = plt.figure() + >>> ax = fig.add_subplot(projection='3d') + >>> ax.scatter(xg.ravel().get(), yg.ravel().get(), data.ravel().get(), + ... s=60, c='k', label='data') + + Evaluate and plot the interpolator on a finer grid + + >>> xx = cp.linspace(-4, 9, 31) + >>> yy = cp.linspace(-4, 9, 31) + >>> X, Y = cp.meshgrid(xx, yy, indexing='ij') + + >>> # interpolator + >>> ax.plot_wireframe(X.get(), Y.get(), interp((X, Y)).get(), + rstride=3, cstride=3, alpha=0.4, color='m', + label='linear interp') + + >>> # ground truth + >>> ax.plot_wireframe(X.get(), Y.get(), ff(X, Y).get(), + rstride=3, cstride=3, + ... alpha=0.4, label='ground truth') + >>> plt.legend() + >>> plt.show() + + See Also + -------- + interpn : a convenience function which wraps `RegularGridInterpolator` + + scipy.ndimage.map_coordinates : interpolation on grids with equal spacing + (suitable for e.g., N-D image resampling) + + References + ---------- + [1] Python package *regulargrid* by Johannes Buchner, see + https://pypi.python.org/pypi/regulargrid/ + [2] Wikipedia, "Trilinear interpolation", + https://en.wikipedia.org/wiki/Trilinear_interpolation + [3] Weiser, Alan, and Sergio E. Zarantonello. "A note on piecewise + linear and multilinear table interpolation in many dimensions." + MATH. COMPUT. 50.181 (1988): 189-196. + https://www.ams.org/journals/mcom/1988-50-181/S0025-5718-1988-0917826-0/S0025-5718-1988-0917826-0.pdf + """ + # this class is based on code originally programmed by Johannes Buchner, + # see https://github.com/JohannesBuchner/regulargrid + + _ALL_METHODS = ["linear", "nearest"] + + def __init__(self, points, values, method="linear", bounds_error=True, + fill_value=cp.nan): + if method not in self._ALL_METHODS: + raise ValueError("Method '%s' is not defined" % method) + + self.method = method + self.bounds_error = bounds_error + self.grid, self._descending_dimensions = _check_points(points) + self.values = self._check_values(values) + self._check_dimensionality(self.grid, self.values) + self.fill_value = self._check_fill_value(self.values, fill_value) + if self._descending_dimensions: + self.values = cp.flip(values, axis=self._descending_dimensions) + + def _check_dimensionality(self, grid, values): + _check_dimensionality(grid, values) + + def _check_points(self, points): + return _check_points(points) + + def _check_values(self, values): + if not cp.issubdtype(values.dtype, cp.inexact): + values = values.astype(float) + + return values + + def _check_fill_value(self, values, fill_value): + if fill_value is not None: + fill_value_dtype = cp.asarray(fill_value).dtype + if (hasattr(values, 'dtype') and + not cp.can_cast(fill_value_dtype, values.dtype, + casting='same_kind')): + raise ValueError("fill_value must be either 'None' or " + "of a type compatible with values") + return fill_value + + def __call__(self, xi, method=None): + """ + Interpolation at coordinates. + + Parameters + ---------- + xi : cupy.ndarray of shape (..., ndim) + The coordinates to evaluate the interpolator at. + + method : str, optional + The method of interpolation to perform. Supported are "linear" and + "nearest". Default is the method chosen when the interpolator was + created. + + Returns + ------- + values_x : cupy.ndarray, shape xi.shape[:-1] + values.shape[ndim:] + Interpolated values at `xi`. See notes for behaviour when + ``xi.ndim == 1``. + + Notes + ----- + In the case that ``xi.ndim == 1`` a new axis is inserted into + the 0 position of the returned array, values_x, so its shape is + instead ``(1,) + values.shape[ndim:]``. + + Examples + -------- + Here we define a nearest-neighbor interpolator of a simple function + + >>> import cupy as cp + >>> x, y = cp.array([0, 1, 2]), cp.array([1, 3, 7]) + >>> def f(x, y): + ... return x**2 + y**2 + >>> data = f(*cp.meshgrid(x, y, indexing='ij', sparse=True)) + >>> from cupyx.scipy.interpolate import RegularGridInterpolator + >>> interp = RegularGridInterpolator((x, y), data, method='nearest') + + By construction, the interpolator uses the nearest-neighbor + interpolation + + >>> interp([[1.5, 1.3], [0.3, 4.5]]) + array([2., 9.]) + + We can however evaluate the linear interpolant by overriding the + `method` parameter + + >>> interp([[1.5, 1.3], [0.3, 4.5]], method='linear') + array([ 4.7, 24.3]) + """ + method = self.method if method is None else method + if method not in self._ALL_METHODS: + raise ValueError("Method '%s' is not defined" % method) + + xi, xi_shape, ndim, nans, out_of_bounds = self._prepare_xi(xi) + + if method == "linear": + indices, norm_distances = self._find_indices(xi.T) + result = self._evaluate_linear(indices, norm_distances) + elif method == "nearest": + indices, norm_distances = self._find_indices(xi.T) + result = self._evaluate_nearest(indices, norm_distances) + + if not self.bounds_error and self.fill_value is not None: + result[out_of_bounds] = self.fill_value + + if nans.ndim < result.ndim: + nans = nans[..., None] + result = cp.where(nans, cp.nan, result) + return result.reshape(xi_shape[:-1] + self.values.shape[ndim:]) + + def _prepare_xi(self, xi): + ndim = len(self.grid) + xi = _ndim_coords_from_arrays(xi, ndim=ndim) + if xi.shape[-1] != len(self.grid): + raise ValueError("The requested sample points xi have dimension " + f"{xi.shape[-1]} but this " + f"RegularGridInterpolator has dimension {ndim}") + + xi_shape = xi.shape + xi = xi.reshape(-1, xi_shape[-1]) + xi = cp.asarray(xi, dtype=float) + + # find nans in input + is_nans = cp.isnan(xi).T + nans = is_nans[0].copy() + for is_nan in is_nans[1:]: + cp.logical_or(nans, is_nan, nans) + + if self.bounds_error: + for i, p in enumerate(xi.T): + if not cp.logical_and(cp.all(self.grid[i][0] <= p), + cp.all(p <= self.grid[i][-1])): + raise ValueError("One of the requested xi is out of bounds" + " in dimension %d" % i) + out_of_bounds = None + else: + out_of_bounds = self._find_out_of_bounds(xi.T) + + return xi, xi_shape, ndim, nans, out_of_bounds + + def _evaluate_linear(self, indices, norm_distances): + # slice for broadcasting over trailing dimensions in self.values + vslice = (slice(None),) + (None,)*(self.values.ndim - len(indices)) + + # Compute shifting up front before zipping everything together + shift_norm_distances = [1 - yi for yi in norm_distances] + shift_indices = [i + 1 for i in indices] + + # The formula for linear interpolation in 2d takes the form: + # values = self.values[(i0, i1)] * (1 - y0) * (1 - y1) + \ + # self.values[(i0, i1 + 1)] * (1 - y0) * y1 + \ + # self.values[(i0 + 1, i1)] * y0 * (1 - y1) + \ + # self.values[(i0 + 1, i1 + 1)] * y0 * y1 + # We pair i with 1 - yi (zipped1) and i + 1 with yi (zipped2) + zipped1 = zip(indices, shift_norm_distances) + zipped2 = zip(shift_indices, norm_distances) + + # Take all products of zipped1 and zipped2 and iterate over them + # to get the terms in the above formula. This corresponds to iterating + # over the vertices of a hypercube. + hypercube = itertools.product(*zip(zipped1, zipped2)) + value = cp.array([0.]) + for h in hypercube: + edge_indices, weights = zip(*h) + term = cp.asarray(self.values[edge_indices]) + for w in weights: + term *= w[vslice] + value = value + term # cannot use += because broadcasting + return value + + def _evaluate_nearest(self, indices, norm_distances): + idx_res = [cp.where(yi <= .5, i, i + 1) + for i, yi in zip(indices, norm_distances)] + return self.values[tuple(idx_res)] + + def _find_indices(self, xi): + # find relevant edges between which xi are situated + indices = [] + # compute distance to lower edge in unity units + norm_distances = [] + # iterate through dimensions + for x, grid in zip(xi, self.grid): + i = cp.searchsorted(grid, x) - 1 + cp.clip(i, 0, grid.size - 2, i) + indices.append(i) + + # compute norm_distances, incl length-1 grids, + # where `grid[i+1] == grid[i]` + denom = grid[i + 1] - grid[i] + norm_dist = cp.where(denom != 0, (x - grid[i]) / denom, 0) + norm_distances.append(norm_dist) + + return indices, norm_distances + + def _find_out_of_bounds(self, xi): + # check for out of bounds xi + out_of_bounds = cp.zeros((xi.shape[1]), dtype=bool) + # iterate through dimensions + for x, grid in zip(xi, self.grid): + out_of_bounds += x < grid[0] + out_of_bounds += x > grid[-1] + return out_of_bounds + + +def interpn(points, values, xi, method="linear", bounds_error=True, + fill_value=cp.nan): + """ + Multidimensional interpolation on regular or rectilinear grids. + + Strictly speaking, not all regular grids are supported - this function + works on *rectilinear* grids, that is, a rectangular grid with even or + uneven spacing. + + Parameters + ---------- + points : tuple of cupy.ndarray of float, with shapes (m1, ), ..., (mn, ) + The points defining the regular grid in n dimensions. The points in + each dimension (i.e. every elements of the points tuple) must be + strictly ascending or descending. + + values : cupy.ndarray of shape (m1, ..., mn, ...) + The data on the regular grid in n dimensions. Complex data can be + acceptable. + + xi : cupy.ndarray of shape (..., ndim) + The coordinates to sample the gridded data at + + method : str, optional + The method of interpolation to perform. Supported are "linear" and + "nearest". + + bounds_error : bool, optional + If True, when interpolated values are requested outside of the + domain of the input data, a ValueError is raised. + If False, then `fill_value` is used. + + fill_value : number, optional + If provided, the value to use for points outside of the + interpolation domain. If None, values outside + the domain are extrapolated. + + Returns + ------- + values_x : ndarray, shape xi.shape[:-1] + values.shape[ndim:] + Interpolated values at `xi`. See notes for behaviour when + ``xi.ndim == 1``. + + Notes + ----- + + In the case that ``xi.ndim == 1`` a new axis is inserted into + the 0 position of the returned array, values_x, so its shape is + instead ``(1,) + values.shape[ndim:]``. + + If the input data is such that input dimensions have incommensurate + units and differ by many orders of magnitude, the interpolant may have + numerical artifacts. Consider rescaling the data before interpolation. + + Examples + -------- + Evaluate a simple example function on the points of a regular 3-D grid: + + >>> import cupy as cp + >>> from cupyx.scipy.interpolate import interpn + >>> def value_func_3d(x, y, z): + ... return 2 * x + 3 * y - z + >>> x = cp.linspace(0, 4, 5) + >>> y = cp.linspace(0, 5, 6) + >>> z = cp.linspace(0, 6, 7) + >>> points = (x, y, z) + >>> values = value_func_3d(*cp.meshgrid(*points, indexing='ij')) + + Evaluate the interpolating function at a point + + >>> point = cp.array([2.21, 3.12, 1.15]) + >>> print(interpn(points, values, point)) + [12.63] + + See Also + -------- + RegularGridInterpolator : interpolation on a regular or rectilinear grid + in arbitrary dimensions (`interpn` wraps this + class). + + cupyx.scipy.ndimage.map_coordinates : interpolation on grids with equal + spacing (suitable for e.g., N-D image + resampling) + """ + # sanity check 'method' kwarg + if method not in ["linear", "nearest"]: + raise ValueError( + "interpn only understands the methods 'linear' and 'nearest'. " + "You provided {method}.") + + ndim = values.ndim + + # sanity check consistency of input dimensions + if len(points) > ndim: + raise ValueError("There are %d point arrays, but values has %d " + "dimensions" % (len(points), ndim)) + + grid, descending_dimensions = _check_points(points) + _check_dimensionality(grid, values) + + # sanity check requested xi + xi = _ndim_coords_from_arrays(xi, ndim=len(grid)) + if xi.shape[-1] != len(grid): + raise ValueError("The requested sample points xi have dimension " + "%d, but this RegularGridInterpolator has " + "dimension %d" % (xi.shape[-1], len(grid))) + + if bounds_error: + for i, p in enumerate(xi.T): + if not cp.logical_and(cp.all(grid[i][0] <= p), + cp.all(p <= grid[i][-1])): + raise ValueError("One of the requested xi is out of bounds " + "in dimension %d" % i) + + # perform interpolation + if method in ["linear", "nearest"]: + interp = RegularGridInterpolator(points, values, method=method, + bounds_error=bounds_error, + fill_value=fill_value) + return interp(xi) diff --git a/src/httomolib/prep/stripe.py b/src/httomolib/prep/stripe.py index 83f4cafa..893923b4 100644 --- a/src/httomolib/prep/stripe.py +++ b/src/httomolib/prep/stripe.py @@ -23,7 +23,8 @@ import cupy as cp import numpy as np from cupy import abs, mean, ndarray -from cupyx.scipy.ndimage import median_filter, binary_dilation +from cupyx.scipy.ndimage import median_filter, binary_dilation, uniform_filter1d +from ._rgi import interpn __all__ = [ 'detect_stripes', @@ -32,7 +33,7 @@ 'remove_stripes_titarenko_cupy', ] -# TODO: port 'remove_all_stripe' and 'remove_dead_stripe' from +# TODO: port 'remove_all_stripe' from # https://github.com/tomopy/tomopy/blob/master/source/tomopy/prep/stripe.py @@ -148,6 +149,80 @@ def _detect_stripe(listdata: ndarray, snr: float) -> ndarray: return listmask +## %%%%%%%%%%%%%%%%%%%%% remove_dead_stripe_cupy %%%%%%%%%%%%%%%%%%%%%%%%% ## +## Naive CuPy port of the NumPy implementation in TomoPy +def remove_dead_stripe_cupy(tomo: ndarray, snr: float=3, size: int=51, + norm: bool=True) -> ndarray: + """ + Remove unresponsive and fluctuating stripe artifacts from sinogram using + Nghia Vo's approach :cite:`Vo:18` (algorithm 6). + + Parameters + ---------- + tomo : ndarray + 3D tomographic data. + + snr : float + Ratio used to detect locations of large stripes. + Greater is less sensitive. + + size : int + Window size of the median filter. + + norm : bool, optional + Remove residual stripes if True. + + Returns + ------- + ndarray + Corrected 3D tomographic data. + """ + matindex = _create_matindex(tomo.shape[2], tomo.shape[0]) + for m in range(tomo.shape[1]): + sino = tomo[:, m, :] + tomo[:, m, :] = _rs_dead(sino, snr, size, matindex, norm) + return tomo + + +def _rs_dead(sinogram, snr, size, matindex, norm=True): + """ + Remove unresponsive and fluctuating stripes. + """ + (nrow, _) = sinogram.shape + sinosmooth = cp.apply_along_axis(uniform_filter1d, 0, sinogram, 10) + listdiff = cp.sum(abs(sinogram - sinosmooth), axis=0) + listdiffbck = median_filter(listdiff, size) + # TODO: Same situation as the analagous part in `_rs_large()`, see that + # function's comment when using `cp.divide()`. + listfact = cp.divide(listdiff, listdiffbck, + #out=np.ones_like(listdiff), + #where=listdiffbck != 0 + ) + listmask = _detect_stripe(listfact, snr) + listmask = binary_dilation(listmask, iterations=1).astype(listmask.dtype) + listmask[0:2] = 0.0 + listmask[-2:] = 0.0 + # TODO: Same situation as the analagous part in `_rs_large()`, see that + # function's comment when using `cp.where()`. + listx = cp.where(listmask < 1.0)[0] + listy = cp.arange(nrow) + matz = sinogram[:, listx] + points = (listy, listx) + # Uses N-dimensional interpolation function copied from CuPy source code of + # v12.0.0b3 pre-release + finter = interpn(points, matz, tuple(cp.meshgrid(listy, listx)), + method='linear') + # TODO: Same situation as the analagous part in `_rs_large()`, see that + # function's comment when using `cp.where()`. + listxmiss = cp.where(listmask > 0.0)[0] + if len(listxmiss) > 0: + sinogram[:, listxmiss] = finter(listxmiss, listy) + # Remove residual stripes + if norm is True: + sinogram = _rs_large(sinogram, snr, size, matindex) + return sinogram + + ## %%%%%%%%%%%%%%%%% remove_stripe_based_sorting_cupy %%%%%%%%%%%%%%%%%%%%% ## ## Naive CuPy port of the NumPy implementation in TomoPy def remove_stripe_based_sorting_cupy( From 0f4e060b4510efd42450aa199cf4217475d21896 Mon Sep 17 00:00:00 2001 From: Yousef Moazzam Date: Fri, 20 Jan 2023 10:52:52 +0000 Subject: [PATCH 3/3] Port all stripe removal from TomoPy to CuPy This method is simply running different stripe removal methods in sequence: - removal of unresponsive and fluctuating stripes, via `_rs_dead()` - removal of full and partial stripes, via `_rs_sort()` and is doing nothing new. Therefore, the non-zero residual between the TomoPy/CPU and CuPy/GPU implementations of this particular method is due to the non-zero residuals present for the dependent stripe removal methods, rather than some new processing being introduced by `remove_all_stripe_cupy()` which adds to the residual. --- src/httomolib/prep/stripe.py | 40 ++++++++++++++++++++++++++++++++++-- 1 file changed, 38 insertions(+), 2 deletions(-) diff --git a/src/httomolib/prep/stripe.py b/src/httomolib/prep/stripe.py index 893923b4..06ca77e7 100644 --- a/src/httomolib/prep/stripe.py +++ b/src/httomolib/prep/stripe.py @@ -33,8 +33,44 @@ 'remove_stripes_titarenko_cupy', ] -# TODO: port 'remove_all_stripe' from -# https://github.com/tomopy/tomopy/blob/master/source/tomopy/prep/stripe.py +## %%%%%%%%%%%%%%%%%%%%% remove_all_stripe_cupy %%%%%%%%%%%%%%%%%%%%%%%%% ## +## Naive CuPy port of the NumPy implementation in TomoPy +def remove_all_stripe_cupy(tomo: ndarray, snr: float=3, la_size: int=61, + sm_size: int=21, dim: int=1): + """ + Remove all types of stripe artifacts from sinogram using Nghia Vo's + approach :cite:`Vo:18` (combination of algorithm 3,4,5, and 6). + + Parameters + ---------- + tomo : ndarray + 3D tomographic data. + + snr : float + Ratio used to locate large stripes. + Greater is less sensitive. + + la_size : int + Window size of the median filter to remove large stripes. + + sm_size : int + Window size of the median filter to remove small-to-medium stripes. + + dim : {1, 2}, optional + Dimension of the window. + + Returns + ------- + ndarray + Corrected 3D tomographic data. + """ + matindex = _create_matindex(tomo.shape[2], tomo.shape[0]) + for m in range(tomo.shape[1]): + sino = tomo[:, m, :] + sino = _rs_dead(sino, snr, la_size, matindex) + sino = _rs_sort(sino, sm_size, matindex, dim) + tomo[:, m, :] = sino + return tomo ## %%%%%%%%%%%%%%%%%%%%% remove_large_stripe_cupy %%%%%%%%%%%%%%%%%%%%%%%%% ##