Commit a99e3a84 authored by Nicola Vigano's avatar Nicola Vigano
Browse files

Alignment: adapted 1d max fitting to the same style as the 2d version



Only difference: here we can use Polynomial.fit from numpy.

Signed-off-by: Nicola Vigano's avatarNicola VIGANÒ <nicola.vigano@esrf.fr>
parent e9022e41
import numpy as np
import math
from numpy.polynomial.polynomial import polyroots
from numpy.polynomial.polynomial import Polynomial
from nabu.utils import previouspow2
......@@ -117,50 +117,59 @@ class AlignmentBase(object):
return vertex_yx
@staticmethod
def refine_max_position_1d(fx, f_vals):
def refine_max_position_1d(f_vals, fx=None):
"""Computes the sub-pixel max position of the given function sampling.
Parameters
----------
fx: numpy.ndarray
Coordinates of the sampled points
f_vals: numpy.ndarray
Function values of the sampled points
fx: numpy.ndarray, optional
Coordinates of the sampled points
Raises
------
ValueError
In case position and values do not have the same size.
In case position and values do not have the same size, or in case
the fitted maximum is outside the fitting region.
Returns
-------
float
Estimated function max, according to the coordinates in fx.
"""
npix = fx.size
if not npix == f_vals.size:
if not len(f_vals.shape) == 1:
raise ValueError(
"Coordinates and values should have the same length. Sizes of fx: %d, f_vals: %d" % (fx.size, f_vals.size)
"The fitted values should form a 1-dimensional array. Array of shape: [%s] was given."
% (" ".join(("%d" % s for s in f_vals.shape)))
)
if fx is None:
fx_half_size = (f_vals.size - 1) / 2
fx = np.linspace(-fx_half_size, fx_half_size, f_vals.size)
elif not (len(fx.shape) == 1 and np.all(fx.size == f_vals.size)):
raise ValueError(
"Base coordinates should have the same length as values array. Sizes of fx: %d, f_vals: %d"
% (fx.size, f_vals.size)
)
matr = fx[:, None] ** (npix - np.arange(1, npix + 1))
# kf is the coeffs of the polynom passing trough the npix pixels
kf = np.linalg.inv(matr).dot(f_vals[:, None])
# compute coefficients of the 1st order derivative of the polynom
cf = ((npix - 1) - np.arange(npix - 1)) * kf[:-1, 0]
# revert for polyroots
cf = np.flip(cf)
# using Polynomial.fit, because supposed to be more numerically stable
# than previous solutions (according to numpy).
poly = Polynomial.fit(fx, f_vals, deg=2)
coeffs = poly.convert().coef
# maximum is the 0 of the 1st derivative of the polynom
zeros_pol = polyroots(cf)
# For a 1D parabola `f(x) = c + bx + ax^2`, the vertex position is:
# x_v = -b / 2a.
vertex_x = - coeffs[1] / (2 * coeffs[2])
# The pixel within the roots of the derivative of the polynom is
# the one near the approximate center of rotation (CoR)
approx_max = fx[np.argmax(f_vals)]
closest_pix_ind = np.argmin((zeros_pol - approx_max) ** 2)
return np.real(zeros_pol[closest_pix_ind])
vertex_min_x = np.min(fx)
vertex_max_x = np.max(fx)
if not (vertex_min_x < vertex_x < vertex_max_x):
raise ValueError(
"Fitted x: {} position is outide the margins of input: x: [{}, {}]".format(
vertex_x, vertex_min_x, vertex_max_x
)
)
return vertex_x
@staticmethod
def _determine_roi(img_shape, roi_yxhw, do_truncate_horz_pow2):
......
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