Skip to content
GitLab
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
tomotools
Nabu
Commits
bb7390cd
Commit
bb7390cd
authored
Sep 28, 2020
by
Pierre Paleo
Browse files
Apply black -l 128 on files
parent
bed0fd0f
Changes
2
Hide whitespace changes
Inline
Side-by-side
nabu/preproc/alignment.py
View file @
bb7390cd
...
...
@@ -38,7 +38,11 @@ except ImportError:
class
AlignmentBase
(
object
):
def
__init__
(
self
,
vert_fft_width
=
False
,
horz_fft_width
=
False
,
verbose
=
False
,
data_type
=
np
.
float32
,
self
,
vert_fft_width
=
False
,
horz_fft_width
=
False
,
verbose
=
False
,
data_type
=
np
.
float32
,
):
"""
Alignment basic functions.
...
...
@@ -89,7 +93,10 @@ class AlignmentBase(object):
raise
ValueError
(
"The same number of images and positions is required."
+
" Shape of stack: %s, shape of positions variable: %s"
%
(
" "
.
join
((
"%d"
%
x
for
x
in
shape_stack
)),
" "
.
join
((
"%d"
%
x
for
x
in
shape_pos
)),)
%
(
" "
.
join
((
"%d"
%
x
for
x
in
shape_stack
)),
" "
.
join
((
"%d"
%
x
for
x
in
shape_pos
)),
)
)
@
staticmethod
...
...
@@ -103,7 +110,10 @@ class AlignmentBase(object):
if
not
np
.
all
(
shape_1
==
shape_2
):
raise
ValueError
(
"Images need to be of the same shape. Shape of image #1: %s, image #2: %s"
%
(
" "
.
join
((
"%d"
%
x
for
x
in
shape_1
)),
" "
.
join
((
"%d"
%
x
for
x
in
shape_2
)),)
%
(
" "
.
join
((
"%d"
%
x
for
x
in
shape_1
)),
" "
.
join
((
"%d"
%
x
for
x
in
shape_2
)),
)
)
@
staticmethod
...
...
@@ -172,7 +182,12 @@ class AlignmentBase(object):
if
np
.
any
(
vertex_yx
<
vertex_min_yx
)
or
np
.
any
(
vertex_yx
>
vertex_max_yx
):
raise
ValueError
(
"Fitted (y: {}, x: {}) positions are outide the input margins y: [{}, {}], and x: [{}, {}]"
.
format
(
vertex_yx
[
0
],
vertex_yx
[
1
],
vertex_min_yx
[
0
],
vertex_max_yx
[
0
],
vertex_min_yx
[
1
],
vertex_max_yx
[
1
],
vertex_yx
[
0
],
vertex_yx
[
1
],
vertex_min_yx
[
0
],
vertex_max_yx
[
0
],
vertex_min_yx
[
1
],
vertex_max_yx
[
1
],
)
)
return
vertex_yx
...
...
@@ -241,7 +256,10 @@ class AlignmentBase(object):
message
=
"Fitted position {} is outide the input margins [{}, {}]"
.
format
(
vertex_x
,
vertex_min_x
,
vertex_max_x
)
else
:
message
=
"Fitted positions outide the input margins [{}, {}]: %d below and %d above"
.
format
(
vertex_min_x
,
vertex_max_x
,
np
.
sum
(
1
-
lower_bound_ok
),
np
.
sum
(
1
-
upper_bound_ok
),
vertex_min_x
,
vertex_max_x
,
np
.
sum
(
1
-
lower_bound_ok
),
np
.
sum
(
1
-
upper_bound_ok
),
)
raise
ValueError
(
message
)
if
return_vertex_val
:
...
...
@@ -362,7 +380,13 @@ class AlignmentBase(object):
return
roi_yxhw
def
_prepare_image
(
self
,
img
,
invalid_val
=
1e-5
,
roi_yxhw
=
None
,
median_filt_shape
=
None
,
low_pass
=
None
,
high_pass
=
None
,
self
,
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.
...
...
@@ -390,7 +414,9 @@ class AlignmentBase(object):
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
()
...
...
@@ -415,7 +441,10 @@ class AlignmentBase(object):
# expanding filter shape with ones, to cover the stack of images
# but disabling inter-image filtering
median_filt_shape
=
np
.
concatenate
(
(
np
.
ones
((
len
(
img_shape
)
-
len
(
median_filt_shape
),),
dtype
=
np
.
int
),
median_filt_shape
,)
(
np
.
ones
((
len
(
img_shape
)
-
len
(
median_filt_shape
),),
dtype
=
np
.
int
),
median_filt_shape
,
)
)
img
=
median_filter
(
img
,
size
=
median_filt_shape
)
else
:
...
...
@@ -510,7 +539,12 @@ class CenterOfRotation(AlignmentBase):
while
Xcor
<
lim2
:
Xcor_rel
=
Xcor
-
(
img_1
.
shape
[
1
]
//
2
)
cor_position
,
merit
,
energy
=
self
.
find_shift
(
img_1
,
img_2
,
low_pass
=
low_pass
,
high_pass
=
high_pass
,
half_tomo_cor_guess
=
Xcor_rel
,
half_tomo_return_stats
=
True
,
img_1
,
img_2
,
low_pass
=
low_pass
,
high_pass
=
high_pass
,
half_tomo_cor_guess
=
Xcor_rel
,
half_tomo_return_stats
=
True
,
)
if
not
np
.
isnan
(
merit
):
found_centers
.
append
([
merit
,
abs
(
Xcor_rel
-
cor_position
),
cor_position
,
energy
])
...
...
@@ -637,7 +671,10 @@ class CenterOfRotation(AlignmentBase):
if
half_tomo_cor_guess
is
not
None
:
cor_in_img
=
img_1
.
shape
[
1
]
//
2
+
half_tomo_cor_guess
tmpsigma
=
min
((
img_1
.
shape
[
1
]
-
cor_in_img
)
/
4.0
,
(
cor_in_img
)
/
4.0
,)
tmpsigma
=
min
(
(
img_1
.
shape
[
1
]
-
cor_in_img
)
/
4.0
,
(
cor_in_img
)
/
4.0
,
)
tmpx
=
(
np
.
arange
(
img_1
.
shape
[
1
])
-
cor_in_img
)
/
tmpsigma
apodis
=
np
.
exp
(
-
tmpx
*
tmpx
/
2.0
)
...
...
@@ -671,7 +708,10 @@ class CenterOfRotation(AlignmentBase):
if
half_tomo_return_stats
:
cor_in_img
=
img_1
.
shape
[
1
]
//
2
+
return_value
tmpsigma
=
min
((
img_1
.
shape
[
1
]
-
cor_in_img
)
/
4.0
,
(
cor_in_img
)
/
4.0
,)
tmpsigma
=
min
(
(
img_1
.
shape
[
1
]
-
cor_in_img
)
/
4.0
,
(
cor_in_img
)
/
4.0
,
)
M1
=
int
(
round
(
return_value
+
img_1
.
shape
[
1
]
//
2
))
-
int
(
round
(
tmpsigma
))
M2
=
int
(
round
(
return_value
+
img_1
.
shape
[
1
]
//
2
))
+
int
(
round
(
tmpsigma
))
...
...
@@ -964,7 +1004,14 @@ class CameraTilt(CenterOfRotation):
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
,
axes
=
(
-
1
,),
high_pass
=
high_pass
,
low_pass
=
low_pass
,)
cc
=
self
.
_compute_correlation_fft
(
img_1
,
img_2
,
padding_mode
,
axes
=
(
-
1
,),
high_pass
=
high_pass
,
low_pass
=
low_pass
,
)
img_shape
=
img_2
.
shape
cc_h_coords
=
np
.
fft
.
fftfreq
(
img_shape
[
-
1
],
1
/
img_shape
[
-
1
])
...
...
@@ -988,7 +1035,10 @@ class CameraTilt(CenterOfRotation):
if
self
.
verbose
:
print
(
"Fitted center of rotation (pixels):"
,
cor_offset_pix
,
"and camera tilt (degrees):"
,
tilt_deg
,
"Fitted center of rotation (pixels):"
,
cor_offset_pix
,
"and camera tilt (degrees):"
,
tilt_deg
,
)
f
,
ax
=
plt
.
subplots
(
1
,
1
)
ax
.
scatter
(
cc_v_coords
,
fitted_shifts_h
)
...
...
@@ -1083,7 +1133,11 @@ class CameraFocus(CenterOfRotation):
roi_yxhw
=
self
.
_determine_roi
(
img_shape
,
roi_yxhw
)
img_stack
=
self
.
_prepare_image
(
img_stack
,
roi_yxhw
=
roi_yxhw
,
median_filt_shape
=
median_filt_shape
,
low_pass
=
low_pass
,
high_pass
=
high_pass
,
img_stack
,
roi_yxhw
=
roi_yxhw
,
median_filt_shape
=
median_filt_shape
,
low_pass
=
low_pass
,
high_pass
=
high_pass
,
)
img_stds
=
np
.
std
(
img_stack
,
axis
=
(
-
2
,
-
1
))
/
np
.
mean
(
img_stack
,
axis
=
(
-
2
,
-
1
))
...
...
@@ -1101,7 +1155,10 @@ class CameraFocus(CenterOfRotation):
if
self
.
verbose
:
print
(
"Fitted focus motor position:"
,
focus_pos
,
"and corresponding image position:"
,
focus_ind
,
"Fitted focus motor position:"
,
focus_pos
,
"and corresponding image position:"
,
focus_ind
,
)
f
,
ax
=
plt
.
subplots
(
1
,
1
)
ax
.
stem
(
img_pos
,
img_stds
)
...
...
@@ -1282,12 +1339,19 @@ class CameraFocus(CenterOfRotation):
self
.
_check_img_block_size
(
roi_yxhw
[
2
:],
regions_number
,
suggest_new_shape
=
False
)
img_stack
=
self
.
_prepare_image
(
img_stack
,
roi_yxhw
=
roi_yxhw
,
median_filt_shape
=
median_filt_shape
,
low_pass
=
low_pass
,
high_pass
=
high_pass
,
img_stack
,
roi_yxhw
=
roi_yxhw
,
median_filt_shape
=
median_filt_shape
,
low_pass
=
low_pass
,
high_pass
=
high_pass
,
)
img_shape
=
img_stack
.
shape
[
-
2
:]
block_size
=
np
.
array
(
img_shape
)
/
regions_number
block_stack_size
=
np
.
array
([
num_imgs
,
regions_number
,
block_size
[
-
2
],
regions_number
,
block_size
[
-
1
]],
dtype
=
np
.
int
,)
block_stack_size
=
np
.
array
(
[
num_imgs
,
regions_number
,
block_size
[
-
2
],
regions_number
,
block_size
[
-
1
]],
dtype
=
np
.
int
,
)
img_stack
=
np
.
reshape
(
img_stack
,
block_stack_size
)
img_stds
=
np
.
std
(
img_stack
,
axis
=
(
-
3
,
-
1
))
/
np
.
mean
(
img_stack
,
axis
=
(
-
3
,
-
1
))
...
...
@@ -1315,7 +1379,10 @@ class CameraFocus(CenterOfRotation):
if
self
.
verbose
:
print
(
"Fitted focus motor position:"
,
focus_pos
,
"and corresponding image position:"
,
focus_ind
,
"Fitted focus motor position:"
,
focus_pos
,
"and corresponding image position:"
,
focus_ind
,
)
print
(
"Fitted tilts (to be divided by pixel size, and converted to deg): (v, h) %s"
%
tilts_vh
)
fig
=
plt
.
figure
()
...
...
@@ -1324,10 +1391,16 @@ class CameraFocus(CenterOfRotation):
regions_half_shape
=
(
regions_number
-
1
)
/
2
base_points
=
np
.
linspace
(
-
regions_half_shape
,
regions_half_shape
,
regions_number
)
ax
.
plot
(
np
.
zeros
((
regions_number
,)),
base_points
,
np
.
polyval
([
tg_v
,
focus_pos
],
base_points
),
"C2"
,
np
.
zeros
((
regions_number
,)),
base_points
,
np
.
polyval
([
tg_v
,
focus_pos
],
base_points
),
"C2"
,
)
ax
.
plot
(
base_points
,
np
.
zeros
((
regions_number
,)),
np
.
polyval
([
tg_h
,
focus_pos
],
base_points
),
"C2"
,
base_points
,
np
.
zeros
((
regions_number
,)),
np
.
polyval
([
tg_h
,
focus_pos
],
base_points
),
"C2"
,
)
ax
.
scatter
(
0
,
0
,
focus_pos
,
marker
=
"o"
,
c
=
"C1"
)
ax
.
set_title
(
"Images std"
)
...
...
nabu/preproc/tests/test_alignment.py
View file @
bb7390cd
...
...
@@ -48,8 +48,18 @@ def bootstrap_fcs(request):
cls
.
abs_tol_dist
=
1e-2
cls
.
abs_tol_tilt
=
2.5e-4
(
cls
.
data
,
cls
.
img_pos
,
cls
.
pixel_size
,
(
calib_data_std
,
calib_data_angle
),)
=
get_focus_data
(
"test_alignment_focus.h5"
)
(
cls
.
angle_best_ind
,
cls
.
angle_best_pos
,
cls
.
angle_tilt_v
,
cls
.
angle_tilt_h
,)
=
calib_data_angle
(
cls
.
data
,
cls
.
img_pos
,
cls
.
pixel_size
,
(
calib_data_std
,
calib_data_angle
),
)
=
get_focus_data
(
"test_alignment_focus.h5"
)
(
cls
.
angle_best_ind
,
cls
.
angle_best_pos
,
cls
.
angle_tilt_v
,
cls
.
angle_tilt_h
,
)
=
calib_data_angle
cls
.
std_best_ind
,
cls
.
std_best_pos
=
calib_data_std
...
...
@@ -73,8 +83,7 @@ def get_cor_data_h5(*dataset_path):
def
get_cor_data_half_tomo
():
""" Obtains two weakly overlapping images with features plus spurious noise for a challenging test of cor retrieval for half tomo.
"""
"""Obtains two weakly overlapping images with features plus spurious noise for a challenging test of cor retrieval for half tomo."""
datasource
=
ExternalResources
(
project
=
"nabu"
,
url_base
=
None
)
myfile
=
os
.
path
.
join
(
datasource
.
data_home
,
"data_for_ht_cor.h5"
)
...
...
@@ -130,7 +139,10 @@ def get_focus_data(*dataset_path):
img_pos
=
hf
[
"/entry/instrument/detector/distance"
][()]
pixel_size
=
np
.
mean
(
[
hf
[
"/entry/instrument/detector/x_pixel_size"
][()],
hf
[
"/entry/instrument/detector/y_pixel_size"
][()],]
[
hf
[
"/entry/instrument/detector/x_pixel_size"
][()],
hf
[
"/entry/instrument/detector/y_pixel_size"
][()],
]
)
angle_best_ind
=
hf
[
"/calibration/focus/angle/best_img"
][()]
...
...
@@ -147,8 +159,7 @@ def get_focus_data(*dataset_path):
def
_peaks2im
(
possY
,
possX
,
widths
,
NY
,
NX
):
""" Given a set of positions, widths and the image shape creates an image with spots
"""
"""Given a set of positions, widths and the image shape creates an image with spots"""
res
=
np
.
zeros
([
NY
,
NX
],
"f"
)
ys
=
np
.
arange
(
NY
)
xs
=
np
.
arange
(
NX
)
...
...
@@ -161,8 +172,7 @@ def _peaks2im(possY, possX, widths, NY, NX):
def
_get_challenging_ImsCouple_for_halftomo_cor
():
""" Creates two weakly overlapping images with features plus spurious noise for a challenging test of cor retrieval for half tomo.
"""
"""Creates two weakly overlapping images with features plus spurious noise for a challenging test of cor retrieval for half tomo."""
np
.
random
.
seed
(
0
)
...
...
@@ -258,7 +268,10 @@ class TestAlignmentBase(object):
cc_coords
=
np
.
arange
(
0
,
8
)
(
found_peaks_val
,
found_peaks_pos
,)
=
alignment
.
AlignmentBase
.
extract_peak_regions_1d
(
img
,
axis
=-
1
,
cc_coords
=
cc_coords
)
(
found_peaks_val
,
found_peaks_pos
,
)
=
alignment
.
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
,
:],
...
...
@@ -293,8 +306,7 @@ class TestCor(object):
@
pytest
.
mark
.
skipif
(
not
(
__has_scipy__
),
reason
=
"need scipy for this test"
)
def
test_noisyHF_cor_posx
(
self
):
""" test with noise at high frequencies
"""
"""test with noise at high frequencies"""
radio1
=
self
.
data
[
0
,
:,
:]
radio2
=
np
.
fliplr
(
self
.
data
[
1
,
:,
:])
...
...
@@ -317,8 +329,7 @@ class TestCor(object):
assert
np
.
isclose
(
self
.
cor_gl_pix
,
cor_position
,
atol
=
self
.
abs_tol
),
message
def
test_half_tomo_cor_Err
(
self
):
""" test that the test is challenging enough
"""
"""test that the test is challenging enough"""
radio1
=
self
.
ht_im1
radio2
=
np
.
fliplr
(
self
.
ht_im2
)
cor_pos
=
self
.
ht_cor
...
...
@@ -334,8 +345,7 @@ class TestCor(object):
assert
not
np
.
isclose
(
cor_pos
,
cor_position
,
atol
=
self
.
abs_tol
),
message
def
test_half_tomo_cor
(
self
):
""" test that the half tomo algorithm is precise enough
"""
"""test that the half tomo algorithm is precise enough"""
radio1
=
self
.
ht_im1
radio2
=
np
.
fliplr
(
self
.
ht_im2
)
cor_pos
=
self
.
ht_cor
...
...
@@ -350,8 +360,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 and global search
"""
"""test the hal_tomo algorithm on experimental data and global search"""
radios
=
nabu_get_data
(
"ha_autocor_radios.npz"
)
radio1
=
radios
[
"radio1"
]
...
...
@@ -374,8 +383,7 @@ class TestCor(object):
assert
np
.
isclose
(
cor_pos
,
cor_position
,
atol
=
self
.
abs_tol
),
message
def
test_half_tomo_cor_exp_limited
(
self
):
""" test the hal_tomo algorithm on experimental data and global search with limits
"""
"""test the hal_tomo algorithm on experimental data and global search with limits"""
radios
=
nabu_get_data
(
"ha_autocor_radios.npz"
)
radio1
=
radios
[
"radio1"
]
...
...
@@ -478,7 +486,10 @@ class TestDetectorTranslation(object):
)
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
.
all_shifts_vh
,)
message
=
"Computed shifts %s and expected %s do not coincide"
%
(
found_shifts_list
,
self
.
all_shifts_vh
,
)
assert
np
.
all
(
np
.
isclose
(
found_shifts_list
,
self
.
all_shifts_vh
,
atol
=
self
.
abs_tol
)),
message
@
pytest
.
mark
.
skipif
(
not
(
__has_scipy__
),
reason
=
"need scipy for this test"
)
...
...
@@ -498,7 +509,10 @@ class TestDetectorTranslation(object):
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
),)
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
...
...
@@ -527,7 +541,10 @@ class TestFocus(object):
expected_tilts_vh
=
np
.
squeeze
(
np
.
array
([
self
.
angle_tilt_v
,
self
.
angle_tilt_h
]))
computed_tilts_vh
=
-
tilts_vh
/
(
self
.
pixel_size
/
1000
)
message
=
"Computed tilts %s and expected %s do not coincide"
%
(
computed_tilts_vh
,
expected_tilts_vh
,)
message
=
"Computed tilts %s and expected %s do not coincide"
%
(
computed_tilts_vh
,
expected_tilts_vh
,
)
assert
np
.
all
(
np
.
isclose
(
computed_tilts_vh
,
expected_tilts_vh
,
atol
=
self
.
abs_tol_tilt
)),
message
def
test_size_determination
(
self
):
...
...
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment