#!/usr/bin/env python
# -*- coding: utf-8 -*-
# #########################################################################
# Copyright (c) 2015-2019, UChicago Argonne, LLC. All rights reserved. #
# #
# Copyright 2015-2019. UChicago Argonne, LLC. This software was produced #
# under U.S. Government contract DE-AC02-06CH11357 for Argonne National #
# Laboratory (ANL), which is operated by UChicago Argonne, LLC for the #
# U.S. Department of Energy. The U.S. Government has rights to use, #
# reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR #
# UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR #
# ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is #
# modified to produce derivative works, such modified software should #
# be clearly marked, so as not to confuse it with the version available #
# from ANL. #
# #
# Additionally, redistribution and use in source and binary forms, with #
# or without modification, are permitted provided that the following #
# conditions are met: #
# #
# * Redistributions of source code must retain the above copyright #
# notice, this list of conditions and the following disclaimer. #
# #
# * Redistributions in binary form must reproduce the above copyright #
# notice, this list of conditions and the following disclaimer in #
# the documentation and/or other materials provided with the #
# distribution. #
# #
# * Neither the name of UChicago Argonne, LLC, Argonne National #
# Laboratory, ANL, the U.S. Government, nor the names of its #
# contributors may be used to endorse or promote products derived #
# from this software without specific prior written permission. #
# #
# THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS #
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT #
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS #
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago #
# Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, #
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, #
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; #
# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER #
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT #
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN #
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE #
# POSSIBILITY OF SUCH DAMAGE. #
# #########################################################################
"""
Module for pre-processing tasks.
"""
from __future__ import (absolute_import, division, print_function,
unicode_literals)
import numpy as np
import pywt
import tomopy.util.extern as extern
import tomopy.util.mproc as mproc
import tomopy.util.dtype as dtype
from tomopy.util.misc import (fft, ifft, fft2, ifft2)
from scipy.ndimage import median_filter
from scipy import signal
from scipy.signal import savgol_filter
from scipy.ndimage import binary_dilation
from scipy.ndimage import uniform_filter1d
from scipy import interpolate
import logging
logger = logging.getLogger(__name__)
__author__ = "Doga Gursoy, Eduardo X. Miqueles, Nghia Vo"
__credits__ = "Juan V. Bermudez, Hugo H. Slepicka"
__copyright__ = "Copyright (c) 2015, UChicago Argonne, LLC."
__docformat__ = 'restructuredtext en'
__all__ = ['remove_stripe_fw',
'remove_stripe_ti',
'remove_stripe_sf',
'remove_stripe_based_sorting',
'remove_stripe_based_filtering',
'remove_stripe_based_fitting',
'remove_large_stripe',
'remove_dead_stripe',
'remove_all_stripe']
[docs]def remove_stripe_fw(
tomo, level=None, wname='db5', sigma=2,
pad=True, ncore=None, nchunk=None):
"""
Remove horizontal stripes from sinogram using the Fourier-Wavelet (FW)
based method :cite:`Munch:09`.
Parameters
----------
tomo : ndarray
3D tomographic data.
level : int, optional
Number of discrete wavelet transform levels.
wname : str, optional
Type of the wavelet filter. 'haar', 'db5', sym5', etc.
sigma : float, optional
Damping parameter in Fourier space.
pad : bool, optional
If True, extend the size of the sinogram by padding with zeros.
ncore : int, optional
Number of cores that will be assigned to jobs.
nchunk : int, optional
Chunk size for each core.
Returns
-------
ndarray
Corrected 3D tomographic data.
"""
if level is None:
size = np.max(tomo.shape)
level = int(np.ceil(np.log2(size)))
arr = mproc.distribute_jobs(
tomo,
func=_remove_stripe_fw,
args=(level, wname, sigma, pad),
axis=1,
ncore=ncore,
nchunk=nchunk)
return arr
def _remove_stripe_fw(tomo, level, wname, sigma, pad):
dx, dy, dz = tomo.shape
nx = dx
if pad:
nx = dx + dx // 8
xshift = int((nx - dx) // 2)
num_jobs = tomo.shape[1]
for m in range(num_jobs):
sli = np.zeros((nx, dz), dtype='float32')
sli[xshift:dx + xshift] = tomo[:, m, :]
# Wavelet decomposition.
cH = []
cV = []
cD = []
for n in range(level):
sli, (cHt, cVt, cDt) = pywt.dwt2(sli, wname)
cH.append(cHt)
cV.append(cVt)
cD.append(cDt)
# FFT transform of horizontal frequency bands.
for n in range(level):
# FFT
fcV = np.fft.fftshift(fft(cV[n], axis=0, extra_info=num_jobs))
my, mx = fcV.shape
# Damping of ring artifact information.
y_hat = (np.arange(-my, my, 2, dtype='float32') + 1) / 2
damp = -np.expm1(-np.square(y_hat) / (2 * np.square(sigma)))
fcV *= np.transpose(np.tile(damp, (mx, 1)))
# Inverse FFT.
cV[n] = np.real(ifft(np.fft.ifftshift(
fcV), axis=0, extra_info=num_jobs))
# Wavelet reconstruction.
for n in range(level)[::-1]:
sli = sli[0:cH[n].shape[0], 0:cH[n].shape[1]]
sli = pywt.idwt2((sli, (cH[n], cV[n], cD[n])), wname)
tomo[:, m, :] = sli[xshift:dx + xshift, 0:dz]
[docs]def remove_stripe_ti(tomo, nblock=0, alpha=1.5, ncore=None, nchunk=None):
"""
Remove horizontal stripes from sinogram using Titarenko's
approach :cite:`Miqueles:14`.
Parameters
----------
tomo : ndarray
3D tomographic data.
nblock : int, optional
Number of blocks.
alpha : int, optional
Damping factor.
ncore : int, optional
Number of cores that will be assigned to jobs.
nchunk : int, optional
Chunk size for each core.
Returns
-------
ndarray
Corrected 3D tomographic data.
"""
arr = mproc.distribute_jobs(
tomo,
func=_remove_stripe_ti,
args=(nblock, alpha),
axis=1,
ncore=ncore,
nchunk=nchunk)
return arr
def _remove_stripe_ti(tomo, nblock, alpha):
for m in range(tomo.shape[1]):
sino = tomo[:, m, :]
if nblock == 0:
d1 = _ring(sino, 1, 1)
d2 = _ring(sino, 2, 1)
p = d1 * d2
d = np.sqrt(p + alpha * np.abs(p.min()))
else:
size = int(sino.shape[0] / nblock)
d1 = _ringb(sino, 1, 1, size)
d2 = _ringb(sino, 2, 1, size)
p = d1 * d2
d = np.sqrt(p + alpha * np.fabs(p.min()))
tomo[:, m, :] = d
def _kernel(m, n):
v = [[np.array([1, -1]),
np.array([-3 / 2, 2, -1 / 2]),
np.array([-11 / 6, 3, -3 / 2, 1 / 3])],
[np.array([-1, 2, -1]),
np.array([2, -5, 4, -1])],
[np.array([-1, 3, -3, 1])]]
return v[m - 1][n - 1]
def _ringMatXvec(h, x):
s = np.convolve(x, np.flipud(h))
u = s[np.size(h) - 1:np.size(x)]
y = np.convolve(u, h)
return y
def _ringCGM(h, alpha, f):
x0 = np.zeros(np.size(f))
r = f - (_ringMatXvec(h, x0) + alpha * x0)
w = -r
z = _ringMatXvec(h, w) + alpha * w
a = np.dot(r, w) / np.dot(w, z)
x = x0 + np.dot(a, w)
B = 0
for i in range(1000000):
r = r - np.dot(a, z)
if np.linalg.norm(r) < 0.0000001:
break
B = np.dot(r, z) / np.dot(w, z)
w = -r + np.dot(B, w)
z = _ringMatXvec(h, w) + alpha * w
a = np.dot(r, w) / np.dot(w, z)
x = x + np.dot(a, w)
return x
def _get_parameter(x):
return 1 / (2 * (x.sum(0).max() - x.sum(0).min()))
def _ring(sino, m, n):
mysino = np.transpose(sino)
R = np.size(mysino, 0)
N = np.size(mysino, 1)
# Remove NaN.
pos = np.where(np.isnan(mysino) is True)
mysino[pos] = 0
# Parameter.
alpha = _get_parameter(mysino)
# Mathematical correction.
pp = mysino.mean(1)
h = _kernel(m, n)
f = -_ringMatXvec(h, pp)
q = _ringCGM(h, alpha, f)
# Update sinogram.
q.shape = (R, 1)
K = np.kron(q, np.ones((1, N)))
new = np.add(mysino, K)
newsino = new.astype(np.float32)
return np.transpose(newsino)
def _ringb(sino, m, n, step):
mysino = np.transpose(sino)
R = np.size(mysino, 0)
N = np.size(mysino, 1)
# Remove NaN.
pos = np.where(np.isnan(mysino) is True)
mysino[pos] = 0
# Kernel & regularization parameter.
h = _kernel(m, n)
# Mathematical correction by blocks.
nblock = int(N / step)
new = np.ones((R, N))
for k in range(0, nblock):
sino_block = mysino[:, k * step:(k + 1) * step]
alpha = _get_parameter(sino_block)
pp = sino_block.mean(1)
f = -_ringMatXvec(h, pp)
q = _ringCGM(h, alpha, f)
# Update sinogram.
q.shape = (R, 1)
K = np.kron(q, np.ones((1, step)))
new[:, k * step:(k + 1) * step] = np.add(sino_block, K)
newsino = new.astype(np.float32)
return np.transpose(newsino)
[docs]def remove_stripe_sf(tomo, size=5, ncore=None, nchunk=None):
"""
Normalize raw projection data using a smoothing filter approach.
Parameters
----------
tomo : ndarray
3D tomographic data.
size : int, optional
Size of the smoothing filter.
ncore : int, optional
Number of cores that will be assigned to jobs.
nchunk : int, optional
Chunk size for each core.
Returns
-------
ndarray
Corrected 3D tomographic data.
"""
tomo = dtype.as_float32(tomo)
arr = mproc.distribute_jobs(
tomo,
func=extern.c_remove_stripe_sf,
args=(size,),
axis=1,
ncore=ncore,
nchunk=nchunk)
return arr
[docs]def remove_stripe_based_sorting(tomo, size=None, ncore=None, nchunk=None):
"""
Remove stripe artifacts from sinogram using Nghia Vo's
approach :cite:`Vo:18`
Algorithm 3 in the paper. Remove stripes using the sorting technique.
Work particularly well for removing partial stripes.
Parameters
----------
tomo : ndarray
3D tomographic data.
size : int
Window size of the median filter.
ncore : int, optional
Number of cores that will be assigned to jobs.
nchunk : int, optional
Chunk size for each core.
Returns
-------
ndarray
Corrected 3D tomographic data.
"""
arr = mproc.distribute_jobs(
tomo,
func=_remove_stripe_based_sorting,
args=(size,),
axis=1,
ncore=ncore,
nchunk=nchunk)
return arr
def _create_matindex(nrow, ncol):
"""
Create a 2D array of indexes used for sorting technique
"""
listindex = np.arange(0.0, ncol, 1.0)
matindex = np.tile(listindex, (nrow, 1))
return matindex
def _rs_sort(sinogram, size, matindex):
"""
Remove stripes using the sorting technique.
"""
sinogram = np.transpose(sinogram)
matcomb = np.asarray(np.dstack((matindex, sinogram)))
matsort = np.asarray(
[row[row[:, 1].argsort()] for row in matcomb])
matsort[:, :, 1] = median_filter(matsort[:, :, 1], (size, 1))
matsortback = np.asarray(
[row[row[:, 0].argsort()] for row in matsort])
sino_corrected = matsortback[:, :, 1]
return np.transpose(sino_corrected)
def _remove_stripe_based_sorting(tomo, size):
matindex = _create_matindex(tomo.shape[2], tomo.shape[0])
if size is None:
if tomo.shape[2] > 2000:
size = 21
else:
size = max(5, int(0.01 * tomo.shape[2]))
for m in range(tomo.shape[1]):
sino = tomo[:, m, :]
tomo[:, m, :] = _rs_sort(sino, size, matindex)
[docs]def remove_stripe_based_filtering(
tomo, sigma=3, size=None, ncore=None, nchunk=None):
"""
Remove stripe artifacts from sinogram using Nghia Vo's
approach :cite:`Vo:18`
Algorithm 2 in the paper. Remove stripes using the filtering technique.
Parameters
----------
tomo : ndarray
3D tomographic data.
sigma : float
Sigma of the Gaussian window which is used to separate
the low-pass and high-pass components of the intensity
profiles of each column.
size : int
Window size of the median filter.
ncore : int, optional
Number of cores that will be assigned to jobs.
nchunk : int, optional
Chunk size for each core.
Returns
-------
ndarray
Corrected 3D tomographic data.
"""
arr = mproc.distribute_jobs(
tomo,
func=_remove_stripe_based_filtering,
args=(sigma, size),
axis=1,
ncore=ncore,
nchunk=nchunk)
return arr
def _create_1d_window(ncol, sigma, pad):
window = signal.gaussian(ncol + 2 * pad, std=sigma)
return window
def _create_listsign(ncol, pad):
listsign = np.power(-1.0, np.arange(ncol + 2 * pad))
return listsign
def _rs_filter(sinogram, window, listsign, sigma, size, pad):
"""
Remove stripes using the filtering technique.
"""
sinogram = np.transpose(sinogram)
padded_sino = np.pad(sinogram, ((0, 0), (pad, pad)), mode='reflect')
(_, ncol) = padded_sino.shape
sinosmooth = np.zeros_like(sinogram)
for i, sinolist in enumerate(padded_sino):
sinosmooth[i] = np.real(
ifft(fft(sinolist * listsign) * window) * listsign)[pad:ncol - pad]
sinosharp = sinogram - sinosmooth
sinosmooth_cor = median_filter(sinosmooth, (size, 1))
return np.transpose(sinosmooth_cor + sinosharp)
def _remove_stripe_based_filtering(tomo, sigma, size):
pad = 150 # To reduce artifacts caused by FFT
window = _create_1d_window(tomo.shape[0], sigma, pad)
listsign = _create_listsign(tomo.shape[0], pad)
if size is None:
if tomo.shape[2] > 2000:
size = 21
else:
size = max(5, int(0.01 * tomo.shape[2]))
for m in range(tomo.shape[1]):
sino = tomo[:, m, :]
tomo[:, m, :] = _rs_filter(
sino, window, listsign, sigma, size, pad)
[docs]def remove_stripe_based_fitting(
tomo, order=3, sigma=(5, 20), ncore=None, nchunk=None):
"""
Remove horizontal stripes from sinogram using Nghia Vo's
approach :cite:`Vo:18`
Algorithm 1 in the paper. Remove stripes using the fitting technique.
Parameters
----------
tomo : ndarray
3D tomographic data.
order : int
Polynomial fit order.
sigma : tuple of 2 floats
Sigmas of a 2D Gaussian window in x and y direction.
ncore : int, optional
Number of cores that will be assigned to jobs.
nchunk : int, optional
Chunk size for each core.
Returns
-------
ndarray
Corrected 3D tomographic data.
"""
arr = mproc.distribute_jobs(
tomo,
func=_remove_stripe_based_fitting,
args=(order, sigma),
axis=1,
ncore=ncore,
nchunk=nchunk)
return arr
def _create_2d_window(nrow, ncol, sigma, pad):
"""
Create a 2D Gaussian window.
Parameters
----------
nrow : int
Height of the window.
ncol: int
Width of the window.
sigma: tuple of 2 floats
Sigmas of the window.
pad : int
Padding.
Returns
-------
2D array.
"""
(sigmax, sigmay) = sigma
height = nrow + 2 * pad
width = ncol + 2 * pad
centerx = (width - 1.0) / 2.0
centery = (height - 1.0) / 2.0
y, x = np.ogrid[-centery:height - centery, -centerx:width - centerx]
numx = 2.0 * sigmax * sigmax
numy = 2.0 * sigmay * sigmay
win2d = np.exp(-(x * x / numx + y * y / numy))
return win2d
def _create_matsign(nrow, ncol, pad):
listx = 1.0 * np.arange(0, ncol + 2 * pad)
listy = 1.0 * np.arange(0, nrow + 2 * pad)
x, y = np.meshgrid(listx, listy)
matsign = np.power(-1.0, x + y)
return matsign
def _2d_filter(mat, win2d, matsign, pad):
"""
Filtering an image using a 2D window.
Parameters
----------
mat : 2D array of floats
nrow : int
Height of the window.
ncol: int
Width of the window.
sigma: tuple of 2 floats
Sigmas of the window.
pad : int
Padding.
Returns
-------
Filtered image.
"""
matpad = np.pad(mat, ((0, 0), (pad, pad)), mode='edge')
matpad = np.pad(matpad, ((pad, pad), (0, 0)), mode='mean')
(nrow, ncol) = matpad.shape
matfilter = np.real(ifft2(fft2(matpad * matsign) * win2d) * matsign)
return matfilter[pad:nrow - pad, pad:ncol - pad]
def _rs_fit(sinogram, order, win2d, matsign, pad):
"""
Remove stripes using the fitting technique.
"""
(nrow, _) = sinogram.shape
nrow1 = nrow
if nrow1 % 2 == 0:
nrow1 = nrow1 - 1
if order >= nrow1:
order = nrow1 - 1
sinofit = savgol_filter(sinogram, nrow1, order, axis=0, mode='mirror')
sinofitsmooth = _2d_filter(sinofit, win2d, matsign, pad)
num1 = np.mean(sinofit)
num2 = np.mean(sinofitsmooth)
sinofitsmooth = num1 * sinofitsmooth / num2
return sinogram / sinofit * sinofitsmooth
def _remove_stripe_based_fitting(tomo, order, sigma):
pad = 50 # To reduce artifacts caused by FFT
nrow = tomo.shape[0]
ncol = tomo.shape[2]
win2d = _create_2d_window(nrow, ncol, sigma, pad)
matsign = _create_matsign(nrow, ncol, pad)
for m in range(tomo.shape[1]):
sino = tomo[:, m, :]
tomo[:, m, :] = _rs_fit(sino, order, win2d, matsign, pad)
[docs]def remove_large_stripe(tomo, snr=3, size=51, ncore=None, nchunk=None):
"""
Remove stripe artifacts from sinogram using Nghia Vo's
approach :cite:`Vo:18`
Algorithm 5 in the paper. Remove large stripes.
Parameters
----------
tomo : ndarray
3D tomographic data.
snr : float
Ratio used to locate of large stripes.
Greater is less sensitive.
size : int
Window size of the median filter.
ncore : int, optional
Number of cores that will be assigned to jobs.
nchunk : int, optional
Chunk size for each core.
Returns
-------
ndarray
Corrected 3D tomographic data.
"""
arr = mproc.distribute_jobs(
tomo,
func=_remove_large_stripe,
args=(snr, size),
axis=1,
ncore=ncore,
nchunk=nchunk)
return arr
def _detect_stripe(listdata, snr):
"""
Algorithm 4 in the paper. Used to locate stripes.
"""
numdata = len(listdata)
listsorted = np.sort(listdata)[::-1]
xlist = np.arange(0, numdata, 1.0)
ndrop = np.int16(0.25 * numdata)
(_slope, _intercept) = np.polyfit(
xlist[ndrop:-ndrop - 1], listsorted[ndrop:-ndrop - 1], 1)
numt1 = _intercept + _slope * xlist[-1]
noiselevel = np.abs(numt1 - _intercept)
val1 = np.abs(listsorted[0] - _intercept) / noiselevel
val2 = np.abs(listsorted[-1] - numt1) / noiselevel
listmask = np.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
def _rs_large(sinogram, snr, size, matindex):
"""
Remove large stripes by: locating stripes, normalizing to remove
full stripes, using the sorting technique to remove partial stripes.
"""
badpixelratio = 0.05
(nrow, ncol) = sinogram.shape
ndrop = np.int16(badpixelratio * nrow)
sinosorted = np.sort(sinogram, axis=0)
sinosmoothed = median_filter(sinosorted, (1, size))
list1 = np.mean(sinosorted[ndrop:nrow - ndrop], axis=0)
list2 = np.mean(sinosmoothed[ndrop:nrow - ndrop], axis=0)
listfact = list1 / list2
# Locate stripes
listmask = _detect_stripe(listfact, snr)
listmask = binary_dilation(listmask, iterations=1).astype(listmask.dtype)
matfact = np.tile(listfact, (nrow, 1))
# Normalize
sinogram = sinogram / matfact
sinogram1 = np.transpose(sinogram)
matcombine = np.asarray(np.dstack((matindex, sinogram1)))
matsort = np.asarray(
[row[row[:, 1].argsort()] for row in matcombine])
matsort[:, :, 1] = np.transpose(sinosmoothed)
matsortback = np.asarray(
[row[row[:, 0].argsort()] for row in matsort])
sino_corrected = np.transpose(matsortback[:, :, 1])
listxmiss = np.where(listmask > 0.0)[0]
sinogram[:, listxmiss] = sino_corrected[:, listxmiss]
return sinogram
def _remove_large_stripe(tomo, snr, size):
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)
[docs]def remove_dead_stripe(tomo, snr=3, size=51, ncore=None, nchunk=None):
"""
Remove stripe artifacts from sinogram using Nghia Vo's
approach :cite:`Vo:18`
Algorithm 6 in the paper. Remove unresponsive and fluctuating stripes.
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.
ncore : int, optional
Number of cores that will be assigned to jobs.
nchunk : int, optional
Chunk size for each core.
Returns
-------
ndarray
Corrected 3D tomographic data.
"""
arr = mproc.distribute_jobs(
tomo,
func=_remove_dead_stripe,
args=(snr, size),
axis=1,
ncore=ncore,
nchunk=nchunk)
return arr
def _rs_dead(sinogram, snr, size, matindex):
"""
Remove unresponsive and fluctuating stripes.
"""
sinogram = np.copy(sinogram) # Make it mutable
(nrow, _) = sinogram.shape
sinosmoothed = np.apply_along_axis(uniform_filter1d, 0, sinogram, 10)
listdiff = np.sum(np.abs(sinogram - sinosmoothed), axis=0)
nmean = np.mean(listdiff)
listdiffbck = median_filter(listdiff, size)
listdiffbck[listdiffbck == 0.0] = nmean
listfact = listdiff / listdiffbck
listmask = _detect_stripe(listfact, snr)
listmask = binary_dilation(listmask, iterations=1).astype(listmask.dtype)
listmask[0:2] = 0.0
listmask[-2:] = 0.0
listx = np.where(listmask < 1.0)[0]
listy = np.arange(nrow)
matz = sinogram[:, listx]
finter = interpolate.interp2d(listx, listy, matz, kind='linear')
listxmiss = np.where(listmask > 0.0)[0]
if len(listxmiss) > 0:
matzmiss = finter(listxmiss, listy)
sinogram[:, listxmiss] = matzmiss
# Use algorithm 5 to remove residual stripes
sinogram = _rs_large(sinogram, snr, size, matindex)
return sinogram
def _remove_dead_stripe(tomo, snr, size):
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)
[docs]def remove_all_stripe(tomo, snr=3, la_size=61, sm_size=21, ncore=None, nchunk=None):
"""
Remove stripe artifacts from sinogram using Nghia Vo's
approach :cite:`Vo:18`
Combine algorithms 6,5,4,3 to remove all types of stripes.
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.
ncore : int, optional
Number of cores that will be assigned to jobs.
nchunk : int, optional
Chunk size for each core.
Returns
-------
ndarray
Corrected 3D tomographic data.
"""
arr = mproc.distribute_jobs(
tomo,
func=_remove_all_stripe,
args=(snr, la_size, sm_size),
axis=1,
ncore=ncore,
nchunk=nchunk)
return arr
def _remove_all_stripe(tomo, snr, la_size, sm_size):
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)
tomo[:, m, :] = sino