"""Raup's model."""
# Copyright 2026 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.
import numpy as np
import numpy.typing as npt
from sklearn.base import (
BaseEstimator,
ClassNamePrefixFeaturesOutMixin,
TransformerMixin,
)
from ._generating_curve import _pad_orientation, _surfaces_to_frame, whorl_theta_range
def _validate_raup_params(w_r: float, t_r: float, d_r: float, r0: float) -> None:
if not w_r > 1.0:
raise ValueError(f"w_r (whorl expansion rate) must be > 1, got {w_r}")
if not (-1.0 < d_r < 1.0):
raise ValueError(f"d_r must be in (-1, 1), got {d_r}")
if not r0 > 0.0:
raise ValueError(f"r0 must be > 0, got {r0}")
def _raup_discriminant(w_r: float, t_r: float, d_r: float) -> float:
r"""Trajectory discriminant :math:`\Lambda` of Raup's model.
.. math::
\Lambda = 4\pi^2 (1 + D_R)^2
+ (\ln W_R)^2 \bigl[(1 + D_R)^2 + 4 T_R^2\bigr].
:math:`\Lambda` is the radicand of the arc-length relation and
the denominator of the Raup-to-growing tube conversion (Noshita 2014).
The published ``4 T_R`` term is corrected to ``4 T_R**2``.
"""
log_w = np.log(w_r)
one_p_d = 1.0 + d_r
return 4.0 * np.pi**2 * one_p_d**2 + log_w**2 * (one_p_d**2 + 4.0 * t_r**2)
def _rotation_x(angle: float) -> npt.NDArray[np.float64]:
c, s = np.cos(angle), np.sin(angle)
return np.array([[1.0, 0.0, 0.0], [0.0, c, -s], [0.0, s, c]])
def _rotation_z(angle: float) -> npt.NDArray[np.float64]:
c, s = np.cos(angle), np.sin(angle)
return np.array([[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]])
def _raup_surface(
w_r: float,
t_r: float,
d_r: float,
delta_r: float = 0.0,
gamma_r: float = 0.0,
r0: float = 1.0,
theta_range: npt.ArrayLike | None = None,
phi_range: npt.ArrayLike | None = None,
aperture=None,
) -> npt.NDArray[np.float64]:
r"""Surface of Raup's model (see :func:`raup`).
Raup's model [Raup_1965]_ [Raup_1966]_, built in the global frame
following the definition in [Noshita_2014]_:
.. math::
p(\theta, \phi) = r_0\, w_r^{\theta / 2\pi}\, R_z(\theta) \cdot
\bigl( R_z(\gamma_r)\, R_x(\delta_r)\, c(\phi) + o \bigr),
with offset :math:`o = ((1+d_r)/(1-d_r),\ 0,\ 2 t_r/(1-d_r))` and coiling
axis :math:`z`.
References
----------
.. [Raup_1965] Raup, D.M., Michelson, A., 1965. Theoretical Morphology of
the Coiled Shell. Science 147, 1294–1295.
.. [Raup_1966] Raup, D.M., 1966. Geometric analysis of shell coiling: general
problems. Journal of Paleontology 40, 1178–1190.
.. [Noshita_2014] Noshita, K., 2014. Quantification and geometric analysis of
coiling patterns in gastropod shells based on 3D and 2D image data.
Journal of Theoretical Biology 363, 93–104.
"""
_validate_raup_params(w_r, t_r, d_r, r0)
if aperture is not None:
raise NotImplementedError("general aperture shapes are not supported yet")
if theta_range is None:
theta_range = whorl_theta_range(4.0)
if phi_range is None:
phi_range = np.linspace(0.0, 2.0 * np.pi, 90)
theta = np.asarray(theta_range, dtype=float)
phi = np.asarray(phi_range, dtype=float)
vx = 2.0 * d_r / (1.0 - d_r) + 1.0
vz = 2.0 * t_r * (d_r / (1.0 - d_r) + 1.0)
offset = np.array([vx, 0.0, vz])
circle = np.column_stack([np.cos(phi), np.zeros_like(phi), np.sin(phi)])
rotation = _rotation_z(gamma_r) @ _rotation_x(delta_r)
point_local = circle @ rotation.T + offset # (n_phi, 3)
scale = r0 * w_r ** (theta / (2.0 * np.pi)) # (n_theta,)
plx, ply, plz = point_local[:, 0], point_local[:, 1], point_local[:, 2]
cos_t = np.cos(theta)[:, None]
sin_t = np.sin(theta)[:, None]
x = scale[:, None] * (plx[None, :] * cos_t - ply[None, :] * sin_t)
y = scale[:, None] * (plx[None, :] * sin_t + ply[None, :] * cos_t)
z = scale[:, None] * plz[None, :]
return np.stack([x, y, z], axis=-1)
[docs]
def raup(
w_r: float,
t_r: float,
d_r: float,
delta_r: float = 0.0,
gamma_r: float = 0.0,
r0: float = 1.0,
theta_range: npt.ArrayLike | None = None,
phi_range: npt.ArrayLike | None = None,
aperture=None,
output: str = "surface",
) -> npt.NDArray[np.float64]:
r"""Generate a form from Raup's model.
Raup’s logarithmic shell coiling model [Raup_1965]_ [Raup_1966]_ describes a shell
by a trajectory of a generating curve that expands, rotates, and translates
along a fixed coiling axis.
Parameters
----------
w_r : float
Whorl expansion rate :math:`W_R` (> 1).
t_r : float
Translation rate :math:`T_R`. ``t_r = 0`` gives a planispiral.
d_r : float
Relative position of generating curve :math:`D_R`, in ``(-1, 1)``.
delta_r, gamma_r : float, default = 0.0
Aperture orientation :math:`(\Delta, \Gamma)` (the global rotation is
``Rz(gamma_r) Rx(delta_r)``). ``(0, 0)`` is the classical radial-axial
aperture plane.
r0 : float, default = 1.0
Initial tube radius (scale).
theta_range : array-like of shape (n_theta,), optional
Coiling-angle (radians). Defaults to four whorls.
phi_range : array-like of shape (n_phi,), optional
Aperture-angle samples. Defaults to ``np.linspace(0, 2*pi, 90)``.
aperture : None
Aperture shape.
output : {"surface"}, default = "surface"
Form representation to return.
Only ``"surface"`` is implemented; other representations are reserved.
Returns
-------
X : ndarray of shape (n_theta, n_phi, 3)
Surface coordinates.
References
----------
.. [Raup_1965] Raup, D.M., Michelson, A., 1965. Theoretical Morphology of
the Coiled Shell. Science 147, 1294–1295.
.. [Raup_1966] Raup, D.M., 1966. Geometric analysis of shell coiling: general
problems. Journal of Paleontology 40, 1178–1190.
"""
if output != "surface":
raise NotImplementedError(
f"output={output!r} is reserved; only 'surface' is implemented"
)
return _raup_surface(
w_r, t_r, d_r, delta_r, gamma_r, r0, theta_range, phi_range, aperture
)
[docs]
def l_r(theta, w_r, t_r, d_r, r0=1.0):
r"""Arc length of growth trajectory at coiling angle ``theta``.
Maps the coiling angle :math:`\theta` of the Raup's model to
the arc length :math:`l_R` of the reference-point trajectory ([Noshita_2014]_):
.. math::
l_R(\theta) = r_0\,(W_R^{\theta / 2\pi} - 1)\,
\frac{\sqrt{\Lambda}}{(1 - D_R)\,\ln W_R},
where :math:`\Lambda` is the trajectory discriminant. The relation is
closed-form for constant parameters; ``theta`` may be array-like.
Parameters
----------
theta : array-like
Coiling angle :math:`\theta` (radians).
w_r, t_r, d_r : float
Raup parameters (``w_r > 1``, ``-1 < d_r < 1``).
r0 : float, default = 1.0
Initial tube radius (arc length scales with ``r0``).
Returns
-------
l_r : float or ndarray
Trajectory arc length at ``theta``.
See Also
--------
theta_r : Inverse function (arc length to coiling angle).
References
----------
.. [Noshita_2014] Noshita, K., 2014. Quantification and geometric analysis of
coiling patterns in gastropod shells based on 3D and 2D image data.
Journal of Theoretical Biology 363, 93–104.
"""
_validate_raup_params(w_r, t_r, d_r, r0)
theta = np.asarray(theta, dtype=float)
log_w = np.log(w_r)
sqrt_lambda = np.sqrt(_raup_discriminant(w_r, t_r, d_r))
out = (
r0
* (w_r ** (theta / (2.0 * np.pi)) - 1.0)
* sqrt_lambda
/ ((1.0 - d_r) * log_w)
)
return float(out) if out.ndim == 0 else out
[docs]
def theta_r(l_r, w_r, t_r, d_r, r0=1.0):
r"""Coiling angle of growth trajectory at arc length ``l_r``.
Inverse of :func:`l_r` (analytic, since the arc length is affine in
:math:`W_R^{\theta/2\pi}`):
.. math::
\theta(l_R) = \frac{2\pi}{\ln W_R}\,
\ln\!\left(1 + \frac{l_R\,(1 - D_R)\,\ln W_R}{r_0\,\sqrt{\Lambda}}\right).
Parameters
----------
l_r : array-like
Trajectory arc length.
w_r, t_r, d_r : float
Raup parameters (``w_r > 1``, ``-1 < d_r < 1``).
r0 : float, default = 1.0
Initial tube radius.
Returns
-------
theta : float or ndarray
Coiling angle :math:`\theta` (radians).
See Also
--------
l_r : Inverse function (coiling angle to arc length).
References
----------
.. [Noshita_2014] Noshita, K., 2014. Quantification and geometric analysis of
coiling patterns in gastropod shells based on 3D and 2D image data.
Journal of Theoretical Biology 363, 93–104.
"""
_validate_raup_params(w_r, t_r, d_r, r0)
l_r = np.asarray(l_r, dtype=float)
log_w = np.log(w_r)
sqrt_lambda = np.sqrt(_raup_discriminant(w_r, t_r, d_r))
out = (2.0 * np.pi / log_w) * np.log1p(
l_r * (1.0 - d_r) * log_w / (r0 * sqrt_lambda)
)
return float(out) if out.ndim == 0 else out
[docs]
class RaupModel(ClassNamePrefixFeaturesOutMixin, TransformerMixin, BaseEstimator):
"""Raup's model.
Raup’s logarithmic shell coiling model [Raup_1965]_ [Raup_1966]_.
``inverse_transform`` is the generative map
``Phi: (w_r, t_r, d_r, delta_r, gamma_r) -> form``.
``transform`` (parameter estimation from measurement data) is not implemented
yet.
Parameters
----------
r0 : float, default = 1.0
Initial tube radius (scale) used for generation.
estimator : str, default = "ml_2d"
Fitting method used by ``transform`` (not yet implemented).
n_jobs : int, optional
Reserved for parallelism.
verbose : int, default = 0
Verbosity level.
References
----------
.. [Raup_1965] Raup, D.M., Michelson, A., 1965. Theoretical Morphology of
the Coiled Shell. Science 147, 1294–1295.
.. [Raup_1966] Raup, D.M., 1966. Geometric analysis of shell coiling: general
problems. Journal of Paleontology 40, 1178–1190.
"""
def __init__(
self,
r0: float = 1.0,
estimator: str = "ml_2d",
n_jobs: int | None = None,
verbose: int = 0,
):
self.r0 = r0
self.estimator = estimator
self.n_jobs = n_jobs
self.verbose = verbose
[docs]
def fit(self, X, y=None):
"""No-op (stateless). Returns self."""
return self
[docs]
def get_feature_names_out(self, input_features=None) -> np.ndarray:
"""Parameter names ``(w_r, t_r, d_r, delta_r, gamma_r)``."""
return np.asarray(["w_r", "t_r", "d_r", "delta_r", "gamma_r"], dtype=str)
