Commit 0535ee49 authored by Pierre Paleo's avatar Pierre Paleo
Browse files

Merge branch 'translate_alignxc' into 'master'

Shift finding for translation stage along beam

See merge request !41
parents ca596003 118f2120
Pipeline #28227 passed with stages
in 5 minutes and 15 seconds
nabu.misc.fourier_filters module
========================
.. automodule:: nabu.misc.fourier_filters
:members:
:undoc-members:
:show-inheritance:
# -*- coding: utf-8 -*-
"""
Fourier filters.
"""
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=None, use_rfft=False, data_type=np.float64):
"""Computes a low pass filter using the erfc function.
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
use_rfft: boolean, optional
Creates a filter to be used with the result of a rfft type of Fourier transform. Defaults to False.
data_type: `numpy.dtype`, optional
Specifies the data type of the computed filter. It defaults to `numpy.float64`
Raises
------
ValueError
In case of malformed cutoff_par
Returns
-------
numpy.array_like
The computed filter
"""
if cutoff_par is None:
return 1
elif 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")
r = np.sqrt(np.sum(np.array(coords, dtype=data_type) ** 2, axis=0))
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.exp(-(np.pi ** 2) * (r ** 2) * (cutoff_pix ** 2) * 2)
# Making sure to force result to chosen data type
res = res.astype(data_type)
if use_rfft:
slicelist = [slice(None)] * (len(res.shape) - 1) + [slice(0, res.shape[-1] // 2 + 1)]
return res[tuple(slicelist)]
else:
return res
def get_highpass_filter(img_shape, cutoff_par=None, use_rfft=False, data_type=np.float64):
"""Computes a high pass filter using the erfc function.
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
use_rfft: boolean, optional
Creates a filter to be used with the result of a rfft type of Fourier transform. Defaults to False.
data_type: `numpy.dtype`, optional
Specifies the data type of the computed filter. It defaults to `numpy.float64`
Raises
------
ValueError
In case of malformed cutoff_par
Returns
-------
numpy.array_like
The computed filter
"""
if cutoff_par is None:
return 1
else:
return 1 - get_lowpass_filter(img_shape, cutoff_par, use_rfft=use_rfft, data_type=data_type)
def get_bandpass_filter(img_shape, cutoff_lowpass=None, cutoff_highpass=None, use_rfft=False, data_type=np.float64):
"""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
use_rfft: boolean, optional
Creates a filter to be used with the result of a rfft type of Fourier transform. Defaults to False.
data_type: `numpy.dtype`, optional
Specifies the data type of the computed filter. It defaults to `numpy.float64`
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, use_rfft=use_rfft, data_type=data_type
) * get_highpass_filter(img_shape, cutoff_par=cutoff_highpass, use_rfft=use_rfft, data_type=data_type)
This diff is collapsed.
import pytest
import numpy as np
import os
import h5py
import scipy.ndimage
from silx.third_party.EdfFile import EdfFile
from nabu.preproc.alignment import CenterOfRotation, DetectorTranslationAlongBeam, AlignmentBase
from nabu.testutils import utilstest
@pytest.fixture(scope="class")
def bootstrap_base(request):
cls = request.cls
cls.abs_tol = 2.5e-2
@pytest.fixture(scope="class")
def bootstrap_cor(request):
cls = request.cls
cls.abs_tol = 0.2
cls.data, cls.px = get_data_h5("tworadios.h5")
@pytest.fixture(scope="class")
def bootstrap_dtr(request):
cls = request.cls
cls.abs_tol = 1e-1
# loading alignxc edf files from Christian. The last one is the dark. The first 6 are images for translation 0,1...6
images = np.array([EdfFile(utilstest.getfile("alignxc%04d.edf" % i)).GetData(0, DataType="FloatValue") for i in range(7)])
align_imgs, dark_img = images[:-1], images[-1]
# removing dark
cls.align_images = align_imgs - dark_img
cls.img_pos = (1 + np.arange(6)) * 0.01
cls.expected_shifts_vh = np.array((129.93, 353.18))
cls.reference_shifts_list = [
[-9.39, 11.29],
[-5.02, 3.81],
[-3.67, 9.73],
[-4.72, 16.71],
[6.02, 20.28],
]
def get_data_h5(*dataset_path):
"""
Get a dataset file from silx.org/pub/nabu/data
dataset_args is a list describing a nested folder structures, ex.
["path", "to", "my", "dataset.h5"]
"""
dataset_relpath = os.path.join(*dataset_path)
dataset_downloaded_path = utilstest.getfile(dataset_relpath)
with h5py.File(dataset_downloaded_path, "r") as hf:
nxentry = "entry/instrument/detector"
px = hf[nxentry + "/x_rotation_axis_pixel_position"][()]
data = hf[nxentry + "/data"][()]
return data, px
@pytest.mark.usefixtures("bootstrap_base")
class TestAlignmentBase(object):
def test_peak_fitting_2d_3x3(self):
# Fit a 3 x 3 grid
fy = np.linspace(-1, 1, 3)
fx = np.linspace(-1, 1, 3)
yy, xx = np.meshgrid(fy, fx, indexing="ij")
peak_pos_yx = np.random.rand(2) * 1.6 - 0.8
f_vals = np.exp(-((yy - peak_pos_yx[0]) ** 2 + (xx - peak_pos_yx[1]) ** 2) / 100)
fitted_peak_pos_yx = AlignmentBase.refine_max_position_2d(f_vals, fy, fx)
message = (
"Computed peak position: (%f, %f) " % (*fitted_peak_pos_yx,)
+ " and real peak position (%f, %f) do not coincide." % (*peak_pos_yx,)
+ " Difference: (%f, %f)," % (*(fitted_peak_pos_yx - peak_pos_yx),)
+ " tolerance: %f" % self.abs_tol
)
assert np.all(np.isclose(peak_pos_yx, fitted_peak_pos_yx, atol=self.abs_tol)), message
def test_peak_fitting_2d_error_checking(self):
# Fit a 3 x 3 grid
fy = np.linspace(-1, 1, 3)
fx = np.linspace(-1, 1, 3)
yy, xx = np.meshgrid(fy, fx, indexing="ij")
peak_pos_yx = np.random.rand(2) + 1.5
f_vals = np.exp(-((yy - peak_pos_yx[0]) ** 2 + (xx - peak_pos_yx[1]) ** 2) / 100)
with pytest.raises(ValueError) as ex:
AlignmentBase.refine_max_position_2d(f_vals, fy, fx)
message = (
"Error should have been raised about the peak being fitted outside margins, "
+ "other error raised instead:\n%s" % str(ex.value)
)
assert "positions are outide the input margins" in str(ex.value), message
def test_extract_peak_regions_1d(self):
img = np.random.randint(0, 10, size=(8, 8))
peaks_pos = np.argmax(img, axis=-1)
peaks_val = np.max(img, axis=-1)
cc_coords = np.arange(0, 8)
(found_peaks_val, found_peaks_pos) = AlignmentBase.extract_peak_regions_1d(img, axis=-1, cc_coords=cc_coords)
message = "The found peak positions do not correspond to the expected peak positions:\n Expected: %s\n Found: %s" % (
peaks_pos,
found_peaks_pos[1, :],
)
assert np.all(peaks_val == found_peaks_val[1, :]), message
@pytest.mark.usefixtures("bootstrap_cor")
class TestCor(object):
def test_cor_posx(self):
radio1 = self.data[0, :, :]
radio2 = np.fliplr(self.data[1, :, :])
CoR_calc = CenterOfRotation()
cor_position = CoR_calc.find_shift(radio1, radio2)
message = "Computed CoR %f " % cor_position + " and real CoR %f do not coincide" % self.px
assert np.abs(self.px - cor_position) < self.abs_tol, message
def test_noisy_cor_posx(self):
radio1 = np.fmax(self.data[0, :, :], 0)
radio2 = np.fmax(self.data[1, :, :], 0)
radio1 = np.random.poisson(radio1 * 400)
radio2 = np.random.poisson(np.fliplr(radio2) * 400)
CoR_calc = CenterOfRotation()
cor_position = CoR_calc.find_shift(radio1, radio2, median_filt_shape=(3, 3))
message = "Computed CoR %f " % cor_position + " and real CoR %f do not coincide" % self.px
assert np.abs(self.px - cor_position) < self.abs_tol, message
def test_noisyHF_cor_posx(self):
""" test with noise at high frequencies
"""
radio1 = self.data[0, :, :]
radio2 = np.fliplr(self.data[1, :, :])
noise_level = radio1.max() / 16.0
noise_ima1 = np.random.normal(0.0, size=radio1.shape) * noise_level
noise_ima2 = np.random.normal(0.0, size=radio2.shape) * noise_level
noise_ima1 = noise_ima1 - scipy.ndimage.filters.gaussian_filter(noise_ima1, 2.0)
noise_ima2 = noise_ima2 - scipy.ndimage.filters.gaussian_filter(noise_ima2, 2.0)
radio1 = radio1 + noise_ima1
radio2 = radio2 + noise_ima2
CoR_calc = CenterOfRotation()
# cor_position = CoR_calc.find_shift(radio1, radio2)
cor_position = CoR_calc.find_shift(radio1, radio2, low_pass=(6.0, 0.3))
message = "Computed CoR %f " % cor_position + " and real CoR %f do not coincide" % self.px
assert np.abs(self.px - cor_position) < self.abs_tol, message
def test_cor_posx_linear(self):
radio1 = self.data[0, :, :]
radio2 = np.fliplr(self.data[1, :, :])
CoR_calc = CenterOfRotation()
cor_position = CoR_calc.find_shift(radio1, radio2, padding_mode="constant")
message = "Computed CoR %f " % cor_position + " and real CoR %f do not coincide" % self.px
assert np.abs(self.px - cor_position) < self.abs_tol, message
def test_error_checking_001(self):
CoR_calc = CenterOfRotation()
radio1 = self.data[0, :, :1:]
radio2 = self.data[1, :, :]
with pytest.raises(ValueError) as ex:
CoR_calc.find_shift(radio1, radio2)
message = "Error should have been raised about img #1 shape, other error raised instead:\n%s" % str(ex.value)
assert "Images need to be 2-dimensional. Shape of image #1" in str(ex.value), message
def test_error_checking_002(self):
CoR_calc = CenterOfRotation()
radio1 = self.data[0, :, :]
radio2 = self.data
with pytest.raises(ValueError) as ex:
CoR_calc.find_shift(radio1, radio2)
message = "Error should have been raised about img #2 shape, other error raised instead:\n%s" % str(ex.value)
assert "Images need to be 2-dimensional. Shape of image #2" in str(ex.value), message
def test_error_checking_003(self):
CoR_calc = CenterOfRotation()
radio1 = self.data[0, :, :]
radio2 = self.data[1, :, 0:10]
with pytest.raises(ValueError) as ex:
CoR_calc.find_shift(radio1, radio2)
message = "Error should have been raised about different image shapes, " + "other error raised instead:\n%s" % str(
ex.value
)
assert "Images need to be of the same shape" in str(ex.value), message
@pytest.mark.usefixtures("bootstrap_dtr")
class TestDetectorTranslation(object):
def test_alignxc(self):
T_calc = DetectorTranslationAlongBeam()
shifts_v, shifts_h, found_shifts_list = T_calc.find_shift(self.align_images, self.img_pos, return_shifts=True)
message = "Computed shifts coefficients %s and expected %s do not coincide" % (
(shifts_v, shifts_h),
self.expected_shifts_vh,
)
assert np.all(np.isclose(self.expected_shifts_vh, [shifts_v, shifts_h], atol=self.abs_tol)), message
message = "Computed shifts %s and expected %s do not coincide" % (found_shifts_list, self.reference_shifts_list)
assert np.all(np.isclose(found_shifts_list, self.reference_shifts_list, atol=self.abs_tol)), message
def test_alignxc_synth(self):
T_calc = DetectorTranslationAlongBeam()
stack = np.zeros([4, 512, 512], "d")
for i in range(4):
stack[i, 200 - i * 10, 200 - i * 10] = 1
stack = scipy.ndimage.filters.gaussian_filter(stack, [0, 10, 10.0]) * 100
x, y = np.meshgrid(np.arange(stack.shape[-1]), np.arange(stack.shape[-2]))
for i in range(4):
xc = x - (250 + i * 1.234)
yc = y - (250 + i * 1.234 * 2)
stack[i] += np.exp(-(xc * xc + yc * yc) * 0.5)
shifts_v, shifts_h, found_shifts_list = T_calc.find_shift(
stack, np.array([0.0, 1, 2, 3]), high_pass=1.0, return_shifts=True
)
message = "Found shifts per units %s and reference %s do not coincide" % ((shifts_v, shifts_h), (-1.234 * 2, -1.234))
assert np.all(np.isclose((shifts_v, shifts_h), (-1.234 * 2, -1.234), atol=self.abs_tol)), message
import pytest
import numpy as np
import os
import h5py
from nabu.preproc.alignment import CenterOfRotation
from nabu.testutils import utilstest
@pytest.fixture(scope='class')
def bootstrap(request):
cls = request.cls
cls.data, cls.px = get_data_h5("tworadios.h5")
cls.cor_tol = 0.2
def get_data_h5(*dataset_path):
"""
Get a dataset file from silx.org/pub/nabu/data
dataset_args is a list describing a nested folder structures, ex.
["path", "to", "my", "dataset.h5"]
"""
dataset_relpath = os.path.join(*dataset_path)
dataset_downloaded_path = utilstest.getfile(dataset_relpath)
with h5py.File(dataset_downloaded_path, 'r') as hf:
nxentry = 'entry/instrument/detector'
px = hf[nxentry+'/x_rotation_axis_pixel_position'][()]
data = hf[nxentry+'/data'][()]
return data, px
@pytest.mark.usefixtures('bootstrap')
class TestCor:
def test_cor_posx(self):
radio1 = self.data[0, :, :]
radio2 = np.fliplr(self.data[1, :, :])
CoR_calc = CenterOfRotation(poly_deg=2)
cor_position = CoR_calc.find_shift(radio1, radio2)
message = "Computed CoR %f " % cor_position + " and real CoR %f do not coincide" % self.px
assert np.abs(self.px - cor_position) < self.cor_tol, message
def test_noisy_cor_posx(self):
radio1 = np.fmax(self.data[0, :, :], 0)
radio2 = np.fmax(self.data[1, :, :], 0)
radio1 = np.random.poisson(radio1 * 400)
radio2 = np.random.poisson(np.fliplr(radio2) * 400)
CoR_calc = CenterOfRotation(poly_deg=2)
cor_position = CoR_calc.find_shift(radio1, radio2, median_filt_shape=(3, 3))
message = "Computed CoR %f " % cor_position + " and real CoR %f do not coincide" % self.px
assert np.abs(self.px - cor_position) < self.cor_tol, message
def test_error_checking_001(self):
CoR_calc = CenterOfRotation(poly_deg=2)
radio1 = self.data[0, :, :1:]
radio2 = self.data[1, :, :]
with pytest.raises(ValueError) as ex:
CoR_calc.find_shift(radio1, radio2)
message = 'Error should have been raised about img #1 shape, other error raised instead:\n%s' % str(ex.value)
assert 'Images need to be 2-dimensional. Shape of image #1' in str(ex.value), message
def test_error_checking_002(self):
CoR_calc = CenterOfRotation(poly_deg=2)
radio1 = self.data[0, :, :]
radio2 = self.data
with pytest.raises(ValueError) as ex:
CoR_calc.find_shift(radio1, radio2)
message = 'Error should have been raised about img #2 shape, other error raised instead:\n%s' % str(ex.value)
assert 'Images need to be 2-dimensional. Shape of image #2' in str(ex.value), message
def test_error_checking_003(self):
CoR_calc = CenterOfRotation(poly_deg=2)
radio1 = self.data[0, :, :]
radio2 = self.data[1, :, 0:10]
with pytest.raises(ValueError) as ex:
CoR_calc.find_shift(radio1, radio2)
message = (
'Error should have been raised about different image shapes, '
+ 'other error raised instead:\n%s' % str(ex.value))
assert 'Images need to be of the same shape' in str(ex.value), 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