Skip to content
GitLab
Commit 7f5957db authored by myron's avatar myron Committed by Pierre Paleo
Browse files

Ready for MR

parent 3e1056f2
Hide whitespace changes
Inline Side-by-side
nabu/preproc/alignment.py
View file @ 7f5957db
......@@ -8,18 +8,21 @@ from ..misc import fourier_filters
try:
from scipy.ndimage.filters import median_filter
import scipy.fft
local_fftn = scipy.fft.rfftn
local_ifftn = scipy.fft.irfftn
__have_scipy__ = True
except ImportError:
from silx.math.medianfilter import medfilt2d as median_filter
from silx.math.medianfilter import medfilt1d as median_filter_1d
local_fftn = np.fft.fftn
local_ifftn = np.fft.ifftn
__have_scipy__ = False
try:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
__have_matplotlib__ = True
except ImportError:
logging.getLogger(__name__).warning("Matplotlib not available. Plotting disabled")
......@@ -590,7 +593,6 @@ class CenterOfRotation(AlignmentBase):
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, 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])
......@@ -599,30 +601,33 @@ class CenterOfRotation(AlignmentBase):
(f_vals, fv, fh) = self.extract_peak_region_2d(cc, peak_radius=peak_fit_radius, cc_vs=cc_vs, cc_hs=cc_hs)
fitted_shifts_vh = self.refine_max_position_2d(f_vals, fv, fh)
return fitted_shifts_vh[-1] / 2.0
return fitted_shifts_vh[-1] / 2.0
class CenterOfRotationAdaptiveSearch(CenterOfRotation):
""" This adaptive method works by applying a gaussian which highlights, by apodisation, a region
which can possibly contain the good center of rotation.
The whole image is spanned during several applications of the apodisation. At each application
the apodisation function, which is a gaussian, is moved to a new guess position.
The lenght of the step, by which the gaussian is moved, and its sigma are
obtained by multiplying the shortest distance from the left or right border with
a self.step_fraction and self.sigma_fraction factors which ensure global overlapping
"""This adaptive method works by applying a gaussian which highlights, by apodisation, a region
which can possibly contain the good center of rotation.
The whole image is spanned during several applications of the apodisation. At each application
the apodisation function, which is a gaussian, is moved to a new guess position.
The lenght of the step, by which the gaussian is moved, and its sigma are
obtained by multiplying the shortest distance from the left or right border with
a self.step_fraction and self.sigma_fraction factors which ensure global overlapping
"""
sigma_fraction = 1.0/4.0
step_fraction = 1.0/6.0
sigma_fraction = 1.0 / 4.0
step_fraction = 1.0 / 6.0
def find_shift(
self,
img_1: np.ndarray,
img_2: np.ndarray,
roi_yxhw=None,
median_filt_shape=None,
padding_mode=None,
peak_fit_radius=1,
high_pass=None,
low_pass=None,
margins=None,
self,
img_1: np.ndarray,
img_2: np.ndarray,
roi_yxhw=None,
median_filt_shape=None,
padding_mode=None,
peak_fit_radius=1,
high_pass=None,
low_pass=None,
margins=None,
):
"""Find the Center of Rotation (CoR), given two images.
......@@ -672,7 +677,7 @@ class CenterOfRotationAdaptiveSearch(CenterOfRotation):
margins: None or a couple of floats or ints
if margins is None or in the form of (margin1,margin2) the search is done between margin1 and dim_x-1-margin2.
If left to None then by default (margin1,margin2) = ( 10, 10 )
Raises
------
ValueError
......@@ -698,13 +703,13 @@ class CenterOfRotationAdaptiveSearch(CenterOfRotation):
>>> cor_position = CoR_calc.find_shift(radio1, radio2, median_filt_shape=(3, 3), high_pass=20, low_pass=1 )
"""
if peak_fit_radius < 1:
logging.getLogger(__name__).warning(
"Parameter peak_fit_radius should be at least 1, given: %d instead." % peak_fit_radius
)
peak_fit_radius = 1
self._check_img_pair_sizes(img_1, img_2)
used_type = img_1.dtype
......@@ -717,78 +722,85 @@ class CenterOfRotationAdaptiveSearch(CenterOfRotation):
dim_radio = img_1.shape[1]
if margins is None:
lim1, lim2 = 10, dim_radio-1 - 10
lim_1, lim_2 = 10, dim_radio - 1 - 10
else:
lim1, lim2 = margins
lim2 = dim_radio - 1 - lim2
lim_1, lim_2 = margins
lim_2 = dim_radio - 1 - lim_2
if lim1 < 1:
lim1 = 1
if lim2 > dim_radio - 2:
lim2 = dim_radio - 2
if lim_1 < 1:
lim_1 = 1
if lim_2 > dim_radio - 2:
lim_2 = dim_radio - 2
if lim2 <= lim1 :
if lim_2 <= lim_1:
message = (
"Image shape or cropped selection too small for global search. After removal of the margins the search limit collide. The cropped size is %d\n" % (dim_radio)
"Image shape or cropped selection too small for global search. After removal of the margins the search limit collide. The cropped size is %d\n"
% (dim_radio)
)
raise ValueError(message)
found_centers = []
Xcor = lim1
while Xcor < lim2:
x_cor = lim_1
while x_cor < lim_2:
tmp_sigma = (
min(
(img_1.shape[1] - x_cor),
(x_cor),
)
* self.sigma_fraction
)
tmp_sigma = min(
(img_1.shape[1] - Xcor) ,
(Xcor) ,
)* self.sigma_fraction
tmp_x = (np.arange(img_1.shape[1]) - x_cor) / tmp_sigma
apodis = np.exp(-tmp_x * tmp_x / 2.0)
tmpx = (np.arange(img_1.shape[1]) - Xcor) / tmp_sigma
apodis = np.exp(-tmpx * tmpx / 2.0)
Xcor_rel = Xcor - (img_1.shape[1] // 2)
x_cor_rel = x_cor - (img_1.shape[1] // 2)
img_1_apodised = img_1 * apodis
cor_position = super( self.__class__, self).find_shift(
cor_position = super(self.__class__, self).find_shift(
img_1_apodised.astype(used_type),
img_2 .astype(used_type),
img_2.astype(used_type),
low_pass=low_pass,
high_pass=high_pass,
roi_yxhw= roi_yxhw,
roi_yxhw=roi_yxhw,
)
p1 = cor_position*2
if cor_position < 0:
p2 = img_2.shape[1] + cor_position*2
p_1 = cor_position * 2
if cor_position < 0:
p_2 = img_2.shape[1] + cor_position * 2
else:
p2 = -img_2.shape[1] + cor_position*2
if abs(Xcor_rel - p1 / 2) < abs(Xcor_rel - p2 / 2):
cor_position = p1 / 2
p_2 = -img_2.shape[1] + cor_position * 2
if abs(x_cor_rel - p_1 / 2) < abs(x_cor_rel - p_2 / 2):
cor_position = p_1 / 2
else:
cor_position = p2 / 2
cor_position = p_2 / 2
cor_in_img = img_1.shape[1] // 2 + cor_position
tmp_sigma = min(
(img_1.shape[1] - cor_in_img) ,
(cor_in_img) ,
)* self.sigma_fraction
tmp_sigma = (
min(
(img_1.shape[1] - cor_in_img),
(cor_in_img),
)
* self.sigma_fraction
)
M1 = int(round(cor_position + img_1.shape[1] // 2)) - int(round(tmp_sigma))
M2 = int(round(cor_position + img_1.shape[1] // 2)) + int(round(tmp_sigma))
piece1 = img_1[:, M1:M2]
piece2 = img_2[:, img_1.shape[1] - M2 : img_1.shape[1] - M1]
energy = np.array(piece1 * piece1 + piece2 * piece2, "d").sum()
diff_energy = np.array((piece1 - piece2) * (piece1 - piece2), "d").sum()
piece_1 = img_1[:, M1:M2]
piece_2 = img_2[:, img_1.shape[1] - M2 : img_1.shape[1] - M1]
piece_1 = piece_1 - piece_1.mean()
piece_2 = piece_2 - piece_2.mean()
energy = np.array(piece_1 * piece_1 + piece_2 * piece_2, "d").sum()
diff_energy = np.array((piece_1 - piece_2) * (piece_1 - piece_2), "d").sum()
cost = diff_energy / energy
print( cost, abs(Xcor_rel - cor_position), cor_position, energy )
if not np.isnan(cost):
found_centers.append([cost, abs(Xcor_rel - cor_position), cor_position, energy])
found_centers.append([cost, abs(x_cor_rel - cor_position), cor_position, energy])
Xcor = min(Xcor + Xcor* self.step_fraction , Xcor + (dim_radio - Xcor)* self.step_fraction )
x_cor = min(x_cor + x_cor * self.step_fraction, x_cor + (dim_radio - x_cor) * self.step_fraction)
found_centers.sort()
cor_position = found_centers[0][2]
......
nabu/preproc/tests/test_alignment.py
View file @ 7f5957db
......@@ -176,49 +176,47 @@ def _get_challenging_ImsCouple_for_halftomo_cor():
np.random.seed(0)
ny = 400
nx = 500
n_y = 400
n_x = 500
npoints = 2000
poss_x = -nx / 2.0 + 2 * nx * np.random.random(npoints)
poss_y = -ny / 2.0 + 2 * ny * np.random.random(npoints)
poss_x = -n_x / 2.0 + 2 * n_x * np.random.random(npoints)
poss_y = -n_y / 2.0 + 2 * n_y * np.random.random(npoints)
widths = 3 * np.random.random(npoints)
im1 = _peaks2im(poss_y, poss_x, widths, ny, nx)
im1 = _peaks2im(poss_y, poss_x, widths, n_y, n_x)
CHT = 201.67
center_x = 400.123
shift_y = 3.1415
CenterX = 400.123
shiftY = 3.1415
poss_x_2 = 2 * center_x - poss_x
poss_y_2 = shift_y + poss_y
poss_x_2 = 2 * CenterX - poss_x
poss_y_2 = shiftY + poss_y
im2 = _peaks2im(poss_y_2, poss_x_2, widths, ny, nx)
im2 = _peaks2im(poss_y_2, poss_x_2, widths, n_y, n_x)
npoints_spurious = 200
spurious_poss_x = nx * np.random.random(npoints_spurious)
spurious_poss_y = ny * np.random.random(npoints_spurious)
spurious_poss_x = n_x * np.random.random(npoints_spurious)
spurious_poss_y = n_y * np.random.random(npoints_spurious)
widths = [20.0] * npoints_spurious
spurious_im = _peaks2im(spurious_poss_y, spurious_poss_x, widths, ny, nx)
spurious_im = _peaks2im(spurious_poss_y, spurious_poss_x, widths, n_y, n_x)
im1[:] += spurious_im * 0.1
spurious_poss_x = nx * np.random.random(npoints_spurious)
spurious_poss_y = ny * np.random.random(npoints_spurious)
spurious_poss_x = n_x * np.random.random(npoints_spurious)
spurious_poss_y = n_y * np.random.random(npoints_spurious)
widths = [20.0] * npoints_spurious
spurious_im = _peaks2im(spurious_poss_y, spurious_poss_x, widths, ny, nx)
spurious_im = _peaks2im(spurious_poss_y, spurious_poss_x, widths, n_y, n_x)
im2[:] += spurious_im * 0.1
noise_level = 0.2
noise_ima1 = np.random.normal(0.0, size=[ny, nx]) * noise_level
noise_ima2 = np.random.normal(0.0, size=[ny, nx]) * noise_level
noise_ima1 = np.random.normal(0.0, size=[n_y, n_x]) * noise_level
noise_ima2 = np.random.normal(0.0, size=[n_y, n_x]) * noise_level
im1[:] += noise_ima1
im2[:] += noise_ima2
return im1, im2, (CenterX - (nx - 1) / 2.0)
return im1, im2, (center_x - (n_x - 1) / 2.0)
@pytest.mark.usefixtures("bootstrap_base")
......@@ -360,7 +358,7 @@ class TestCor(object):
assert np.isclose(cor_pos, cor_position, atol=self.abs_tol), message
def test_half_tomo_cor_exp(self):
"""test the hal_tomo algorithm on experimental data """
"""test the half_tomo algorithm on experimental data """
radios = nabu_get_data("ha_autocor_radios.npz")
radio1 = radios["radio1"]
......@@ -369,7 +367,6 @@ class TestCor(object):
radio2 = np.fliplr(radio2)
CoR_calc = alignment.CenterOfRotationAdaptiveSearch()
cor_position = CoR_calc.find_shift(radio1, radio2, low_pass=1, high_pass=20)
......@@ -388,20 +385,13 @@ class TestCor(object):
radios = nabu_get_data("ha_autocor_radios.npz")
radio1 = radios["radio1"]
radio2 = radios["radio2"]
cor_pos =radios["cor_pos"]
cor_pos = radios["cor_pos"]
radio2 = np.fliplr(radio2)
CoR_calc = alignment.CenterOfRotationAdaptiveSearch()
CoR_calc = alignment.CenterOfRotation()
cor_position = CoR_calc.find_shift(
radio1,
radio2,
low_pass=1,
high_pass=20,
global_search=((radio1.shape[1] // 2) + cor_pos - 100, (radio1.shape[1] // 2) + cor_pos + 30),
)
cor_position = CoR_calc.find_shift(radio1, radio2, low_pass=1, high_pass=20, margins=(100, 10))
print("Found cor_position", cor_position)
......
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