"""Elliptic Fourier Analysis"""
# Copyright 2020 Koji Noshita
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
from __future__ import annotations
import numpy as np
import numpy.typing as npt
import pandas as pd
from sklearn.base import (
BaseEstimator,
ClassNamePrefixFeaturesOutMixin,
TransformerMixin,
)
from sklearn.utils.parallel import Parallel, delayed
from ._registration import moment_register
# Tolerance for detecting degenerate (near-zero) geometric quantities
# (arc length, semi-axes, phase-angle denominator).
_DEGENERACY_TOL = 1e-15
# Floor value for near-zero arc-length segments when
# duplicated_points="infinitesimal".
_INFINITESIMAL_DT = 1e-10
# Tolerance for gimbal-lock detection in ZXZ Euler angle extraction.
_GIMBAL_TOL = 1e-10
[docs]
class EllipticFourierAnalysis(
ClassNamePrefixFeaturesOutMixin, TransformerMixin, BaseEstimator
):
r"""
Elliptic Fourier Analysis (EFA)
Parameters
----------
n_harmonics: int, default=20
Number of harmonics
n_dim: int, default=2
Dimension of the codomain (number of components on the
curve). Any positive integer is supported for the raw expansion
(``registration=None``).
Arc-length parameterization is computed automatically only for
``n_dim`` in ``(2, 3)`` (spatial shape coordinates). For ``n_dim=1``
or ``n_dim>3`` the codomain is treated as non-shape data, so the
parameter ``t`` is required and must be supplied explicitly (e.g.
derived from the corresponding shape); it is not inferred from the
codomain values.
registration : {"auto", None, "first_order", "moment"}, default="auto"
Shape-registration method. ``"auto"`` (default) registers 2D/3D shape
data with ``"first_order"`` and leaves other dimensions unregistered
(``None``). ``None`` returns raw coefficients. ``"first_order"`` uses
the 1st-harmonic ellipse; ``"moment"`` uses the inertia tensor of all
coefficients. Registration applies to 2D/3D shape data only.
`registration` must be ``None`` for``n_dim`` not in ``(2, 3)``.
scale : bool, default=True
Whether registration removes size. ``False`` keeps size.
Only used when ``registration != None``. (``scale=False``
is currently implemented for ``"moment"`` only.)
scale_method : str or None, default=None
Size measure when ``scale=True``; one of ``"semi_major_axis"``,
``"ellipse_area"``, ``"centroid_size"``, or ``None``.
``None`` resolves to the method's
default (``"first_order"``: 2D ``"semi_major_axis"``, 3D
``"ellipse_area"``; ``"moment"``: ``"centroid_size"``).
``"semi_major_axis"`` and ``"ellipse_area"`` require
``registration="first_order"``; ``"centroid_size"`` requires
``registration="moment"``.
align_parameter : bool, default=True
Whether to align the parameter domain (phase) in ``"first_order"``.
Only ``True`` is implemented (it is always applied); ``False`` raises
``NotImplementedError``.
reflect : bool, default=False
Whether to also remove reflection (chirality). Honored by
``"moment"``. For ``"first_order"`` only ``reflect=False`` (orientation
preserved) is implemented; ``reflect=True`` raises
``NotImplementedError``.
return_transform : bool, default=False
Append the registration parameters as extra output columns.
Currently only for ``"first_order"`` with ``n_dim`` in ``(2, 3)``;
other methods raise ``NotImplementedError``.
- 2D: Appends ``[psi, scale]`` to the end of the coefficient vector.
- 3D: Appends ``[alpha, beta, gamma, phi, scale]`` to the end.
n_jobs: int, default=None
The number of jobs to run in parallel. None means 1 unless in a
joblib.parallel_backend context. -1 means using all processors.
verbose: int, default=0
The verbosity level.
norm : bool, default=True
Deprecated alias. ``True`` maps to ``registration="first_order"``,
``False`` to ``registration=None``. Use ``registration`` instead.
norm_method : {None, "area", "semi_major_axis"}, default=None
Deprecated alias for ``scale_method`` (``"area"`` ->
``"ellipse_area"``). When ``None`` the dimension-appropriate
convention is used: ``"semi_major_axis"`` for 2D and ``"area"`` for 3D.
- ``"area"``: Scale by ``sqrt(pi * a1 * b1)`` where ``a1`` and ``b1``
are the semi-major and semi-minor axis lengths of the 1st harmonic
ellipse (Godefroy et al. 2012).
- ``"semi_major_axis"``: Scale by the semi-major axis length ``a1``
of the 1st harmonic ellipse (Kuhl & Giardina 1982).
return_orientation_scale : bool, default=False
Deprecated alias for ``return_transform`` producing identical output
(requires ``registration != None``).
Notes
-----
EFA is widely applied for outline shape analysis
in two-dimensional space [Kuhl_Giardina_1982]_.
.. math::
\begin{align}
x(l) &=
\frac{a_0}{2} + \sum_{i=1}^{n}
\left[ a_i \cos\left(\frac{2\pi i t}{T}\right)
+ b_i \sin\left(\frac{2\pi i t}{T}\right) \right]\\
y(l) &=
\frac{c_0}{2} + \sum_{i=1}^{n}
\left[ c_i \cos\left(\frac{2\pi i t}{T}\right)
+ d_i \sin\left(\frac{2\pi i t}{T}\right) \right]\\
\end{align}
EFA is also applied for a closed curve in the three-dimensional space
(e.g., [Lestrel_1997]_, [Lestrel_et_al_1997]_, and [Godefroy_et_al_2012]_).
For 3D data (``n_dim=3``), ``"first_order"`` registration follows
Godefroy et al. (2012) §3.1: rescaling by the 1st harmonic ellipse area,
reorientation using ZXZ Euler angles, phase shift, and direction correction.
When ``return_transform=True`` with registered data, extra values
are appended to the output:
- 2D: ``[psi, scale]`` where ``psi`` is the orientation angle and
``scale`` is the normalization factor.
- 3D: ``[alpha, beta, gamma, phi, scale]`` where ``(alpha, beta, gamma)``
are ZXZ Euler angles (in radians) of the 1st harmonic ellipse
orientation, ``phi`` is the phase angle, and ``scale`` is the
normalization factor.
The ``scale`` value depends on ``scale_method``: ``sqrt(pi * a1 * b1)``
for ``"ellipse_area"``, or ``a1`` for ``"semi_major_axis"``.
The legacy parameters ``norm`` / ``norm_method`` /
``return_orientation_scale`` remain as deprecated aliases of
``registration`` / ``scale_method`` / ``return_transform``.
References
----------
.. [Kuhl_Giardina_1982] Kuhl, F.P., Giardina, C.R. (1982) Elliptic Fourier
features of a closed contour. Comput. Graph. Image Process. 18: 236–258.
https://doi.org/10.1016/0146-664X(82)90034-X
.. [Lestrel_1997] Lestrel, P.E., 1997. Introduction and overview of Fourier
descriptors, in: Fourier Descriptors and Their Applications in Biology.
Cambridge University Press, pp. 22–44.
https://doi.org/10.1017/cbo9780511529870.003
.. [Lestrel_et_al_1997] Lestrel, P.E., Read, D.W., Wolfe, C., 1997. Size and
shape of the rabbit orbit: 3-D Fourier descriptors, in: Lestrel, P.E.
(Ed.), Fourier Descriptors and Their Applications in Biology. Cambridge
University Press, pp. 359–378. https://doi.org/10.1017/cbo9780511529870.017
.. [Godefroy_et_al_2012] Godefroy, J.E., Bornert, F., Gros, C.I.,
Constantinesco, A., 2012. Elliptical Fourier descriptors for contours in
three dimensions: A new tool for morphometrical analysis in biology.
C. R. Biol. 335, 205–213. https://doi.org/10.1016/j.crvi.2011.12.004
"""
_VALID_NORM_METHODS = {None, "area", "semi_major_axis"}
_VALID_REGISTRATIONS = {
None,
"first_order",
"moment",
"landmark",
"rotational_match",
}
_VALID_SCALE_METHODS = {
None,
"semi_major_axis",
"ellipse_area",
"centroid_size",
}
# Size measures permitted per registration method.
_SCALE_METHODS_BY_REGISTRATION = {
"first_order": {None, "semi_major_axis", "ellipse_area"},
"moment": {None, "centroid_size"},
}
# Deprecated norm_method -> new scale_method names.
_LEGACY_SCALE_METHOD_MAP = {
"area": "ellipse_area",
"semi_major_axis": "semi_major_axis",
}
# New scale_method -> legacy values understood by _normalize_2d/_normalize_3d.
_SCALE_METHOD_TO_LEGACY = {
"ellipse_area": "area",
"semi_major_axis": "semi_major_axis",
}
def __init__(
self,
n_harmonics: int = 20,
n_dim: int = 2,
registration: str | None = "auto",
scale: bool = True,
scale_method: str | None = None,
align_parameter: bool = True,
reflect: bool = False,
return_transform: bool = False,
n_jobs: int | None = None,
verbose: int = 0,
norm: bool = True,
norm_method: str | None = None,
return_orientation_scale: bool = False,
):
self.n_harmonics = n_harmonics
self.n_dim = n_dim
self.registration = registration
self.scale = scale
self.scale_method = scale_method
self.align_parameter = align_parameter
self.reflect = reflect
self.return_transform = return_transform
self.n_jobs = n_jobs
self.verbose = verbose
# Deprecated aliases (kept for backward compatibility).
self.norm = norm
self.norm_method = norm_method
self.return_orientation_scale = return_orientation_scale
def _resolve_registration(self):
"""Resolve effective registration settings from new + legacy params.
Legacy params (``norm`` / ``norm_method`` /
``return_orientation_scale``) map onto the new API. Setting a legacy
param to a non-default value together with its new counterpart raises
``ValueError``.
Returns
-------
tuple
``(method, scale, scale_method, return_extras)``.
"""
# registration <- norm (norm=False is the only active legacy signal;
# norm=True equals the default and does not override).
if self.norm is False:
if self.registration != "auto":
raise ValueError(
"Cannot set both `norm` (deprecated) and `registration`. "
"Use `registration` only."
)
method = None
scale = True
else:
method = self.registration
scale = self.scale
# Resolve "auto": register 2D/3D shape data, leave others unregistered.
if method == "auto":
method = "first_order" if self.n_dim in (2, 3) else None
# scale_method <- norm_method
if self.norm_method is not None:
if self.norm_method not in self._VALID_NORM_METHODS:
raise ValueError(
f"norm_method must be None, 'area', or 'semi_major_axis', "
f"got '{self.norm_method}'"
)
if self.scale_method is not None:
raise ValueError(
"Cannot set both `norm_method` (deprecated) and "
"`scale_method`. Use `scale_method` only."
)
scale_method = self._LEGACY_SCALE_METHOD_MAP[self.norm_method]
else:
scale_method = self.scale_method
# return columns <- return_transform / return_orientation_scale.
# Both flags request the same method-specific appended layout
# (2D `[psi, scale]`; 3D `[alpha, beta, gamma, phi, scale]`);
# `return_orientation_scale` is a deprecated alias of `return_transform`
# producing identical output. The registration transform is 2D/3D only,
# so there is no separate n-D / matrix "general" layout.
if self.return_orientation_scale and self.return_transform:
raise ValueError(
"Cannot set both `return_orientation_scale` (deprecated) "
"and `return_transform`. Use `return_transform` only."
)
return_extras = self.return_transform or self.return_orientation_scale
return method, scale, scale_method, return_extras
def _validate_registration(self, method, scale, scale_method, return_extras):
"""Validate effective registration settings; raise on bad combos."""
if method not in self._VALID_REGISTRATIONS:
valid = sorted(str(m) for m in self._VALID_REGISTRATIONS)
raise ValueError(f"registration must be one of {valid}, got '{method}'")
if method in ("landmark", "rotational_match"):
raise NotImplementedError(
f"registration='{method}' is reserved and not implemented yet."
)
if scale_method not in self._VALID_SCALE_METHODS:
raise ValueError(
f"scale_method must be one of "
f"{sorted(str(m) for m in self._VALID_SCALE_METHODS)}, "
f"got '{scale_method}'"
)
if method is None:
if return_extras:
raise ValueError(
"return_transform/return_orientation_scale requires "
"registration != None."
)
return
if self.n_dim not in (2, 3):
raise ValueError(
f"registration='{method}' applies to 2D/3D shape data only; "
f"got n_dim={self.n_dim}. Use registration=None for non-shape "
"data (normalization of n-D fields belongs to a separate "
"interface)."
)
allowed = self._SCALE_METHODS_BY_REGISTRATION[method]
if scale_method not in allowed:
raise ValueError(
f"scale_method='{scale_method}' is not valid for "
f"registration='{method}'; valid options: "
f"{sorted(str(m) for m in allowed)}."
)
if return_extras and not (method == "first_order" and self.n_dim in (2, 3)):
raise NotImplementedError(
"return_transform/return_orientation_scale is currently "
"implemented only for registration='first_order' with "
"n_dim in (2, 3)."
)
if method == "first_order" and not self.align_parameter:
raise NotImplementedError(
"align_parameter=False is not yet implemented; 'first_order' "
"always aligns the parameter domain (phase). Use "
"align_parameter=True."
)
if method == "first_order" and self.reflect:
raise NotImplementedError(
"reflect=True is not yet implemented for "
"registration='first_order' in EFA; orientation is preserved. "
"Use registration='moment' for reflection removal, or "
"reflect=False."
)
[docs]
def fit(self, X, y=None):
"""Fit the model (no-op for stateless transformer).
Parameters
----------
X : ignored
y : ignored
Returns
-------
self
"""
return self
###########################################################
#
# 2D
#
###########################################################
def _transform_single(
self,
X: np.ndarray,
t: np.ndarray | None = None,
duplicated_points: str = "infinitesimal",
):
"""Fit the model with a single outline.
Computes the raw Fourier coefficients for an ``n_dim``-valued closed
curve and, when ``self.norm`` is set, applies the dimension-specific
normalization.
Parameters
----------
X : ndarray of shape (n_coords, n_dim)
Coordinate values of the outline.
t : ndarray of shape (n_coords,), optional
A parameter indicating the position on the outline. For
``n_dim`` in ``(2, 3)`` (spatial shape coordinates), ``None``
triggers automatic arc-length parameterization. For ``n_dim=1``
or ``n_dim>3`` the codomain is treated as non-shape data, so
``t`` is required and is not inferred from the codomain values.
duplicated_points : str, default="infinitesimal"
Strategy for zero-length segments
(``"infinitesimal"`` or ``"deletion"``).
Returns
-------
X_transformed : ndarray
Flat Fourier coefficient vector. Per axis the layout is
``[cos_0..cos_n, sin_0..sin_n]``, concatenated over the axes.
"""
n_harmonics = self.n_harmonics
if t is None and self.n_dim not in (2, 3):
raise ValueError(
"Automatic arc-length parameterization is only available for "
f"n_dim in (2, 3); got n_dim={self.n_dim}. For n_dim=1 or "
"n_dim>3 the codomain is treated as non-shape data, so `t` "
"must be supplied explicitly (e.g. a parameterization derived "
"from the corresponding shape)."
)
X_arr, diffs, dt = _preprocess_outline(X, t, duplicated_points)
n_dim = X_arr.shape[1]
if n_dim != self.n_dim:
raise ValueError(
f"Each sample must have n_dim={self.n_dim} columns; got {n_dim}."
)
# Fourier series expansion
T = np.sum(dt)
# DC offset per axis (weighted mean over the contour)
offsets = 2.0 * (dt @ X_arr[1:]) / T
cos_rows = [
np.append(offsets[d], _cse(diffs[:, d], dt, n_harmonics))
for d in range(n_dim)
]
sin_rows = [
np.append(0.0, _sse(diffs[:, d], dt, n_harmonics)) for d in range(n_dim)
]
method, scale, scale_method, return_extras = self._resolve_registration()
self._validate_registration(method, scale, scale_method, return_extras)
def _raw_rows():
out = []
for d in range(n_dim):
out.append(cos_rows[d])
out.append(sin_rows[d])
return out
if method is None:
return np.hstack(_raw_rows())
if method == "moment":
if return_extras:
raise NotImplementedError(
"return_transform is not yet implemented for registration='moment'."
)
raw_flat = np.hstack(_raw_rows())
return moment_register(raw_flat, n_dim, scale=scale, reflect=self.reflect)
# method == "first_order"
if not scale:
raise NotImplementedError(
"scale=False (form space) is not yet implemented for "
"registration='first_order'."
)
# New scale_method names map onto the legacy values understood by
# _normalize_2d/_normalize_3d ("ellipse_area" -> "area").
legacy_sm = (
None if scale_method is None else self._SCALE_METHOD_TO_LEGACY[scale_method]
)
if n_dim == 2:
an, bn, cn, dn, psi, s = self._normalize_2d(
cos_rows[0],
sin_rows[0],
cos_rows[1],
sin_rows[1],
scale_method=legacy_sm,
)
rows = [an, bn, cn, dn]
extras = [psi, s]
elif n_dim == 3:
(
an,
bn,
cn,
dn,
en,
fn,
alpha,
beta,
gamma,
phi,
s,
) = self._normalize_3d(
cos_rows[0],
sin_rows[0],
cos_rows[1],
sin_rows[1],
cos_rows[2],
sin_rows[2],
scale_method=legacy_sm,
)
rows = [an, bn, cn, dn, en, fn]
extras = [alpha, beta, gamma, phi, s]
else: # unreachable: _validate_registration rejects non-2D/3D first.
raise ValueError(
f"registration applies to 2D/3D shape data only; got n_dim={n_dim}."
)
if return_extras:
return np.hstack(rows + extras)
return np.hstack(rows)
def _transform_single_2d(
self,
X: np.ndarray,
t: np.ndarray | None = None,
duplicated_points: str = "infinitesimal",
):
"""Backward-compatible 2D entry point.
Thin wrapper over :meth:`_transform_single`; retained for callers
that target the 2D path explicitly.
"""
return self._transform_single(X, t=t, duplicated_points=duplicated_points)
def _normalize_2d(self, an, bn, cn, dn, keep_start_point=False, scale_method=None):
"""Normalize Fourier coefficients.
``scale_method`` accepts the legacy values ``"area"`` /
``"semi_major_axis"``; ``None`` falls back to ``self.norm_method``.
Todo:
- [ ] Procrustes alignment -> in coordinate values?
Returns
-------
An, Bn, Cn, Dn : np.ndarray
Normalized coefficient arrays (offset + harmonics).
psi : float
Orientation (phase) of the 1st ellipse in radians.
scale : float
Scaling factor. semi-major axis length, or area of the 1st ellipse.
"""
a1 = an[1]
b1 = bn[1]
c1 = cn[1]
d1 = dn[1]
theta = 0.5 * np.arctan2(2 * (a1 * b1 + c1 * d1), a1**2 + c1**2 - b1**2 - d1**2)
[[a_s, b_s], [c_s, d_s]] = np.array([[a1, b1], [c1, d1]]).dot(
rotation_matrix_2d(theta)
)
s1 = a_s**2 + c_s**2
s2 = b_s**2 + d_s**2
if s1 < s2:
if theta < 0:
theta = theta + np.pi / 2
else:
theta = theta - np.pi / 2
cos_th = np.cos(theta)
sin_th = np.sin(theta)
a_s = a1 * cos_th + b1 * sin_th
c_s = c1 * cos_th + d1 * sin_th
semi_major = np.sqrt(a_s**2 + c_s**2)
norm_method = scale_method if scale_method is not None else self.norm_method
if norm_method is None:
norm_method = "semi_major_axis"
if norm_method == "semi_major_axis":
scale = semi_major
else: # "area"
b_s = -a1 * sin_th + b1 * cos_th
d_s = -c1 * sin_th + d1 * cos_th
semi_minor = np.sqrt(b_s**2 + d_s**2)
scale = np.sqrt(np.pi * semi_major * semi_minor)
psi = np.arctan2(c_s, a_s)
if keep_start_point:
theta = 0
coef_norm_list = []
r_psi = rotation_matrix_2d(-psi)
for n in range(1, len(an)):
r_ntheta = rotation_matrix_2d(n * theta)
coef_orig = np.array([[an[n], bn[n]], [cn[n], dn[n]]])
coef_norm_tmp = (1 / scale) * np.dot(np.dot(r_psi, coef_orig), r_ntheta)
coef_norm_list.append(coef_norm_tmp.reshape(-1))
coef_norm = np.stack(coef_norm_list)
An = np.append(an[0], coef_norm[:, 0])
Bn = np.append(bn[0], coef_norm[:, 1])
Cn = np.append(cn[0], coef_norm[:, 2])
Dn = np.append(dn[0], coef_norm[:, 3])
return An, Bn, Cn, Dn, psi, scale
def _inverse_transform_single(self, X_transformed, t_num=100):
coef_array = np.asarray(X_transformed, dtype=float)
n_axes = 2 * self.n_dim
n_extras = {2: 2, 3: 5}.get(self.n_dim, 0)
expected_base = n_axes * (self.n_harmonics + 1)
if coef_array.shape[0] == expected_base:
coef_core = coef_array
elif n_extras and coef_array.shape[0] == expected_base + n_extras:
coef_core = coef_array[:expected_base]
else:
allowed = (
f"{expected_base} or {expected_base + n_extras}"
if n_extras
else f"{expected_base}"
)
raise ValueError(
f"Expected {allowed} coefficients, got {coef_array.shape[0]}."
)
# Reshape to (n_axes, n_harmonics+1).
# Axes are ordered [cos0, sin0, cos1, sin1, ...] per coordinate.
axes = coef_core.reshape([n_axes, -1])
# Offsets sit at index 0 of the cos rows. Registered coefficients are
# translation-free (centered), so the stored offset is dropped.
offsets = axes[::2, 0].copy()
method, _, _, _ = self._resolve_registration()
if method is not None:
offsets[:] = 0.0
# (n_dim, n_harmonics)
cos_coefs = axes[::2, 1:]
sin_coefs = axes[1::2, 1:]
n_max = cos_coefs.shape[1]
theta = np.linspace(2 * np.pi / t_num, 2 * np.pi, t_num)
ns = np.arange(1, n_max + 1)
# (n_max, t_num)
cos_basis = np.cos(np.outer(ns, theta))
sin_basis = np.sin(np.outer(ns, theta))
# Reconstruct coordinates: (n_dim, t_num)
coords = offsets[:, None] / 2 + cos_coefs @ cos_basis + sin_coefs @ sin_basis
return coords.T
###########################################################
#
# 3D
#
###########################################################
def _transform_single_3d(
self,
X: np.ndarray,
t: np.ndarray | None = None,
duplicated_points: str = "infinitesimal",
):
"""Backward-compatible 3D entry point.
Thin wrapper over :meth:`_transform_single`; retained for callers
that target the 3D path explicitly.
"""
return self._transform_single(X, t=t, duplicated_points=duplicated_points)
def _normalize_3d(self, an, bn, cn, dn, en, fn, scale_method=None):
"""Normalize 3D EFA coefficients.
``scale_method`` accepts the legacy values ``"area"`` /
``"semi_major_axis"``; ``None`` falls back to ``self.norm_method``.
Applies the 4-step normalization algorithm:
1. Rescaling by a scale factor determined by ``scale_method``
(``None`` resolves to ``"area"`` for 3D):
- ``"area"``: ``scale = sqrt(pi * a1 * b1)``
- ``"semi_major_axis"``: ``scale = a1``
2. Reorientation using the 1st harmonic's Euler angles
3. Phase shift using the 1st harmonic's phase angle
4. Direction correction (sign of y-sine component)
Parameters
----------
an, bn, cn, dn, en, fn : np.ndarray of shape (n_harmonics+1,)
Raw Fourier coefficient arrays. Index 0 is the offset.
Returns
-------
An, Bn, Cn, Dn, En, Fn : np.ndarray of shape (n_harmonics+1,)
Normalized coefficient arrays.
alpha, beta, gamma : float
ZXZ Euler angles of the 1st harmonic ellipse.
phi : float
Phase angle of the 1st harmonic ellipse.
scale : float
Scaling factor. ``sqrt(pi * a1 * b1)`` when ``norm_method="area"``,
or ``a1`` when ``norm_method="semi_major_axis"``.
Notes
-----
When ``return_orientation_scale=True`` in 3D, these five values are
appended to the transform output in the order:
``[alpha, beta, gamma, phi, scale]``.
"""
# Extract geometric parameters of the 1st harmonic
phi1, a1, b1, alpha1, beta1, gamma1 = _compute_ellipse_geometry_3d(
an[1], bn[1], cn[1], dn[1], en[1], fn[1]
)
# Handle degenerate 1st harmonic
if a1 < _DEGENERACY_TOL:
raise ValueError(
"Degenerate 1st harmonic: the ellipse has near-zero semi-axes. "
"Cannot normalize 3D EFA coefficients."
)
# 1. Rescaling
norm_method = scale_method if scale_method is not None else self.norm_method
if norm_method is None:
norm_method = "area"
if norm_method == "semi_major_axis":
scale = a1
else:
# Area-based (Godefroy et al. 2012)
area1 = np.pi * a1 * b1
scale = np.sqrt(area1)
# 2. Reorientation matrix (Omega1_inv = Omega1^T)
Omega1 = rotation_matrix_3d_euler_zxz(alpha1, beta1, gamma1)
Omega1_inv = Omega1.T
n_harmonics = len(an) - 1
An = np.empty_like(an)
Bn = np.empty_like(bn)
Cn = np.empty_like(cn)
Dn = np.empty_like(dn)
En = np.empty_like(en)
Fn = np.empty_like(fn)
An[0] = an[0]
Bn[0] = bn[0]
Cn[0] = cn[0]
Dn[0] = dn[0]
En[0] = en[0]
Fn[0] = fn[0]
for k in range(1, n_harmonics + 1):
# Build 3x2 coefficient matrix
# C_k = [[an_k, bn_k], [cn_k, dn_k], [en_k, fn_k]]
C_k = np.array(
[
[an[k], bn[k]],
[cn[k], dn[k]],
[en[k], fn[k]],
]
)
# 3. Phase rotation uses k*phi1 for harmonic k
# Removing phase phi1 means substituting t -> t + phi1:
# new_xc = xc*cos(k*phi1) + xs*sin(k*phi1)
# new_xs = -xc*sin(k*phi1) + xs*cos(k*phi1)
# In matrix form: C_k @ R(-k*phi1) where R is the standard rotation matrix
angle_k = k * phi1
cos_k = np.cos(angle_k)
sin_k = np.sin(angle_k)
R_phase_k = np.array(
[
[cos_k, sin_k],
[-sin_k, cos_k],
]
)
# Apply: C'_k = (1/scale) * Omega1_inv @ C_k @ R_phase_k
C_norm = (1.0 / scale) * Omega1_inv @ C_k @ R_phase_k
An[k] = C_norm[0, 0]
Bn[k] = C_norm[0, 1]
Cn[k] = C_norm[1, 0]
Dn[k] = C_norm[1, 1]
En[k] = C_norm[2, 0]
Fn[k] = C_norm[2, 1]
# 4. Direction correction
# If the y-sine coefficient of the 1st harmonic is negative,
# negate all sine columns
if Dn[1] < 0:
Bn = -Bn
Dn = -Dn
Fn = -Fn
return An, Bn, Cn, Dn, En, Fn, alpha1, beta1, gamma1, phi1, scale
###########################################################
#
# set_output API
#
###########################################################
def __sklearn_is_fitted__(self):
"""Return True since this is a stateless transformer."""
return True
[docs]
def get_feature_names_out(
self, input_features: None | npt.ArrayLike = None
) -> np.ndarray:
"""Get output feature names.
Parameters
----------
input_features : ignored
Returns
-------
feature_names_out : ndarray of str objects
Transformed feature names.
"""
method, _, _, return_extras = self._resolve_registration()
include_orientation = (
return_extras and method == "first_order" and self.n_dim in (2, 3)
)
return np.asarray(self._build_feature_names(include_orientation), dtype=str)
@property
def _n_features_out(self):
"""Number of transformed output features."""
method, _, _, return_extras = self._resolve_registration()
base = (self.n_harmonics + 1) * (2 * self.n_dim)
if return_extras and method == "first_order" and self.n_dim in (2, 3):
if self.n_dim == 3:
return base + 5
if self.n_dim == 2:
return base + 2
return base
def _build_feature_names(self, include_orientation: bool) -> list[str]:
n = self.n_harmonics + 1
# Legacy letter names for the common 1D/2D/3D cases; systematic
# ``x{d}_cos_i`` / ``x{d}_sin_i`` names for higher dimensions.
legacy_letters = ["a", "b", "c", "d", "e", "f"]
feature_names: list[str] = []
if self.n_dim <= 3:
for letter in legacy_letters[: 2 * self.n_dim]:
feature_names += [f"{letter}_{i}" for i in range(n)]
else:
for d in range(self.n_dim):
feature_names += [f"x{d}_cos_{i}" for i in range(n)]
feature_names += [f"x{d}_sin_{i}" for i in range(n)]
if include_orientation:
if self.n_dim == 3:
feature_names += ["alpha", "beta", "gamma", "phi", "scale"]
elif self.n_dim == 2:
feature_names += ["psi", "scale"]
return feature_names
###########################################################
#
# utility functions
#
###########################################################
def _coord_names(n_dim: int) -> list[str]:
"""Return coordinate column names for ``n_dim``-valued reconstructions.
Uses ``x``/``y``/``z`` for ``n_dim <= 3`` and systematic ``x0``,
``x1``, ... names otherwise.
"""
base = ["x", "y", "z"]
if n_dim <= len(base):
return base[:n_dim]
return [f"x{d}" for d in range(n_dim)]
def _preprocess_outline(X, t, duplicated_points="infinitesimal"):
"""Prepare an outline for EFA.
Wraps the contour (prepends last point), computes coordinate differences
and arc-length parameterization, validates inputs, and handles duplicated
(zero-length) segments.
Parameters
----------
X : ndarray of shape (n_coords, n_dim)
Coordinate values of an outline.
t : ndarray of shape (n_coords,) or None
Positional parameter.
If None, arc-length parameterization is used.
duplicated_points : str
Strategy for zero-length segments:
``"infinitesimal"`` (default) or ``"deletion"``.
Returns
-------
X_arr : ndarray of shape (n_coords + 1, n_dim)
Wrapped coordinate array (last point prepended).
diffs : ndarray of shape (m, n_dim)
Per-axis coordinate differences (m <= n_coords after deletion).
dt : ndarray of shape (m,)
Parameter increments.
"""
X_arr = np.vstack([X[-1:], np.asarray(X)])
if not np.all(np.isfinite(X_arr)):
raise ValueError("Input coordinates must not contain NaN or Inf values.")
diffs = X_arr[1:] - X_arr[:-1]
if t is None:
dt = np.linalg.norm(diffs, axis=1)
else:
t_ = np.append(0, t)
dt = t_[1:] - t_[:-1]
tp = np.cumsum(dt)
if len(tp) != len(X):
raise ValueError(
"len(t) must have a same len(X), len(t): "
+ str(len(tp))
+ ", len(X): "
+ str(len(X))
)
if tp[-1] < _DEGENERACY_TOL:
raise ValueError(
"Degenerate outline: total arc length is near zero. "
"All points may be identical."
)
if duplicated_points == "infinitesimal":
dt[dt < _INFINITESIMAL_DT] = _INFINITESIMAL_DT
elif duplicated_points == "deletion":
idx_duplicated_points = np.where(dt == 0)[0]
if len(idx_duplicated_points) > 0:
diffs = np.delete(diffs, idx_duplicated_points, axis=0)
dt = np.delete(dt, idx_duplicated_points)
X_arr = np.delete(X_arr, idx_duplicated_points, 0)
else:
raise ValueError("'duplicated_points' must be 'infinitesimal' or 'deletion'")
return X_arr, diffs, dt
[docs]
def rotation_matrix_2d(theta):
rot_mat = np.array(
[[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]
)
return rot_mat
def rotation_matrix_3d_euler_zxz(alpha: float, beta: float, gamma: float) -> np.ndarray:
"""Construct 3x3 rotation matrix from ZXZ Euler angles.
The rotation is composed as Omega = R_gamma @ R_beta @ R_alpha,
following the convention in Godefroy et al. (2012) Fig. 1.
Parameters
----------
alpha, beta, gamma : float
ZXZ Euler angles in radians.
Returns
-------
rotation_matrix : np.ndarray of shape (3, 3)
Orthogonal rotation matrix with determinant +1.
"""
ca, sa = np.cos(alpha), np.sin(alpha)
cb, sb = np.cos(beta), np.sin(beta)
cg, sg = np.cos(gamma), np.sin(gamma)
return np.array(
[
[ca * cg - sa * cb * sg, -ca * sg - sa * cb * cg, sa * sb],
[sa * cg + ca * cb * sg, -sa * sg + ca * cb * cg, -ca * sb],
[sb * sg, sb * cg, cb],
]
)
def _compute_ellipse_geometry_3d(
xc: float,
xs: float,
yc: float,
ys: float,
zc: float,
zs: float,
) -> tuple[float, float, float, float, float, float]:
"""Compute geometric parameters of a 3D ellipse from Fourier coefficients.
Returns the canonical solution with the geometric major axis as the
reference: a >= b > 0 (``a`` is the semi-major axis), beta in [0, pi].
The phase ``phi`` is the major-axis phase and is not constrained to
]-pi/4, pi/4[ (constraining it would re-swap a/b).
Parameters
----------
xc, xs, yc, ys, zc, zs : float
Cosine and sine Fourier coefficients for x, y, z coordinates.
Returns
-------
phi : float
Phase angle of the major axis.
a : float
Semi-major axis length (a > 0).
b : float
Semi-minor axis length (b > 0).
alpha : float
First Euler angle (ZXZ convention).
beta : float
Second Euler angle, beta in [0, pi].
gamma : float
Third Euler angle (ZXZ convention).
"""
sum_c2 = xc**2 + yc**2 + zc**2
sum_s2 = xs**2 + ys**2 + zs**2
dot_cs = xc * xs + yc * ys + zc * zs
# phase angle phi
# 2*dot_cs / denom = -tan(2*phi)
# -> phi_0 = -(1/2) * arctan2(2*dot_cs, denom)
# For a perfect circle both arguments are zero; arctan2(0, 0) = 0.0
denom = sum_c2 - sum_s2
phi_0 = -0.5 * np.arctan2(2 * dot_cs, denom)
cos_phi = np.cos(phi_0)
sin_phi = np.sin(phi_0)
sin_2phi = np.sin(2 * phi_0)
a2 = sum_c2 * cos_phi**2 + sum_s2 * sin_phi**2 - dot_cs * sin_2phi
b2 = sum_c2 * sin_phi**2 + sum_s2 * cos_phi**2 + dot_cs * sin_2phi
# Enforce a >= b; if not, shift phi by pi/2
if a2 < b2:
a2, b2 = b2, a2
phi_0 = phi_0 + np.pi / 2 if phi_0 < 0 else phi_0 - np.pi / 2
# Use the major-axis (a >= b) phase directly. Do NOT renormalize phi into
# ]-pi/4, pi/4[: that re-swaps a/b and ties the orientation to the phase
# branch, flipping the registered shape by pi/2 with the start point.
phi = phi_0
a = np.sqrt(max(a2, 0.0))
b = np.sqrt(max(b2, 0.0))
if a < _DEGENERACY_TOL:
return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
cos_phi = np.cos(phi)
sin_phi = np.sin(phi)
# Recover rotation matrix Omega columns from coefficient vectors.
# The parametric equation relates coefficients to Omega and local-frame values:
# [xc, yc, zc]^T = Omega @ [a*cos(phi), b*sin(phi), 0]^T
# [xs, ys, zs]^T = Omega @ [-a*sin(phi), b*cos(phi), 0]^T
# Solving for the first two columns of Omega:
# col0 = (cos(phi)*[xc,yc,zc] - sin(phi)*[xs,ys,zs]) / a
# col1 = (sin(phi)*[xc,yc,zc] + cos(phi)*[xs,ys,zs]) / b
coef_c = np.array([xc, yc, zc])
coef_s = np.array([xs, ys, zs])
col0 = (cos_phi * coef_c - sin_phi * coef_s) / a
if b < _DEGENERACY_TOL:
col1 = np.zeros(3)
col2 = np.zeros(3)
else:
col1 = (sin_phi * coef_c + cos_phi * coef_s) / b
col2 = np.cross(col0, col1)
# Rotation matrix entries
Omega_11, Omega_21, Omega_31 = col0
Omega_12, Omega_22, Omega_32 = col1
Omega_13, Omega_23, Omega_33 = col2
# Extract Euler angles (ZXZ) from the rotation matrix
cos_beta = np.clip(Omega_33, -1.0, 1.0)
beta = float(np.arccos(cos_beta))
if abs(np.sin(beta)) < _GIMBAL_TOL:
# Gimbal lock
gamma = 0.0
alpha = float(np.arctan2(Omega_21, Omega_11))
else:
sin_beta = np.sin(beta)
gamma = float(np.arctan2(Omega_31 / sin_beta, Omega_32 / sin_beta))
alpha = float(np.arctan2(Omega_13 / sin_beta, -Omega_23 / sin_beta))
return float(phi), float(a), float(b), alpha, beta, float(gamma)
def _cse(dx: np.ndarray, dt: np.ndarray, n_harmonics: int) -> np.ndarray:
"""Cos series expansion n>=1
Parameters
----------
dx : np.ndarray
differences of coordinates
dt : np.ndarray
differences of parameter
n_harmonics : int
number of harmonics
Returns
-------
coef : np.ndarray
coefficients of cos series expansion
"""
t = np.concatenate([[0], np.cumsum(dt)])
T = t[-1]
cn = [
(T / (2 * (np.pi**2) * (n**2)))
* np.sum(
(dx / dt)
* (np.cos(2 * np.pi * n * t[1:] / T) - np.cos(2 * np.pi * n * t[:-1] / T))
)
for n in range(1, n_harmonics + 1, 1)
]
coef = np.array(cn)
return coef
def _sse(dx: np.ndarray, dt: np.ndarray, n_harmonics: int) -> np.ndarray:
"""Sin series expansion n>=1"""
t = np.concatenate([[0], np.cumsum(dt)])
T = t[-1]
cn = [
(T / (2 * (np.pi**2) * (n**2)))
* np.sum(
(dx / dt)
* (np.sin(2 * np.pi * n * t[1:] / T) - np.sin(2 * np.pi * n * t[:-1] / T))
)
for n in range(1, n_harmonics + 1, 1)
]
coef = np.array(cn)
return coef
