# ============================================================================
# ============================================================================
# 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 correction techniques.
# Contributors:
# ============================================================================
"""
Module of correction methods in the preprocessing stage:
- Flat-field correction.
- Distortion correction.
- MTF deconvolution.
- Tilted sinogram generation.
- Tilted 1D intensity-profile generation.
- Beam hardening correction.
- Sinogram upsampling.
"""
import numpy as np
import numpy.fft as fft
from scipy.ndimage import map_coordinates
import algotom.prep.removal as remo
import algotom.prep.filtering as filt
import algotom.prep.phase as ps
import algotom.util.utility as util
import algotom.prep.conversion as conv
[docs]def flat_field_correction(proj, flat, dark, ratio=1.0, use_dark=True,
**options):
"""
Perform flat-field correction with options to remove zinger artifacts
and/or stripe artifacts.
Parameters
----------
proj : array_like
3D or 2D array. Projection images or a sinogram image.
flat : array_like
2D or 1D array. Flat-field image or a single row of it.
dark : array_like
2D or 1D array. Dark-field image or a single row of it.
ratio : float
Ratio between exposure time used for recording projections
and exposure time used for recording flat field.
use_dark : bool
Subtracting dark field if True.
options : dict, optional
Apply a zinger removal method and/or ring removal methods.
E.g option1={"method": "dezinger", "para1": 0.001, "para2": 1},
option2={"method": "remove_stripe_based_sorting",\
"para1": 15, "para2": 1}
Returns
-------
array_like
3D or 2D array. Corrected projections or corrected sinograms.
"""
msg = "\n Please use the dictionary format: options={'method':" \
" 'filter_name', 'para1': parameter_1, 'para2': parameter_2}"
flat = ratio * flat
if use_dark:
flat_dark = flat - dark
if 0.0 in flat_dark:
nmean = np.mean(flat_dark)
if nmean != 0.0:
flat_dark[flat_dark == 0.0] = nmean
else:
flat_dark[flat_dark == 0.0] = 1
proj_corr = (np.float32(proj) - dark) / flat_dark
else:
proj_corr = (np.float32(proj) - dark) / flat_dark
else:
if 0.0 in flat:
nmean = np.mean(flat)
if nmean != 0.0:
flat[flat == 0.0] = nmean
else:
flat[flat == 0.0] = 1
proj_corr = np.float32(proj) / flat
else:
proj_corr = np.float32(proj) / flat
if len(options) != 0:
for opt_name in options:
opt = options[opt_name]
if isinstance(opt, dict):
method = tuple(opt.values())[0]
para = tuple(opt.values())[1:]
if proj_corr.ndim == 2:
if method in dir(remo):
proj_corr = getattr(remo, method)(proj_corr, *para)
elif method in dir(filt):
proj_corr = getattr(filt, method)(proj_corr, *para)
elif method in dir(ps):
proj_corr = getattr(ps, method)(proj_corr, *para)
else:
raise ValueError("Can't find the method: '{}' in"
" the namespace".format(method))
else:
for i in np.arange(proj_corr.shape[1]):
if method in dir(remo):
proj_corr[:, i, :] = getattr(remo, method)(
proj_corr[:, i, :], *para)
elif method in dir(filt):
proj_corr[:, i, :] = getattr(filt, method)(
proj_corr[:, i, :], *para)
elif method in dir(ps):
proj_corr[:, i, :] = getattr(ps, method)(
proj_corr[:, i, :], *para)
else:
raise ValueError("Can't find the method: '{}' in "
"the namespace".format(method))
else:
if opt is not None:
raise ValueError(msg)
return proj_corr
[docs]def unwarp_projection(proj, xcenter, ycenter, list_fact):
"""
Apply distortion correction to a projection image using the polynomial
backward model (Ref. [1]).
Parameters
----------
proj : array_like
2D array. Projection image.
xcenter : float
Center of distortion in x-direction.
ycenter : float
Center of distortion in y-direction.
list_fact : list of float
Polynomial coefficients of the backward model.
Returns
-------
array_like
2D array. Distortion corrected.
References
----------
[1] : https://doi.org/10.1364/OE.23.032859
"""
(height, width) = proj.shape
xu_list = np.arange(width) - xcenter
yu_list = np.arange(height) - ycenter
xu_mat, yu_mat = np.meshgrid(xu_list, yu_list)
ru_mat = np.sqrt(xu_mat ** 2 + yu_mat ** 2)
fact_mat = np.sum(np.asarray(
[factor * ru_mat ** i for i, factor in enumerate(list_fact)]), axis=0)
xd_mat = np.float32(np.clip(xcenter + fact_mat * xu_mat, 0, width - 1))
yd_mat = np.float32(np.clip(ycenter + fact_mat * yu_mat, 0, height - 1))
indices = np.reshape(yd_mat, (-1, 1)), np.reshape(xd_mat, (-1, 1))
proj = map_coordinates(proj, indices, order=1, mode='reflect')
return proj.reshape((height, width))
[docs]def unwarp_sinogram(data, index, xcenter, ycenter, list_fact, **option):
"""
Unwarp sinogram [:,index.:] of a 3D tomographic dataset or
a hdf/nxs object.
Parameters
----------
data : array_like or hdf object
3D array.
index : int
Index of the sinogram.
xcenter : float
Center of distortion in x-direction.
ycenter : float
Center of distortion in y-direction.
list_fact : list of float
Polynomial coefficients of the backward model.
option : list or tuple of int
To extract subset data along axis 0 from a hdf object. E.g option =
(start, stop, step)
Returns
-------
array_like
2D array. Distortion-corrected sinogram.
"""
if data.ndim != 3:
raise ValueError("Input must be a 3D data !!!")
(depth, height, width) = data.shape
if len(option) != 0:
opt = option[list(option.keys())[0]]
start, stop, step = opt
else:
start, stop, step = 0, depth, 1
list_idx = range(start, stop, step)
depth = len(list_idx)
xu_list = np.arange(0, width) - xcenter
yu = index - ycenter
ru_list = np.sqrt(xu_list ** 2 + yu ** 2)
flist = np.sum(np.asarray(
[factor * ru_list ** i for i, factor in enumerate(list_fact)]), axis=0)
xd_list = np.clip(xcenter + flist * xu_list, 0, width - 1)
yd_list = np.clip(ycenter + flist * yu, 0, height - 1)
yd_min = np.int16(np.floor(np.amin(yd_list)))
yd_max = np.int16(np.ceil(np.amax(yd_list))) + 1
yd_list = yd_list - yd_min
indices = yd_list, xd_list
sinogram = np.zeros((depth, width), dtype=np.float32)
for i, idx in enumerate(list_idx):
sinogram[i] = map_coordinates(
data[idx, yd_min:yd_max, :], indices, order=1, mode='reflect')
return sinogram
[docs]def unwarp_sinogram_chunk(data, start_index, stop_index, xcenter, ycenter,
list_fact, **option):
"""
Unwarp chunk of sinograms [:, start_index: stop_index, :]
of a 3D tomographic dataset or a hdf/nxs object.
Parameters
----------
data : array_like or hdf object
3D array.
start_index : int
Starting index of sinograms.
stop_index : int
Stopping index of sinograms.
xcenter : float
Center of distortion in x-direction.
ycenter : float
Center of distortion in y-direction.
list_fact : list of float
Polynomial coefficients of the backward model.
option : list or tuple of int
To extract subset data along axis 0 from a hdf object. E.g option =
[start, stop, step]
Returns
-------
array_like
3D array. Distortion corrected.
"""
if data.ndim != 3:
raise ValueError("Input must be a 3D data !!!")
(depth, height, width) = data.shape
if len(option) != 0:
opt = option[list(option.keys())[0]]
start, stop, step = opt
else:
start, stop, step = 0, depth, 1
list_idx = range(start, stop, step)
if stop_index == -1:
stop_index = height
xu_list = np.arange(0, width) - xcenter
yu1 = start_index - ycenter
ru_list = np.sqrt(xu_list ** 2 + yu1 ** 2)
flist = np.sum(np.asarray(
[factor * ru_list ** i for i, factor in enumerate(list_fact)]), axis=0)
yd_list1 = np.clip(ycenter + flist * yu1, 0, height - 1)
yu2 = stop_index - ycenter
ru_list = np.sqrt(xu_list ** 2 + yu2 ** 2)
flist = np.sum(np.asarray(
[factor * ru_list ** i for i, factor in enumerate(list_fact)]), axis=0)
yd_list2 = np.clip(ycenter + flist * yu2, 0, height - 1)
yd_min = np.int16(np.floor(np.amin(yd_list1)))
yd_max = np.int16(np.ceil(np.amax(yd_list2)))
yu_list = np.arange(start_index, stop_index) - ycenter
xu_mat, yu_mat = np.meshgrid(xu_list, yu_list)
ru_mat = np.sqrt(xu_mat ** 2 + yu_mat ** 2)
fact_mat = np.sum(np.asarray(
[factor * ru_mat ** i for i, factor in enumerate(list_fact)]), axis=0)
xd_mat = np.float32(np.clip(xcenter + fact_mat * xu_mat, 0, width - 1))
yd_mat = np.float32(
np.clip(ycenter + fact_mat * yu_mat, 0, height - 1)) - yd_min
sino_chunk = np.asarray(
[util.mapping(data[i, yd_min:yd_max, :], xd_mat, yd_mat) for i in
list_idx])
return sino_chunk
[docs]def mtf_deconvolution(mat, window, pad):
"""
Deconvolve a projection-image using division in the Fourier domain.
Window can be determined using the approach in Ref. [1].
Parameters
----------
mat : array_like
2D array. Projection image.
window : array_like
2D array. MTF function.
pad : int
Padding width to reduce the side effects of the Fourier transform.
Returns
-------
array_like
2D array. Deconvolved image.
References
----------
[1] : https://doi.org/10.1117/12.2530324
"""
(height1, width1) = mat.shape
(height2, width2) = window.shape
if (height1 > height2) or (width1 > width2):
raise ValueError(
"The sizes of the image are larger than the sizes of the window")
else:
pad_row = height2 - height1
pad_col = width2 - width1
mat = np.pad(mat, ((0, pad_row), (0, pad_col)), mode="edge")
mat_pad = np.pad(mat, pad, mode="reflect")
win_pad = np.pad(window, pad, mode="constant", constant_values=1.0)
mat_dec = np.real(fft.ifft2(fft.fft2(mat_pad) / fft.ifftshift(win_pad)))
mat_dec = mat_dec[pad: height2 + pad, pad: width2 + pad]
return mat_dec[:height1, :width1]
[docs]def generate_tilted_sinogram(data, index, angle, **option):
"""
Generate a tilted sinogram of a 3D tomographic dataset or a hdf/nxs object.
Parameters
----------
data : array_like or hdf object
3D array.
index : int
Index of the sinogram.
angle : float
Tilted angle in degree.
option : list or tuple of int
To extract subset data along axis 0 from a hdf object. E.g option =
(start, stop, step)
Returns
-------
array_like
2D array. Tilted sinogram.
"""
if data.ndim != 3:
raise ValueError("Input must be a 3D data !!!")
(depth, height, width) = data.shape
if len(option) != 0:
opt = option[list(option.keys())[0]]
start, stop, step = opt
else:
start, stop, step = 0, depth, 1
list_idx = range(start, stop, step)
depth = len(list_idx)
x_cen = (width - 1.0) / 2
y_cen = (height - 1.0) / 2
xlist = np.arange(0, width) - x_cen
y = index - y_cen
angle = angle * np.pi / 180.0
xt_list = xlist * np.cos(angle) - y * np.sin(angle)
yt_list = xlist * np.sin(angle) + y * np.cos(angle)
xt_list = np.clip(x_cen + xt_list, 0, width - 1)
yt_list = np.clip(y_cen + yt_list, 0, height - 1)
yt_min = np.int16(np.floor(np.amin(yt_list)))
yt_max = np.int16(np.ceil(np.amax(yt_list))) + 1
yt_list = yt_list - yt_min
indices = yt_list, xt_list
sinogram = np.zeros((depth, width), dtype=np.float32)
for i, idx in enumerate(list_idx):
sinogram[i] = map_coordinates(
data[idx, yt_min:yt_max, :], indices, order=1, mode='reflect')
return sinogram
[docs]def generate_tilted_sinogram_chunk(data, start_index, stop_index, angle,
**option):
"""
Generate a chunk of tilted sinograms of a 3D tomographic dataset or a
hdf/nxs object.
Parameters
----------
data : array_like or hdf object
3D array.
start_index : int
Starting index of sinograms.
stop_index : int
Stopping index of sinograms.
angle : float
Tilted angle in degree.
option : list or tuple of int
To extract subset data along axis 0 from a hdf object. E.g option =
(start, stop, step)
Returns
-------
array_like
3D array. Chunk of tilted sinograms.
"""
if data.ndim != 3:
raise ValueError("Input must be a 3D data !!!")
(depth, height, width) = data.shape
if len(option) != 0:
opt = option[list(option.keys())[0]]
start, stop, step = opt
else:
start, stop, step = 0, depth, 1
list_idx = range(start, stop, step)
x_cen = (width - 1.0) / 2
y_cen = (height - 1.0) / 2
angle = angle * np.pi / 180.0
xlist = np.arange(0, width) - x_cen
ylist = np.arange(start_index, stop_index + 1) - y_cen
x_mat, y_mat = np.meshgrid(xlist, ylist)
x_mat1 = x_mat * np.cos(angle) - y_mat * np.sin(angle)
y_mat1 = x_mat * np.sin(angle) + y_mat * np.cos(angle)
x_mat1 = np.clip(x_cen + x_mat1, 0, width - 1)
y_mat1 = np.clip(y_cen + y_mat1, 0, height - 1)
y_min = np.int16(np.floor(np.amin(y_mat1)))
y_max = np.int16(np.ceil(np.amax(y_mat1))) + 1
y_mat1 = y_mat1 - y_min
sino_chunk = np.asarray(
[util.mapping(data[i, y_min:y_max, :], x_mat1, y_mat1) for i in
list_idx])
return sino_chunk
[docs]def generate_tilted_profile_line(mat, index, angle):
"""
Generate a tilted horizontal intensity-profile of an image.
Parameters
----------
mat : array_like
2D array.
index : int
Index of the line.
angle : float
Tilted angle in degree.
Returns
-------
array_like
1D array.
"""
if mat.ndim != 2:
raise ValueError("Input must be a 2D array !!!")
(height, width) = mat.shape
x_cen = (width - 1.0) / 2
y_cen = (height - 1.0) / 2
xlist = np.arange(0, width) - x_cen
y = index - y_cen
angle = angle * np.pi / 180.0
xt_list = xlist * np.cos(angle) - y * np.sin(angle)
yt_list = xlist * np.sin(angle) + y * np.cos(angle)
xt_list = np.clip(x_cen + xt_list, 0, width - 1)
yt_list = np.clip(y_cen + yt_list, 0, height - 1)
yt_min = np.int16(np.floor(np.amin(yt_list)))
yt_max = np.int16(np.ceil(np.amax(yt_list))) + 1
yt_list = yt_list - yt_min
indices = yt_list, xt_list
profile_line = map_coordinates(mat[yt_min:yt_max, :], indices, order=1,
mode='reflect')
return profile_line
[docs]def generate_tilted_profile_chunk(mat, start_index, stop_index, angle):
"""
Generate a chunk of tilted horizontal intensity-profiles of an image.
Parameters
----------
mat : array_like
2D array.
start_index : int
Starting index of lines.
stop_index : int
Stopping index of lines.
angle : float
Tilted angle in degree.
Returns
-------
array_like
2D array.
"""
if mat.ndim != 2:
raise ValueError("Input must be a 2D data !!!")
(height, width) = mat.shape
x_cen = (width - 1.0) / 2
y_cen = (height - 1.0) / 2
angle = angle * np.pi / 180.0
xlist = np.arange(0, width) - x_cen
ylist = np.arange(start_index, stop_index + 1) - y_cen
x_mat, y_mat = np.meshgrid(xlist, ylist)
x_mat1 = x_mat * np.cos(angle) - y_mat * np.sin(angle)
y_mat1 = x_mat * np.sin(angle) + y_mat * np.cos(angle)
x_mat1 = np.clip(x_cen + x_mat1, 0, width - 1)
y_mat1 = np.clip(y_cen + y_mat1, 0, height - 1)
y_min = np.int16(np.floor(np.amin(y_mat1)))
y_max = np.int16(np.ceil(np.amax(y_mat1))) + 1
y_mat1 = y_mat1 - y_min
profile_chunk = util.mapping(mat[y_min:y_max, :], x_mat1, y_mat1)
return profile_chunk
[docs]def non_linear_function(intensity, q, n, opt=True):
"""
Function used to define the response curve.
Parameters
----------
intensity : float
Values stay in the range of [0; 1]
q : float
Positive number.
n : float
Positive number. Must larger than 1.
opt : bool
True: Curve more to values closer to 1.0.
False: Curve more to values closer to 0.0
Returns
-------
float
"""
if opt:
x = 1.0 - intensity
num = np.log(1.0 - q * (1 - n))
result = 1.0 - (np.log(1.0 - q * (x**n - n * x)) / num)
else:
x = intensity
num = np.log(1.0 - q * (1 - n))
result = (np.log(1.0 - q * (x**n - n * x)) / num)
return result
[docs]def beam_hardening_correction(mat, q, n, opt=True):
"""
Correct the grayscale values of a normalized image using a non-linear
function.
Parameters
----------
mat : array_like
Normalized projection image or sinogram image.
q : float
Positive number. Recommended range [0.005, 50].
n : float
Positive number. Must larger than 1.
opt : bool
True: Curve towards 0.0.
False: Curve towards 1.0.
Returns
-------
array_like
Corrected image.
"""
if np.max(mat) >= 2.0:
raise ValueError("Input image must be normalized, i.e. gray-scales"
" are in the range of [0.0, 1.0]) !!!")
if n < 2.0:
raise ValueError("!!! n must be larger than or equal to 2 !!!")
return np.asarray([non_linear_function(x, q, n, opt) for x in mat])
[docs]def upsample_sinogram(sinogram, scale, center=0, sino_type="180", iteration=1,
pad=50):
"""
Upsample a sinogram-image along angular direction based on the
double-wedge filter (Ref. [1]).
Parameters
----------
sinogram : array_like
2D array. Sinogram image.
scale : int
Upscale 2n_x time. E.g. 2, 4, 6.
center : float, optional
Center-of-rotation. No need for a 360-sinogram.
sino_type : {"180", "360"}
Sinogram type : 180-degree or 360-degree.
iteration : int, optional
Number of iteration for the double-wedge filter.
pad : int, optional
Padding width for FFT.
Returns
-------
array_like
Upsampled sinogram.
"""
scale = np.clip(int(scale), 1, None)
if scale % 2 != 0:
raise ValueError("Scaling value must be an even number and starting "
"from 2!!!")
if sino_type == "180":
if center > 0:
sino_360 = conv.convert_sinogram_180_to_360(sinogram, center)
else:
raise ValueError("Please input the center-of-rotation !!!")
else:
sino_360 = np.copy(sinogram)
num_iter = int(scale / 2)
for i in range(num_iter):
(height, width) = sino_360.shape
db_height = 2 * height - 1
sino_tmp = np.zeros((db_height, width), dtype=np.float32)
sino_tmp[0:2 * height:2] = sino_360
sino_tmp = filt.double_wedge_filter(sino_tmp, sino_type="360", ratio=1,
iteration=iteration, pad=pad)
num = np.mean(sino_tmp)
sino_tmp = sino_tmp * np.mean(sino_360) / num
sino_tmp[0:2 * height:2] = sino_360
sino_tmp = filt.double_wedge_filter(sino_tmp, sino_type="360", ratio=1,
iteration=iteration, pad=pad)
sino_360 = sino_tmp
if sino_type == "180":
height = sino_360.shape[0]
return sino_360[:height // 2 + 1]
else:
return sino_360