Source code for algotom.prep.removal

# ============================================================================
# ============================================================================
# Copyright (c) 2021 Nghia T. Vo. All rights reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ============================================================================
# Author: Nghia T. Vo
# E-mail:  
# Description: Python module of artifact removal techniques.
# Contributors:
# ============================================================================

"""
Module of removal methods in the preprocessing stage:

    -   Many methods for removing stripe artifact in a sinogram
        (<-> ring artifact in a reconstructed image).
    -   A zinger removal method.
    -   Blob removal methods.
"""

import numpy as np
import scipy.ndimage as ndi
from scipy import interpolate
import numpy.fft as fft
import algotom.util.utility as util


[docs]def remove_stripe_based_sorting(sinogram, size=21, dim=1, **options): """ Remove stripe artifacts in a sinogram using the sorting technique, algorithm 3 in Ref. [1]. Angular direction is along the axis 0. Parameters ---------- sinogram : array_like 2D array. Sinogram image. size : int Window size of the median filter. dim : {1, 2}, optional Dimension of the window. options : dict, optional Use another smoothing filter rather than the median filter. E.g. options={"method": "gaussian_filter", "para1": (1,21)} Returns ------- array_like 2D array. Stripe-removed sinogram. References ---------- [1] : https://doi.org/10.1364/OE.26.028396 """ msg = "\n Please use the dictionary format: options={'method':" \ " 'filter_name', 'para1': parameter_1, 'para2': parameter_2}" sino_sort, sino_index = util.sort_forward(np.float32(sinogram), axis=0) if len(options) == 0: if dim == 2: sino_sort = ndi.median_filter(sino_sort, (size, size)) else: sino_sort = ndi.median_filter(sino_sort, (1, size)) else: for opt_name in options: opt = options[opt_name] if not isinstance(opt, dict): raise ValueError(msg) method = tuple(opt.values())[0] para = tuple(opt.values())[1:] if method in dir(ndi): try: sino_sort = getattr(ndi, method)(sino_sort, *para) except Exception: raise ValueError(msg) else: if method in dir(util): try: sino_sort = getattr(util, method)(sino_sort, *para) except Exception: raise ValueError(msg) else: raise ValueError("Can't find the method: '{}' in the" " namespace".format(method)) return util.sort_backward(sino_sort, sino_index, axis=0)
[docs]def remove_stripe_based_filtering(sinogram, sigma=3, size=21, dim=1, sort=True, **options): """ Remove stripe artifacts in a sinogram using the filtering technique, algorithm 2 in Ref. [1]. Angular direction is along the axis 0. Parameters ---------- sinogram : array_like 2D array. Sinogram image sigma : int Sigma of the Gaussian window used to separate the low-pass and high-pass components of the intensity profile of each column. size : int Window size of the median filter. dim : {1, 2}, optional Dimension of the window. sort : bool, optional Apply sorting if True. options : dict, optional Use another smoothing filter rather than the median filter. E.g. options={"method": "gaussian_filter", "para1": (1,21))}. Returns ------- array_like 2D array. Stripe-removed sinogram. References ---------- [1] : https://doi.org/10.1364/OE.26.028396 """ msg = "\n Please use the dictionary format: options={'method':" \ " 'filter_name', 'para1': parameter_1, 'para2': parameter_2}" window = {"name": "gaussian", "sigma": sigma} sino_smooth, sino_sharp = util.separate_frequency_component( np.float32(sinogram), axis=0, window=window) sino_index = None if sort is True: sino_smooth, sino_index = util.sort_forward(sino_smooth, axis=0) if len(options) == 0: if dim == 2: sino_smooth = ndi.median_filter(sino_smooth, (size, size)) else: sino_smooth = ndi.median_filter(sino_smooth, (1, size)) else: for opt_name in options: opt = options[opt_name] if not isinstance(opt, dict): raise ValueError(msg) method = tuple(opt.values())[0] para = tuple(opt.values())[1:] if method in dir(ndi): try: sino_smooth = getattr(ndi, method)(sino_smooth, *para) except Exception: raise ValueError(msg) else: if method in dir(util): try: sino_smooth = getattr(util, method)(sino_smooth, *para) except Exception: raise ValueError(msg) else: raise ValueError("Can't find the method: '{}' in the" " namespace".format(method)) if sort is True: sino_smooth = util.sort_backward(sino_smooth, sino_index, axis=0) return sino_smooth + sino_sharp
[docs]def remove_stripe_based_fitting(sinogram, order=2, sigma=10, sort=False, num_chunk=1, **options): """ Remove stripe artifacts in a sinogram using the fitting technique, algorithm 1 in Ref. [1]. Angular direction is along the axis 0. Parameters ---------- sinogram : array_like 2D array. Sinogram image order : int Polynomial fit order. sigma : int Sigma of the Gaussian window in the x-direction. Smaller is stronger. sort : bool, optional Apply sorting if True. num_chunk : int Number of chunks of rows to apply the fitting. options : dict, optional Use another smoothing filter rather than the Fourier gaussian filter. E.g. options={"method": "gaussian_filter", "para1": (1,21))}. Returns ------- array_like 2D array. Stripe-removed sinogram. References ---------- [1] : https://doi.org/10.1364/OE.26.028396 """ msg = "\n Please use the dictionary format: options={'method':" \ " 'filter_name', 'para1': parameter_1, 'para2': parameter_2}" (nrow, ncol) = sinogram.shape pad = min(150, int(0.1 * nrow)) sigmay = int(np.clip(min(60, int(0.1 * ncol)), 10, None)) sino_index = None if sort is True: sinogram, sino_index = util.sort_forward(sinogram, axis=0) sino_fit = util.generate_fitted_image(sinogram, order, axis=0, num_chunk=num_chunk) if len(options) == 0: sino_filt = util.apply_gaussian_filter(sino_fit, sigma, sigmay, pad) else: sino_filt = np.copy(sino_fit) for opt_name in options: opt = options[opt_name] if not isinstance(opt, dict): raise ValueError(msg) method = tuple(opt.values())[0] para = tuple(opt.values())[1:] if method in dir(ndi): try: sino_filt = getattr(ndi, method)(sino_filt, *para) except Exception: raise ValueError(msg) else: if method in dir(util): try: sino_filt = getattr(util, method)(sino_filt, *para) except Exception: raise ValueError(msg) else: raise ValueError("Can't find the method: '{}' in the" " namespace".format(method)) sino_filt = np.mean(np.abs(sino_fit)) * sino_filt / np.mean( np.abs(sino_filt)) sino_corr = ((sinogram / sino_fit) * sino_filt) if sort is True: sino_corr = util.sort_backward(sino_corr, sino_index, axis=0) return sino_corr
[docs]def remove_large_stripe(sinogram, snr=3.0, size=51, drop_ratio=0.1, norm=True, **options): """ Remove large stripe artifacts in a sinogram, algorithm 5 in Ref. [1]. Angular direction is along the axis 0. Parameters ---------- sinogram : array_like 2D array. Sinogram image snr : float Ratio (>1.0) for stripe detection. Greater is less sensitive. size : int Window size of the median filter. drop_ratio : float, optional Ratio of pixels to be dropped, which is used to reduce the possibility of the false detection of stripes. norm : bool, optional Apply normalization if True. options : dict, optional Use another smoothing filter rather than the median filter. E.g. options={"method": "gaussian_filter", "para1": (1,21))}. Returns ------- array_like 2D array. Stripe-removed sinogram. References ---------- [1] : https://doi.org/10.1364/OE.26.028396 """ msg = "\n Please use the dictionary format: options={'method':" \ " 'filter_name', 'para1': parameter_1, 'para2': parameter_2}" sinogram = np.copy(np.float32(sinogram)) drop_ratio = np.clip(drop_ratio, 0.0, 0.8) (nrow, ncol) = sinogram.shape ndrop = int(0.5 * drop_ratio * nrow) sino_sort, sino_index = util.sort_forward(sinogram, axis=0) if len(options) == 0: sino_smooth = ndi.median_filter(sino_sort, (1, size)) else: sino_smooth = np.copy(sino_sort) for opt_name in options: opt = options[opt_name] if not isinstance(opt, dict): raise ValueError(msg) method = tuple(opt.values())[0] para = tuple(opt.values())[1:] if method in dir(ndi): try: sino_smooth = getattr(ndi, method)(sino_smooth, *para) except Exception: raise ValueError(msg) else: if method in dir(util): try: sino_smooth = getattr(util, method)(sino_smooth, *para) except Exception: raise ValueError(msg) else: raise ValueError("Can't find the method: '{}' in the" " namespace".format(method)) list1 = np.mean(sino_sort[ndrop:nrow - ndrop], axis=0) list2 = np.mean(sino_smooth[ndrop:nrow - ndrop], axis=0) list_fact = np.divide(list1, list2, out=np.ones_like(list1), where=list2 != 0) list_mask = util.detect_stripe(list_fact, snr) list_mask = np.float32(ndi.binary_dilation(list_mask, iterations=1)) if norm is True: sinogram = sinogram / np.tile(list_fact, (nrow, 1)) sino_corr = util.sort_backward(sino_smooth, sino_index, axis=0) xlist_miss = np.where(list_mask > 0.0)[0] sinogram[:, xlist_miss] = sino_corr[:, xlist_miss] return sinogram
[docs]def remove_dead_stripe(sinogram, snr=3.0, size=51, residual=True, smooth_strength=10): """ Remove unresponsive or fluctuating stripe artifacts in a sinogram, algorithm 6 in Ref. [1]. Angular direction is along the axis 0. Parameters ---------- sinogram : array_like 2D array. Sinogram image. snr : float Ratio (>1.0) for stripe detection. Greater is less sensitive. size : int Window size of the median filter. residual : bool, optional Removing residual stripes if True. smooth_strength : int, optional Window size of the uniform filter used to detect stripes. Returns ------- ndarray 2D array. Stripe-removed sinogram. References ---------- [1] : https://doi.org/10.1364/OE.26.028396 """ sinogram = np.copy(sinogram) (nrow, ncol) = sinogram.shape sino_smooth = np.apply_along_axis(ndi.uniform_filter1d, 0, sinogram, smooth_strength) list_diff = np.sum(np.abs(sinogram - sino_smooth), axis=0) list_diff_bck = ndi.median_filter(list_diff, size) nmean = np.mean(np.abs(list_diff_bck)) list_diff_bck[list_diff_bck == 0.0] = nmean list_fact = list_diff / list_diff_bck list_mask = util.detect_stripe(list_fact, snr) list_mask = np.float32(ndi.binary_dilation(list_mask, iterations=1)) list_mask[0:2] = 0.0 list_mask[-2:] = 0.0 xlist = np.where(list_mask < 1.0)[0] ylist = np.arange(nrow) finter = interpolate.RectBivariateSpline(ylist, xlist, sinogram[:, xlist], kx=1, ky=1) xlist_miss = np.where(list_mask > 0.0)[0] if (ncol // 3) > len(xlist_miss) > 0: x_mat_miss, y_mat = np.meshgrid(xlist_miss, ylist) output = finter.ev(np.ndarray.flatten(y_mat), np.ndarray.flatten(x_mat_miss)) sinogram[:, xlist_miss] = output.reshape(x_mat_miss.shape) if residual is True: sinogram = remove_large_stripe(sinogram, snr, size) return sinogram
[docs]def remove_all_stripe(sinogram, snr=3.0, la_size=51, sm_size=21, drop_ratio=0.1, dim=1, **options): """ Remove all types of stripe artifacts in a sinogram by combining algorithm 6, 5, 4, and 3 in Ref. [1]. Angular direction is along the axis 0. Parameters ---------- sinogram : array_like 2D array. Sinogram image. snr : float Ratio (>1.0) for stripe detection. 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. drop_ratio : float, optional Ratio of pixels to be dropped, which is used to reduce the possibility of the false detection of stripes. dim : {1, 2}, optional Dimension of the window. options : dict, optional Use another smoothing filter rather than the median filter. E.g. options={"method": "gaussian_filter", "para1": (1,21))} Returns ------- array_like 2D array. Stripe-removed sinogram. References ---------- [1] : https://doi.org/10.1364/OE.26.028396 """ sinogram = remove_dead_stripe(sinogram, snr, la_size, residual=False) sinogram = remove_large_stripe(sinogram, snr, la_size, drop_ratio, **options) sinogram = remove_stripe_based_sorting(sinogram, sm_size, dim, **options) return sinogram
[docs]def remove_stripe_based_2d_filtering_sorting(sinogram, sigma=3, size=21, dim=1, **options): """ Remove stripes using a 2D low-pass filter and the sorting-based technique, algorithm in section 3.3.4 in Ref. [1]. Angular direction is along the axis 0. Parameters --------- sinogram : array_like 2D array. Sinogram image. sigma : int Sigma of the Gaussian window. size : int Window size of the median filter. dim : {1, 2}, optional Dimension of the window. Returns ------- array_like 2D array. Stripe-removed sinogram. References ---------- [1] : https://doi.org/10.1117/12.2530324 """ (nrow, ncol) = sinogram.shape pad = min(150, int(0.1 * min(nrow, ncol))) sino_smooth = util.apply_gaussian_filter(sinogram, sigma, sigma, pad) sino_sharp = sinogram - sino_smooth sino_sharp = remove_stripe_based_sorting(sino_sharp, size, dim, **options) return sino_smooth + sino_sharp
[docs]def remove_stripe_based_normalization(sinogram, sigma=15, num_chunk=1, sort=True, **options): """ Remove stripes using the method in Ref. [1]. Angular direction is along the axis 0. Parameters ---------- sinogram : array_like 2D array. Sinogram image. sigma : int Sigma of the Gaussian window. num_chunk : int Number of chunks of rows. sort : bool, optional Apply sorting (Ref. [2]) if True. options : dict, optional Use another smoothing 1D-filter rather than the Gaussian filter. E.g. options={"method": "median_filter", "para1": 21)}. Returns ------- array_like 2D array. Stripe-removed sinogram. References ---------- [1] : https://www.mcs.anl.gov/research/projects/X-ray-cmt/rivers/tutorial.html [2] : https://doi.org/10.1364/OE.26.028396 """ msg = "\n Please use the dictionary format: options={'method':" \ " 'filter_name', 'para1': parameter_1, 'para2': parameter_2}" \ "\n Note that the filter must be a 1D-filter." (nrow, _) = sinogram.shape sinogram = np.copy(sinogram) sino_index = None if sort is True: sinogram, sino_index = util.sort_forward(sinogram, axis=0) list_index = np.array_split(np.arange(nrow), num_chunk) for pos in list_index: bindex = pos[0] eindex = pos[-1] + 1 list_mean = np.mean(sinogram[bindex:eindex], axis=0) if len(options) == 0: list_filt = ndi.gaussian_filter(list_mean, sigma) else: list_filt = np.copy(list_mean) for opt_name in options: opt = options[opt_name] if not isinstance(opt, dict): raise ValueError(msg) method = tuple(opt.values())[0] para = tuple(opt.values())[1:] if method in dir(ndi): try: list_filt = getattr(ndi, method)(list_filt, *para) except Exception: raise ValueError(msg) else: if method in dir(util): try: list_filt = getattr(util, method)(list_filt, *para) except Exception: raise ValueError(msg) else: raise ValueError("Can't find the method: '{}' in the" " namespace".format(method)) list_coe = list_filt - list_mean matcoe = np.tile(list_coe, (eindex - bindex, 1)) sinogram[bindex:eindex, :] = sinogram[bindex:eindex, :] + matcoe if sort is True: sinogram = util.sort_backward(sinogram, sino_index, axis=0) return sinogram
[docs]def remove_stripe_based_regularization(sinogram, alpha=0.0005, num_chunk=1, apply_log=True, sort=True): """ Remove stripes using the method in Ref. [1]. Angular direction is along the axis 0. Parameters ---------- sinogram : array_like 2D array. Sinogram image. alpha : float Regularization parameter, e.g. 0.0005. Smaller is stronger. num_chunk : int Number of chunks of rows. apply_log : bool Apply the logarithm function to the sinogram if True. sort : bool, optional Apply sorting (Ref. [2]) if True. Returns ------- array_like 2D array. Stripe-removed sinogram. References ---------- [1] : https://doi.org/10.1016/j.aml.2010.08.022 [2] : https://doi.org/10.1364/OE.26.028396 """ (nrow, ncol) = sinogram.shape sinogram = np.copy(sinogram) sino_index = None if sort is True: sinogram, sino_index = util.sort_forward(sinogram, axis=0) if apply_log is True: if np.any(sinogram <= 0.0): nmean = np.mean(np.abs(sinogram)) sinogram[sinogram <= 0.0] = nmean sinogram = -np.log(sinogram) else: sinogram = -np.log(sinogram) sijmat = util.calculate_regularization_coefficient(ncol, alpha) list_index = np.array_split(np.arange(nrow), num_chunk) list_grad = np.zeros(ncol, dtype=np.float32) mat_grad = np.zeros((ncol, ncol), dtype=np.float32) for pos in list_index: bindex = pos[0] eindex = pos[-1] + 1 list_mean = np.mean(sinogram[bindex:eindex], axis=0) list_grad[1:-1] = (-1) * np.diff(list_mean, 2) list_grad[0] = list_mean[0] - list_mean[1] list_grad[-1] = list_mean[-1] - list_mean[-2] mat_grad[:] = list_grad list_coe = np.sum(mat_grad * sijmat, axis=1) mat_coe = np.tile(list_coe, (eindex - bindex, 1)) sinogram[bindex:eindex, :] = sinogram[bindex:eindex, :] + mat_coe if sort is True: sinogram = util.sort_backward(sinogram, sino_index, axis=0) if apply_log is True: sinogram = np.exp(-sinogram) return sinogram
[docs]def remove_stripe_based_fft(sinogram, u=20, n=8, v=1, sort=False): """ Remove stripes using the method in Ref. [1]. Angular direction is along the axis 0. Parameters ---------- sinogram : array_like 2D array. Sinogram image. u : int Cutoff frequency. n : int Filter order. v : int Number of rows (* 2) to be applied the filter. sort : bool, optional Apply sorting (Ref. [2]) if True. Returns ------- ndarray 2D array. Stripe-removed sinogram. References ---------- [1] : https://doi.org/10.1063/1.1149043 [2] : https://doi.org/10.1364/OE.26.028396 """ sino_index = None if sort is True: sinogram, sino_index = util.sort_forward(sinogram, axis=0) pad = min(150, int(0.1 * np.min(sinogram.shape))) sinogram = np.pad(sinogram, ((pad, pad), (0, 0)), mode='mean') sinogram = np.pad(sinogram, ((0, 0), (pad, pad)), mode='edge') (nrow, ncol) = sinogram.shape window_2d = util.make_2d_butterworth_window(ncol, nrow, u, v, n) sinogram = fft.ifft2( np.fft.ifftshift(np.fft.fftshift(fft.fft2(sinogram)) * window_2d)) sinogram = np.abs(sinogram[pad:nrow - pad, pad:ncol - pad]) if sort is True: sinogram = util.sort_backward(sinogram, sino_index, axis=0) return sinogram
[docs]def remove_stripe_based_wavelet_fft(sinogram, level=5, size=1, wavelet_name="db9", window_name="gaussian", sort=False, **options): """ Remove stripes using the method in Ref. [1]. Angular direction is along the axis 0. Parameters ---------- sinogram : array_like 2D array. Sinogram image. level : int Wavelet decomposition level. size : int Damping parameter. Larger is stronger. wavelet_name : str Name of a wavelet. Search pywavelets API for a full list. window_name : str High-pass window. Two options: "gaussian" or "butter". sort : bool, optional Apply sorting (Ref. [2]) if True. options : dict, optional Use another smoothing filter rather than the fft-gaussian-filter. E.g. options={"method": "gaussian_filter", "para1": (1,11))} Returns ------- array_like 2D array. Stripe-removed sinogram. References ---------- [1] : https://doi.org/10.1364/OE.17.008567 [2] : https://doi.org/10.1364/OE.26.028396 """ msg = "\n Please use the dictionary format: options={'method':" \ " 'filter_name', 'para1': parameter_1, 'para2': parameter_2}" sino_index = None if sort is True: sinogram, sino_index = util.sort_forward(sinogram, axis=0) (nrow, ncol) = sinogram.shape pad = min(150, int(0.1 * min(nrow, ncol))) sinogram = np.pad(sinogram, ((pad, pad), (0, 0)), mode='mean') sinogram = np.pad(sinogram, ((0, 0), (pad, pad)), mode='edge') output_data = util.apply_wavelet_decomposition(sinogram, wavelet_name, level=level) output_data = [list(data) for data in output_data] n_level = len(output_data[1:]) for i in range(1, n_level + 1): if len(options) == 0: (height, width) = output_data[i][1].shape window = np.transpose( util.make_2d_damping_window(height, width, size, window_name)) output_data[i][1] = np.real(np.fft.ifft2(np.fft.ifftshift( np.fft.fftshift(np.fft.fft2(output_data[i][1])) * window))) else: mat_smooth = np.copy(output_data[i][1]) for opt_name in options: opt = options[opt_name] if not isinstance(opt, dict): raise ValueError(msg) method = tuple(opt.values())[0] para = tuple(opt.values())[1:] if method in dir(ndi): try: mat_smooth = getattr(ndi, method)(mat_smooth, *para) except Exception: raise ValueError(msg) elif method in dir(util): try: mat_smooth = getattr(util, method)(mat_smooth, *para) except Exception: raise ValueError(msg) else: raise ValueError("Can't find the method: '{}' in the" " namespace".format(method)) output_data[i][1] = mat_smooth sinogram = util.apply_wavelet_reconstruction(output_data, wavelet_name) sinogram = sinogram[pad:nrow + pad, pad:ncol + pad] if sort is True: sinogram = util.sort_backward(sinogram, sino_index, axis=0) return sinogram
[docs]def remove_stripe_based_interpolation(sinogram, snr=3.0, size=51, drop_ratio=0.1, norm=True, kind="linear", **options): """ Combination of algorithm 4, 5, and 6 in Ref. [1]. Angular direction is along the axis 0. Parameters ---------- sinogram : array_like 2D array. Sinogram image snr : float Ratio (>1.0) for stripe detection. Greater is less sensitive. size : int Window size of the median filter used to detect stripes. drop_ratio : float, optional Ratio of pixels to be dropped, which is used to reduce the possibility of the false detection of stripes. norm : bool, optional Apply normalization if True. kind : {'linear', 'cubic', 'quintic'}, optional The kind of spline interpolation to use. Default is 'linear'. options : dict, optional Use another smoothing filter rather than the median filter. E.g. options={"method": "gaussian_filter", "para1": (1,21))} Returns ------- array_like 2D array. Stripe-removed sinogram. References ---------- [1] : https://doi.org/10.1364/OE.26.028396 """ msg = "\n Please use the dictionary format: options={'method':" \ " 'filter_name', 'para1': parameter_1, 'para2': parameter_2}" drop_ratio = np.clip(drop_ratio, 0.0, 0.8) sinogram = np.copy(sinogram) (nrow, ncol) = sinogram.shape ndrop = int(0.5 * drop_ratio * nrow) sino_sort = np.sort(sinogram, axis=0) if len(options) == 0: sino_smooth = ndi.median_filter(sino_sort, (1, size)) else: sino_smooth = np.copy(sino_sort) for opt_name in options: opt = options[opt_name] if not isinstance(opt, dict): raise ValueError(msg) method = tuple(opt.values())[0] para = tuple(opt.values())[1:] if method in dir(ndi): try: sino_smooth = getattr(ndi, method)(sino_smooth, *para) except Exception: raise ValueError(msg) else: if method in dir(util): try: sino_smooth = getattr(util, method)(sino_smooth, *para) except Exception: raise ValueError(msg) else: raise ValueError("Can't find the method: '{}' in the" " namespace".format(method)) list1 = np.mean(sino_sort[ndrop:nrow - ndrop], axis=0) list2 = np.mean(sino_smooth[ndrop:nrow - ndrop], axis=0) list1[list1 == 0.0] = np.float32(1.0) list_fact = np.divide(list1, list2, out=np.ones_like(list1), where=list2 != 0) list_mask = util.detect_stripe(list_fact, snr) list_mask = np.float32(ndi.binary_dilation(list_mask, iterations=1)) mat_fact = np.tile(list_fact, (nrow, 1)) if norm is True: sinogram = sinogram / mat_fact list_mask[0:2] = 0.0 list_mask[-2:] = 0.0 xlist = np.where(list_mask < 1.0)[0] ylist = np.arange(nrow) if kind == "cubic": finter = interpolate.RectBivariateSpline(ylist, xlist, sinogram[:, xlist], kx=2, ky=2) elif kind == "quintic": finter = interpolate.RectBivariateSpline(ylist, xlist, sinogram[:, xlist], kx=2, ky=2) else: finter = interpolate.RectBivariateSpline(ylist, xlist, sinogram[:, xlist], kx=1, ky=1) xlist_miss = np.where(list_mask > 0.0)[0] if len(xlist_miss) > 0: x_mat_miss, y_mat = np.meshgrid(xlist_miss, ylist) output = finter.ev(np.ndarray.flatten(y_mat), np.ndarray.flatten(x_mat_miss)) sinogram[:, xlist_miss] = output.reshape(x_mat_miss.shape) return sinogram
[docs]def check_zinger_size(mat, max_size): """ Check if the size of a zinger is smaller than a given size. Parameters ---------- mat : array_like 2D array. max_size : int Maximum size. Returns ------- bool """ check = False zinger_size = mat.sum() if zinger_size <= max_size: check = True return check
[docs]def select_zinger(mat, max_size): """ Select zingers smaller than a certain size. Parameters ---------- mat : array_like 2D array. max_size : int Maximum size in pixel. Returns ------- array_like 2D binary array. """ list_zin = ndi.find_objects(ndi.label(mat)[0]) zin_sel = [zin for zin in list_zin if check_zinger_size(mat[zin], max_size)] mat_out = np.zeros_like(mat) for _, j in enumerate(zin_sel): mat_out[j] = mat[j] return mat_out
[docs]def remove_zinger(mat, threshold, size=2, check_size=False): """ Remove zinger using the method in Ref. [1], working on a projection image or sinogram image. Parameters ---------- mat : array_like 2D array. Projection image or sinogram image. threshold : float Threshold to segment zingers. Smaller is more sensitive. Recommended range [0.05, 0.1]. size : int Size of a zinger. check_size : bool Enable/disable size checking before removal. Returns ------- array_like 2D array. Zinger-removed image. References ---------- [1] : https://doi.org/10.1364/OE.418448 """ step = np.clip(size, 1, None) mat = np.copy(mat) mat_ave = [] for i in range(-size, size + 1, step): for j in range(-size, size + 1, step): if (i != 0) or (j != 0): mat_ave.append(np.roll(np.roll(mat, i, axis=0), j, axis=1)) mat_ave = np.mean(np.asarray(mat_ave), axis=0) mat_ave[mat_ave == 0.0] = 1.0 mat_nor = mat / mat_ave - 1.0 mask = np.asarray(mat_nor > threshold, dtype=np.float32) if check_size: mask = select_zinger(mask, size) mat[mask > 0.0] = mat_ave[mask > 0.0] return mat
[docs]def generate_blob_mask(flat, size, snr): """ Generate a binary mask of blobs from a flat-field image (Ref. [1]). Parameters ---------- flat : array_like 2D array. Flat-field image. size : float Estimated size of the largest blob. snr : float Ratio used to segment blobs. Returns ------- array_like 2D array. Binary mask. References ---------- [1] : https://doi.org/10.1364/OE.418448 """ mat = ndi.median_filter(flat, (2, 2)) mask = np.zeros_like(mat) for i in range(mat.shape[0]): line = mat[i] line_fil = ndi.median_filter(line, size) line_fil[line_fil == 0.0] = np.mean(line_fil) line_norm = line / line_fil mask_1d = util.detect_stripe(line_norm, snr) mask_1d[0:2] = 0.0 mask_1d[-2:] = 0.0 mask[i] = np.float32(ndi.binary_dilation(mask_1d, iterations=1)) return mask
[docs]def remove_blob_1d(sino_1d, mask_1d): """ Remove blobs in one row of a sinogram, e.g. for a helical scan as shown in Ref. [1]. Parameters ---------- sino_1d : array_like 1D array. A row of a sinogram. mask_1d : array_like 1D binary mask. Returns ------- array_like 1D array. Notes ----- The method is used to remove streak artifacts caused by blobs in a sinogram generated from a helical-scan data [1]. References ---------- [1] : https://doi.org/10.1364/OE.418448 """ sino_1d = np.copy(sino_1d) xlist = np.where(mask_1d < 1.0)[0] ylist = sino_1d[xlist] finter = interpolate.interp1d(xlist, ylist) mask_1d[:2] = 0.0 mask_1d[-2:] = 0.0 xlist_miss = np.where(mask_1d > 0.0)[0] if len(xlist_miss) > 0: sino_1d[xlist_miss] = finter(xlist_miss) return sino_1d
[docs]def remove_blob(mat, mask): """ Remove blobs in an image. Parameters ---------- mat : array_like 2D array. Projection image or sinogram image. mask : array_like 2D binary mask. Returns ------- array_like 2D array. """ mat = np.copy(mat) if mat.shape != mask.shape: raise ValueError("The image and the mask are not the same shape !!!") for i in range(mat.shape[0]): array_1d = mat[i] mask_1d = mask[i] mask_1d[:2] = 0.0 mask_1d[-2:] = 0.0 xlist = np.where(mask_1d < 1.0)[0] ylist = array_1d[xlist] finter = interpolate.interp1d(xlist, ylist) xlist_miss = np.where(mask_1d > 0.0)[0] if len(xlist_miss) > 0: array_1d[xlist_miss] = finter(xlist_miss) mat[i] = array_1d return mat