Skip to content
GitLab
Commit 9eb7a31c authored by Nicola Vigano's avatar Nicola Vigano
Browse files

Alignment: added band-pass filter computation and perform filtering in correlation 


Signed-off-by: Nicola Vigano's avatarNicola VIGANÒ <nicola.vigano@esrf.fr>
parent 3353b97e
Hide whitespace changes
Inline Side-by-side
nabu/misc/fourier_filters.py
View file @ 9eb7a31c
......@@ -7,14 +7,16 @@ import numpy as np
try:
import scipy.special as spspe
__have_scipy__ = True
except ImportError:
import math
__have_scipy__ = False
def get_lowpass_filter(img_shape, cutoff_par):
"""Computes a low pass filter with the erfc.
def get_lowpass_filter(img_shape, cutoff_par=None):
"""Computes a low pass filter using the erfc function.
Parameters
----------
......@@ -38,7 +40,9 @@ def get_lowpass_filter(img_shape, cutoff_par):
numpy.array_like
The computed filter
"""
if isinstance(cutoff_par, (int, float)):
if cutoff_par is None:
return 1
elif isinstance(cutoff_par, (int, float)):
cutoff_pix = cutoff_par
cutoff_trans_fact = None
else:
......@@ -69,13 +73,13 @@ def get_lowpass_filter(img_shape, cutoff_par):
else:
res = np.array(list(map(math.erfc, k_pos_rescaled))) / 2
else:
res = np.exp(- (np.pi ** 2) * (r ** 2) * (cutoff_pix ** 2) * 2)
res = np.exp(-(np.pi ** 2) * (r ** 2) * (cutoff_pix ** 2) * 2)
return res
def get_highpass_filter(img_shape, cutoff_par):
"""Computes a high pass filter with the erfc.
def get_highpass_filter(img_shape, cutoff_par=None):
"""Computes a high pass filter using the erfc function.
Parameters
----------
......@@ -99,4 +103,41 @@ def get_highpass_filter(img_shape, cutoff_par):
numpy.array_like
The computed filter
"""
return 1 - get_lowpass_filter(img_shape, cutoff_par)
if cutoff_par is None:
return 1
else:
return 1 - get_lowpass_filter(img_shape, cutoff_par)
def get_bandpass_filter(img_shape, cutoff_lowpass=None, cutoff_highpass=None):
"""Computes a band pass filter using the erfc function.
The cutoff structures should be formed as follows:
- tuple of two floats: the first indicates the cutoff frequency, the second \
determines the width of the transition region, as fraction of the cutoff frequency.
- one float -> it represents the sigma of a gaussian which acts as a filter or anti-filter (1 - filter).
Parameters
----------
img_shape: tuple
Shape of the image
cutoff_lowpass: float or sequence of two floats
Cutoff parameters for the low-pass filter
cutoff_highpass: float or sequence of two floats
Cutoff parameters for the high-pass filter
Raises
------
ValueError
In case of malformed cutoff_par
Returns
-------
numpy.array_like
The computed filter
"""
return get_lowpass_filter(img_shape, cutoff_par=cutoff_lowpass) * get_highpass_filter(
img_shape, cutoff_par=cutoff_highpass
)
nabu/preproc/alignment.py
View file @ 9eb7a31c
......@@ -202,7 +202,7 @@ class AlignmentBase(object):
return roi_yxhw
@staticmethod
def _prepare_image(img, invalid_val=1e-5, roi_yxhw=None, median_filt_shape=None, high_pass=None, low_pass=None):
def _prepare_image(img, invalid_val=1e-5, roi_yxhw=None, median_filt_shape=None, low_pass=None, high_pass=None):
"""
Prepare and returns a cropped and filtered image, or array of filtered images if the input is an array of images.
......@@ -214,20 +214,10 @@ class AlignmentBase(object):
value to be used in replacement of nan and inf values
median_filt_shape: int or sequence of int
the width or the widths of the median window
low_pass: float or sequence of two floats
Low-pass filter properties, as described in `nabu.misc.fourier_filters`
high_pass: float or sequence of two floats
Position of the cut off in pixels, if a sequence is given the second float expresses the
width of the transition region which is given as a fraction of the cutoff frequency.
When only one float is given for this argument a gaussian is applied whose sigma is the
parameter, and the result is subtracted from 1 to obtain the high pass filter
When a sequence of two numbers is given then the filter is 1 ( no filtering) above the cutoff
frequency while a smooth transition to zero is done for smaller frequency
low_pass: float or sequence of two floats
Position of the cut off in pixels, if a sequence is given the second float expresses the
width of the transition region which is given as a fraction of the cutoff frequency.
When only one float is given for this argument a gaussian is applied whose sigma is the
parameter.
When a sequence of two numbers is given then the filter is 1 ( no filtering) till the cutoff
frequency and then a smooth erfc transition to zero is done
High-pass filter properties, as described in `nabu.misc.fourier_filters`
Returns
-------
......@@ -248,11 +238,9 @@ class AlignmentBase(object):
img[np.isinf(img)] = invalid_val
if high_pass is not None or low_pass is not None:
img_filter = np.ones(img.shape[-2:], dtype=img.dtype)
if low_pass is not None:
img_filter[:] *= fourier_filters.get_lowpass_filter(img.shape[-2:], low_pass)
if high_pass is not None:
img_filter[:] *= fourier_filters.get_highpass_filter(img.shape[-2:], high_pass)
img_filter = fourier_filters.get_bandpass_filter(
img.shape[-2:], cutoff_lowpass=low_pass, cutoff_highpass=high_pass
)
# fft2 and iff2 use axes=(-2, -1) by default
img = np.fft.ifft2(np.fft.fft2(img) * img_filter).real
......@@ -278,7 +266,7 @@ class AlignmentBase(object):
return img
@staticmethod
def _compute_correlation_fft(img_1, img_2, padding_mode, axes=(-2, -1)):
def _compute_correlation_fft(img_1, img_2, padding_mode, axes=(-2, -1), low_pass=None, high_pass=None):
do_circular_conv = padding_mode is None or padding_mode == "wrap"
if not do_circular_conv:
img_shape = img_2.shape
......@@ -293,15 +281,22 @@ class AlignmentBase(object):
# compute fft's of the 2 images
img_fft_1 = np.fft.fftn(img_1, axes=axes)
img_fft_2 = np.conjugate(np.fft.fftn(img_2, axes=axes))
img_prod = img_fft_1 * img_fft_2
if low_pass is not None or high_pass is not None:
filt = fourier_filters.get_bandpass_filter(img_prod.shape[-2:], cutoff_lowpass=low_pass, cutoff_highpass=high_pass)
img_prod *= filt
# inverse fft of the product to get cross_correlation of the 2 images
cc = np.real(np.fft.ifftn(img_fft_1 * img_fft_2, axes=axes))
cc = np.real(np.fft.ifftn(img_prod, axes=axes))
if not do_circular_conv:
cc = np.fft.fftshift(cc, axes=axes)
slicing = [slice(None)] * len(img_shape)
for a in axes:
slicing[a] = slice(pad_size[a], cc.shape[a]-pad_size[a])
slicing[a] = slice(pad_size[a], cc.shape[a] - pad_size[a])
cc = cc[tuple(slicing)]
cc = np.fft.ifftshift(cc, axes=axes)
......@@ -396,21 +391,10 @@ class CenterOfRotation(AlignmentBase):
peak_fit_radius: int, optional
Radius size around the max correlation pixel, for sub-pixel fitting.
Minimum and default value is 1.
high_pass: float or sequence of two floats
Position of the cut off in pixels, if a sequence is given the second float expresses the
width of the transition region which is given as a fraction of the cutoff frequency.
When only one float is given for this argument a gaussian is applied whose sigma is the
parameter, and the result is subtracted from 1 to obtain the high pass filter
When a sequence of two numbers is given then the filter is 1 ( no filtering) above the cutoff
frequency while a smooth transition to zero is done for smaller frequency
low_pass: float or sequence of two floats
Position of the cut off in pixels, if a sequence is given the second float expresses the
width of the transition region which is given as a fraction of the cutoff frequency.
When only one float is given for this argument a gaussian is applied whose sigma is the
parameter.
When a sequence of two numbers is given then the filter is 1 ( no filtering) till the cutoff
frequency and then a smooth erfc transition to zero is done
Low-pass filter properties, as described in `nabu.misc.fourier_filters`
high_pass: float or sequence of two floats
High-pass filter properties, as described in `nabu.misc.fourier_filters`
Raises
------
......@@ -448,14 +432,10 @@ class CenterOfRotation(AlignmentBase):
img_shape = img_2.shape
roi_yxhw = self._determine_roi(img_shape, roi_yxhw, self.truncate_horz_pow2)
img_1 = self._prepare_image(
img_1, roi_yxhw=roi_yxhw, median_filt_shape=median_filt_shape, high_pass=high_pass, low_pass=low_pass
)
img_2 = self._prepare_image(
img_2, roi_yxhw=roi_yxhw, median_filt_shape=median_filt_shape, high_pass=high_pass, low_pass=low_pass
)
img_1 = self._prepare_image(img_1, roi_yxhw=roi_yxhw, median_filt_shape=median_filt_shape)
img_2 = self._prepare_image(img_2, roi_yxhw=roi_yxhw, median_filt_shape=median_filt_shape)
cc = self._compute_correlation_fft(img_1, img_2, padding_mode)
cc = self._compute_correlation_fft(img_1, img_2, padding_mode, high_pass=high_pass, low_pass=low_pass)
img_shape = img_2.shape
cc_vs = np.fft.fftfreq(img_shape[-2], 1 / img_shape[-2])
cc_hs = np.fft.fftfreq(img_shape[-1], 1 / img_shape[-1])
......@@ -518,6 +498,8 @@ class DetectorTranslationAlongBeam(AlignmentBase):
median_filt_shape=None,
padding_mode=None,
peak_fit_radius=1,
high_pass=None,
low_pass=None,
equispaced_increments=True,
return_shifts=False,
):
......@@ -553,6 +535,10 @@ class DetectorTranslationAlongBeam(AlignmentBase):
peak_fit_radius: int, optional
Radius size around the max correlation pixel, for sub-pixel fitting.
Minimum and default value is 1.
low_pass: float or sequence of two floats
Low-pass filter properties, as described in `nabu.misc.fourier_filters`.
high_pass: float or sequence of two floats
High-pass filter properties, as described in `nabu.misc.fourier_filters`.
equispaced_increments: boolean, optional
Tells whether the position increments are equispaced or not. If
equispaced increments are used, we have to compute the correlation
......@@ -598,7 +584,11 @@ class DetectorTranslationAlongBeam(AlignmentBase):
# do correlations
ccs = [
self._compute_correlation_fft(
img_stack[0 if equispaced_increments else ii - 1, ...], img_stack[ii, ...], padding_mode
img_stack[0 if equispaced_increments else ii - 1, ...],
img_stack[ii, ...],
padding_mode,
high_pass=high_pass,
low_pass=low_pass,
)
for ii in range(1, num_imgs)
]
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment