diff options
Diffstat (limited to 'silx/image/tomography.py')
| -rw-r--r-- | silx/image/tomography.py | 169 |
1 files changed, 154 insertions, 15 deletions
diff --git a/silx/image/tomography.py b/silx/image/tomography.py index 66dc677..c2aedd8 100644 --- a/silx/image/tomography.py +++ b/silx/image/tomography.py @@ -32,27 +32,140 @@ __date__ = "12/09/2017" import numpy as np from math import pi +from itertools import product +from bisect import bisect from silx.math.fit import leastsq +# ------------------------------------------------------------------------------ +# -------------------- Filtering-related functions ----------------------------- +# ------------------------------------------------------------------------------ -def rescale_intensity(img, from_subimg=None, percentiles=None): +def compute_ramlak_filter(dwidth_padded, dtype=np.float32): """ - clamp intensity into the [2, 98] percentiles + Compute the Ramachandran-Lakshminarayanan (Ram-Lak) filter, used in + filtered backprojection. - :param img: - :param from_subimg: - :param percentiles: - :return: the rescale intensity + :param dwidth_padded: width of the 2D sinogram after padding + :param dtype: data type """ - 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 + 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 + + +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): @@ -158,3 +271,29 @@ def calc_center_centroid(sino): 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 + |