diff options
Diffstat (limited to 'src/silx/image/tomography.py')
-rw-r--r-- | src/silx/image/tomography.py | 301 |
1 files changed, 301 insertions, 0 deletions
diff --git a/src/silx/image/tomography.py b/src/silx/image/tomography.py new file mode 100644 index 0000000..53855c1 --- /dev/null +++ b/src/silx/image/tomography.py @@ -0,0 +1,301 @@ +# coding: utf-8 +# /*########################################################################## +# Copyright (C) 2017 European Synchrotron Radiation Facility +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in +# all copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +# THE SOFTWARE. +# +# ############################################################################*/ +""" +This module contains utilitary functions for tomography +""" + +__author__ = ["P. Paleo"] +__license__ = "MIT" +__date__ = "12/09/2017" + + +import numpy as np +from math import pi +from functools import lru_cache +from itertools import product +from bisect import bisect +from silx.math.fit import leastsq + +# ------------------------------------------------------------------------------ +# -------------------- Filtering-related functions ----------------------------- +# ------------------------------------------------------------------------------ + +def compute_ramlak_filter(dwidth_padded, dtype=np.float32): + """ + Compute the Ramachandran-Lakshminarayanan (Ram-Lak) filter, used in + filtered backprojection. + + :param dwidth_padded: width of the 2D sinogram after padding + :param dtype: data type + """ + L = dwidth_padded + h = np.zeros(L, dtype=dtype) + L2 = L//2+1 + h[0] = 1/4. + j = np.linspace(1, L2, L2//2, False).astype(dtype) # np < 1.9.0 + h[1:L2:2] = -1./(pi**2 * j**2) + h[L2:] = np.copy(h[1:L2-1][::-1]) + return h + + +def tukey(N, alpha=0.5): + """ + Compute the Tukey apodization window. + + :param int N: Number of points. + :param float alpha: + """ + apod = np.zeros(N) + x = np.arange(N)/(N-1) + r = alpha + M1 = (0 <= x) * (x < r/2) + M2 = (r/2 <= x) * (x <= 1 - r/2) + M3 = (1 - r/2 < x) * (x <= 1) + apod[M1] = (1 + np.cos(2*pi/r * (x[M1] - r/2)))/2. + apod[M2] = 1. + apod[M3] = (1 + np.cos(2*pi/r * (x[M3] - 1 + r/2)))/2. + return apod + + +def lanczos(N): + """ + Compute the Lanczos window (truncated sinc) of width N. + + :param int N: window width + """ + x = np.arange(N)/(N-1) + return np.sin(pi*(2*x-1))/(pi*(2*x-1)) + + +def compute_fourier_filter(dwidth_padded, filter_name, cutoff=1.): + """ + Compute the filter used for FBP. + + :param dwidth_padded: padded detector width. As the filtering is done by the + Fourier convolution theorem, dwidth_padded should be at least 2*dwidth. + :param filter_name: Name of the filter. Available filters are: + Ram-Lak, Shepp-Logan, Cosine, Hamming, Hann, Tukey, Lanczos. + :param cutoff: Cut-off frequency, if relevant. + """ + Nf = dwidth_padded + #~ filt_f = np.abs(np.fft.fftfreq(Nf)) + rl = compute_ramlak_filter(Nf, dtype=np.float64) + filt_f = np.fft.fft(rl) + + filter_name = filter_name.lower() + if filter_name in ["ram-lak", "ramlak"]: + return filt_f + + w = 2 * pi * np.fft.fftfreq(dwidth_padded) + d = cutoff + apodization = { + # ~OK + "shepp-logan": np.sin(w[1:Nf]/(2*d))/(w[1:Nf]/(2*d)), + # ~OK + "cosine": np.cos(w[1:Nf]/(2*d)), + # OK + "hamming": 0.54*np.ones_like(filt_f)[1:Nf] + .46 * np.cos(w[1:Nf]/d), + # OK + "hann": (np.ones_like(filt_f)[1:Nf] + np.cos(w[1:Nf]/d))/2., + # These one is not compatible with Astra - TODO investigate why + "tukey": np.fft.fftshift(tukey(dwidth_padded, alpha=d/2.))[1:Nf], + "lanczos": np.fft.fftshift(lanczos(dwidth_padded))[1:Nf], + } + if filter_name not in apodization: + raise ValueError("Unknown filter %s. Available filters are %s" % + (filter_name, str(apodization.keys()))) + filt_f[1:Nf] *= apodization[filter_name] + return filt_f + + +@lru_cache(maxsize=1) +def generate_powers(): + """ + Generate a list of powers of [2, 3, 5, 7], + up to (2**15)*(3**9)*(5**6)*(7**5). + """ + primes = [2, 3, 5, 7] + maxpow = {2: 15, 3: 9, 5: 6, 7: 5} + valuations = [] + for prime in primes: + # disallow any odd number (for R2C transform), and any number + # not multiple of 4 (Ram-Lak filter behaves strangely when + # dwidth_padded/2 is not even) + minval = 2 if prime == 2 else 0 + valuations.append(range(minval, maxpow[prime]+1)) + powers = product(*valuations) + res = [] + for pw in powers: + res.append(np.prod(list(map(lambda x : x[0]**x[1], zip(primes, pw))))) + return np.unique(res) + + +def get_next_power(n, powers=None): + """ + Given a number, get the closest (upper) number p such that + p is a power of 2, 3, 5 and 7. + """ + if powers is None: + powers = generate_powers() + idx = bisect(powers, n) + if powers[idx-1] == n: + return n + return powers[idx] + + +# ------------------------------------------------------------------------------ +# ------------- Functions for determining the center of rotation -------------- +# ------------------------------------------------------------------------------ + + + +def calc_center_corr(sino, fullrot=False, props=1): + """ + Compute a guess of the Center of Rotation (CoR) of a given sinogram. + The computation is based on the correlation between the line projections at + angle (theta = 0) and at angle (theta = 180). + + Note that for most scans, the (theta=180) angle is not included, + so the CoR might be underestimated. + In a [0, 360[ scan, the projection angle at (theta=180) is exactly in the + middle for odd number of projections. + + :param numpy.ndarray sino: Sinogram + :param bool fullrot: optional. If False (default), the scan is assumed to + be [0, 180). + If True, the scan is assumed to be [0, 380). + :param int props: optional. Number of propositions for the CoR + """ + + n_a, n_d = sino.shape + first = 0 + last = -1 if not(fullrot) else n_a // 2 + proj1 = sino[first, :] + proj2 = sino[last, :][::-1] + + # Compute the correlation in the Fourier domain + proj1_f = np.fft.fft(proj1, 2 * n_d) + proj2_f = np.fft.fft(proj2, 2 * n_d) + corr = np.abs(np.fft.ifft(proj1_f * proj2_f.conj())) + + if props == 1: + pos = np.argmax(corr) + if pos > n_d // 2: + pos -= n_d + return (n_d + pos) / 2. + else: + corr_argsorted = np.argsort(corr)[:props] + corr_argsorted[corr_argsorted > n_d // 2] -= n_d + return (n_d + corr_argsorted) / 2. + + +def _sine_function(t, offset, amplitude, phase): + """ + Helper function for calc_center_centroid + """ + n_angles = t.shape[0] + res = amplitude * np.sin(2 * pi * (1. / (2 * n_angles)) * t + phase) + return offset + res + + +def _sine_function_derivative(t, params, eval_idx): + """ + Helper function for calc_center_centroid + """ + offset, amplitude, phase = params + n_angles = t.shape[0] + w = 2.0 * pi * (1. / (2.0 * n_angles)) * t + phase + grad = (1.0, np.sin(w), amplitude*np.cos(w)) + return grad[eval_idx] + + +def calc_center_centroid(sino): + """ + Compute a guess of the Center of Rotation (CoR) of a given sinogram. + The computation is based on the computation of the centroid of each + projection line, which should be a sine function according to the + Helgason-Ludwig condition. + This method is unlikely to work in local tomography. + + :param numpy.ndarray sino: Sinogram + """ + + n_a, n_d = sino.shape + # Compute the vector of centroids of the sinogram + i = np.arange(n_d) + centroids = np.sum(sino*i, axis=1)/np.sum(sino, axis=1) + + # Fit with a sine function : phase, amplitude, offset + # Using non-linear Levenberg–Marquardt algorithm + angles = np.linspace(0, n_a, n_a, True) + # Initial parameter vector + cmax, cmin = centroids.max(), centroids.min() + offs = (cmax + cmin) / 2. + amp = (cmax - cmin) / 2. + phi = 1.1 + p0 = (offs, amp, phi) + + constraints = np.zeros((3, 3)) + + popt, _ = leastsq(model=_sine_function, + xdata=angles, + ydata=centroids, + p0=p0, + sigma=None, + constraints=constraints, + model_deriv=None, + epsfcn=None, + deltachi=None, + full_output=0, + check_finite=True, + left_derivative=False, + max_iter=100) + return popt[0] + + + +# ------------------------------------------------------------------------------ +# -------------------- Visualization-related functions ------------------------- +# ------------------------------------------------------------------------------ + + +def rescale_intensity(img, from_subimg=None, percentiles=None): + """ + clamp intensity into the [2, 98] percentiles + + :param img: + :param from_subimg: + :param percentiles: + :return: the rescale intensity + """ + if percentiles is None: + percentiles = [2, 98] + else: + assert type(percentiles) in (tuple, list) + assert(len(percentiles) == 2) + data = from_subimg if from_subimg is not None else img + imin, imax = np.percentile(data, percentiles) + res = np.clip(img, imin, imax) + return res + |