Commit cc030af1 authored by myron's avatar myron
Browse files

Restructured filtering chain and adopted more understandable parameters. Documentation.

parent 0a254fb7
Pipeline #27538 passed with stages
in 2 minutes and 32 seconds
......@@ -16,32 +16,84 @@ except ImportError:
__have_scipy__ = False
def _get_lowpass_filter(img_shape, cutoff_pix, cutoff_trans_pix):
def _get_lowpass_filter(img_shape, cutoff_par):
"""Computes a low pass filter with the erfc.
:param img_shape: Shape of the image
:type img_shape: tuple
:param cutoff_pix: Position of the cut off in pixels
:type cutoff_pix: float
:param cutoff_trans_pix: Size of the cut off transition in pixels
:type cutoff_trans_pix: float
:return: The computes filter
:rtype: `numpy.array_like`
Parameters
----------
img_shape: tuple
Shape of the image
cutoff_par: 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 while a smooth erfc transition to zero is done
Raises
------
ValueError
In case cutoff_par is not properly given
Returns
-------
numpy.array_like
The computed filter
"""
if isinstance(cutoff_par, (int, float)):
cutoff_pix = cutoff_par
cutoff_trans_fact = None
else:
try:
cutoff_pix, cutoff_trans_fact = cutoff_par
except ValueError:
raise ValueError("Argument cutoff_par (which specifies the pass filter shape) must be either a scalar or a"
" sequence of two scalars")
if (not isinstance(cutoff_pix, (int, float))) or (not isinstance(cutoff_trans_fact, (int, float))):
raise ValueError("Argument cutoff_par (which specifies the pass filter shape) must be one number or a sequence"
"of two numbers")
coords = [np.fft.fftfreq(s, 1) for s in img_shape]
coords = np.meshgrid(*coords, indexing="ij")
k_cut = 0.5 / cutoff_pix
k_cut_width = 0.5 / cutoff_trans_pix
r = np.sqrt(np.sum(np.array(coords) ** 2, axis=0))
k_pos_rescaled = (r - k_cut) / k_cut_width
if __have_scipy__:
res = spspe.erfc(k_pos_rescaled) / 2
if cutoff_trans_fact is not None:
k_cut = 0.5 / cutoff_pix
k_cut_width = k_cut * cutoff_trans_fact
k_pos_rescaled = (r - k_cut) / k_cut_width
if __have_scipy__:
res = spspe.erfc(k_pos_rescaled) / 2
else:
res = np.array(list(map(math.erfc, k_pos_rescaled))) / 2
else:
res = np.array(list(map(math.erfc, k_pos_rescaled))) / 2
res = np.exp(- np.pi*np.pi*r*r*cutoff_pix*cutoff_pix*2)
return res
def _get_highpass_filter(img_shape, cutoff_pix, cutoff_par):
"""Computes a high pass filter with the erfc.
Parameters
----------
img_shape: tuple
Shape of the image
cutoff_par: 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 and then a smooth transition to zero is done for smaller frequency
Returns
-------
numpy.array_like
The computed filter
"""
res = 1 - _get_lowpass_filter(img_shape, cutoff_pix, cutoff_par)
return res
......@@ -188,14 +240,46 @@ class AlignmentBase(object):
@staticmethod
def _prepare_image(
img, invalid_val=1e-5, roi_yxhw=None, median_filt_shape=None, cutoff_pix=None, cutoff_trans_pix=None,
img, invalid_val=1e-5, roi_yxhw=None, median_filt_shape=None, high_pass=None, low_pass=None
):
"""Prepare and returns a cropped and filtered image, or array of filtered images if the input is an array of images.
Parameters
----------
img: numpy.ndarray
image or stack of images
invalid_val: float
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
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
Returns
-------
numpy.array_like
The computed filter
"""
img = np.squeeze(img) # Removes singleton dimensions, but does a shallow copy
img = np.ascontiguousarray(img)
if roi_yxhw is not None:
img = img[
..., roi_yxhw[0] : roi_yxhw[0] + roi_yxhw[2], roi_yxhw[1] : roi_yxhw[1] + roi_yxhw[3],
..., roi_yxhw[0]: roi_yxhw[0] + roi_yxhw[2], roi_yxhw[1]: roi_yxhw[1] + roi_yxhw[3],
]
img = img.copy()
......@@ -203,8 +287,12 @@ class AlignmentBase(object):
img[np.isnan(img)] = invalid_val
img[np.isinf(img)] = invalid_val
if cutoff_pix is not None and cutoff_trans_pix is not None:
myfilter = _get_lowpass_filter(img.shape[-2:], cutoff_pix, cutoff_trans_pix)
if high_pass is not None or low_pass is not None:
myfilter = np.ones_like(img)
if low_pass is not None:
myfilter[:] *= _get_lowpass_filter(img.shape[-2:], low_pass)
if high_pass is not None:
myfilter[:] *= _get_highpass_filter(img.shape[-2:], high_pass)
img_shape = img.shape
if len(img_shape) == 2:
img = np.fft.ifft2(np.fft.fft2(img) * myfilter).real
......@@ -252,7 +340,7 @@ class AlignmentBase(object):
if not do_circular_conv:
cc = np.fft.fftshift(cc, axes=(-2, -1))
cc = cc[padding[0] : -padding[0], padding[1] : -padding[1]]
cc = cc[padding[0]: -padding[0], padding[1]: -padding[1]]
cc = np.fft.ifftshift(cc, axes=(-2, -1))
return cc
......@@ -302,8 +390,8 @@ class CenterOfRotation(AlignmentBase):
median_filt_shape=None,
padding_mode=None,
peak_fit_radius=1,
cutoff_pix=None,
cutoff_trans_pix=None,
high_pass=None,
low_pass=None
):
"""Find the Center of Rotation (CoR), given to images.
......@@ -346,6 +434,21 @@ class CenterOfRotation(AlignmentBase):
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
Raises
------
ValueError
......@@ -386,15 +489,15 @@ class CenterOfRotation(AlignmentBase):
img_1,
roi_yxhw=roi_yxhw,
median_filt_shape=median_filt_shape,
cutoff_pix=cutoff_pix,
cutoff_trans_pix=cutoff_trans_pix,
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,
cutoff_pix=cutoff_pix,
cutoff_trans_pix=cutoff_trans_pix,
high_pass=high_pass,
low_pass=low_pass
)
cc = self._compute_correlation_fft(img_1, img_2, padding_mode)
......
......@@ -81,7 +81,7 @@ class TestCor:
# cor_position = CoR_calc.find_shift(radio1, radio2)
cor_position = CoR_calc.find_shift(
radio1, radio2, cutoff_pix=6.0, cutoff_trans_pix=20.0
radio1, radio2, low_pass=(6.0, 0.3)
)
message = (
......
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