Differentials#
Overview#
The differentials module in SplineOps provides a collection of algorithms for the computation of image differentials based on cubic B-spline interpolation [1], [2]. By modeling a grayscale image as a continuous function reconstructed from its discrete samples, the module enables the accurate computation of derivatives. It offers several operations such as
Gradient Magnitude: the local rate of change of the intensity;
Gradient Direction: the direction along which the intensity changes most;
Laplacian:the sum of second-order derivatives;
Largest Hessian Eigenvalue: the maximal curvature;
Smallest Hessian Eigenvalue: the minimal curvature; and
Hessian Orientation: the principal orientation of the curvature.
The numerical behavior is deliberately separate from visualization:
Differentials.run returns raw values, does not replace the source image,
and performs no console output. normalize=True is an explicit convenience
for non-angular display maps. spacing gives one physical sample distance
per array axis; all values default to one.
Image Representation#
A grayscale image is modeled as the continuous function
where \(\varphi(x_1,x_2)\) is defined as the tensor-product cubic B‑spline
This formulation allows one to compute exact derivatives of the image by first determining the spline coefficients \(c[k_1,k_2]\).
Differentiation Operations#
Based on the spline representation, the module computes several differential operators. We define
Gradient Magnitude: Computed as the Euclidean norm of the first derivatives
\[\|\pmb{\nabla}f(x,y)\| = \sqrt{\bigl(f_1\bigr)^2 + \bigl(f_2\bigr)^2}.\]Gradient Direction: The direction of the gradient is given by
\[\theta(x,y) = \operatorname{atan2}\!\bigl(f_\mathrm{row}, f_\mathrm{column}\bigr).\]Laplacian: A second-order operator that highlights regions of rapid intensity change as
\[\Delta f(x,y) = f_{11} + f_{22}.\]Largest Hessian Eigenvalue: The maximal eigenvalue of the Hessian matrix is given by
\[\lambda_{\text{max}} = \tfrac12\Bigl(f_{11} + f_{22} + \sqrt{4f_{12}^2 + (f_{11} - f_{22})^2}\Bigr).\]Smallest Hessian Eigenvalue: The minimal eigenvalue of the Hessian matrix is given by
\[\lambda_{\text{min}} = \tfrac12\Bigl(f_{11} + f_{22} - \sqrt{4f_{12}^2 + (f_{11} - f_{22})^2}\Bigr).\]Hessian Orientation: This operation returns the orientation that corresponds to the maximal second derivative, as
\[\theta_H(x,y) = \pm \tfrac12 \arccos\!\Bigl(\tfrac{f_{11} - f_{22}}{\sqrt{4f_{12}^2 + (f_{11}-f_{22})^2}}\Bigr),\]where the sign is determined by the sign of the cross derivative \(f_{12}\).
Implementation Details#
The Differentials class implements these operations as follows: the input image is provided by its samples. We assume it to be a cubic B-spline and first determine its interpolation coefficients. The differential-based computations that we perform are then perfectly consistent with this continuously defined function. We finally build the output image by sampling the ideal, continuously defined intermediate result.
The preferred public class is Differentials; the historical lowercase
differentials alias remains available. Direct component methods avoid
reconstructing quantities from composite maps:
from splineops.differentials import Differentials
operator = Differentials(image, spacing=(0.7, 1.3))
vertical, horizontal = operator.gradient_components()
vertical2, cross, horizontal2 = operator.hessian_components()
vertical follows increasing row coordinates and horizontal follows
increasing column coordinates. Mirror-boundary derivatives are zero at the
outermost sample for the antisymmetric first-derivative filter. Polynomial
and trigonometric fields are used as analytical tests away from that boundary.
Volumes and multi-output plans#
The same component interface accepts one scalar 3-D volume. Components are
returned in increasing axis order, while Hessians use packed
upper-triangular order (00, 01, 02, 11, 12, 22). Gradient magnitude,
Laplacian, and minimum/maximum Hessian eigenvalue maps generalize to 3-D;
gradient direction and Hessian orientation remain explicitly 2-D quantities.
Use splineops.differentials.DifferentialPlan when several derivative
families are needed together, or when changing volumes share one shape and
spacing:
from splineops.differentials import DifferentialPlan
plan = DifferentialPlan(volume.shape, spacing=(0.7, 0.7, 1.5))
result = plan.apply(volume, gradient=True, hessian=True)
gx, gy, gz = result.gradient
hxx, hxy, hxz, hyy, hyz, hzz = result.hessian
laplacian = result.laplacian
Output families can be selected independently. laplacian=None preserves
the original behavior and follows hessian; request it explicitly to avoid
building gradients, packed mixed Hessians, or their output arrays:
laplacian = plan.apply(
volume,
gradient=False,
hessian=False,
laplacian=True,
).laplacian
For repeated allocation-sensitive calls, pass a
splineops.differentials.DifferentialResult containing exact-shape and
exact-dtype destination arrays through out=. Fields corresponding to
unrequested families must be None; the returned object is the supplied
buffer container. All destinations are checked before computation begins.
Writable strided arrays are accepted, while read-only arrays, overlap with the
source, and overlap between output components are rejected. This prevents
partial writes and preserves the promise that differentiation does not mutate
its input.
One batched per-call workspace is used, so shared coefficient and derivative
intermediates are not recomputed for the requested outputs.
DifferentialPlan also accepts explicit spatial_axes for batch and
channel arrays; components retain the full input shape and follow the selected
axis order. The legacy Differentials object continues to model one scalar
image or volume, keeping its historical component helpers uncomplicated.
See Stability-soak workflows for a batched 3-D affine-to-feature pipeline.
The following figure, taken from Differentials Module, shows some of these differential maps (gradient magnitude, gradient direction, Laplacian and Hessian-based quantities).
