numerical_analysis 0.2.0

//! Differential geometry on continuous and discrete curves and surfaces.
//!
//! Everything here is generic over [`linear_isomorphic`] vector spaces, so the
//! same code works for any conforming vector type.
//!
//! - Derivatives of parametric curves: [`numerical_derivative`], [`tangent`].
//! - Gradients of scalar fields: [`numerical_gradient`].
//! - Curvature of implicit surfaces (level sets of a scalar field):
//!   [`implicit_curvature_features`] and the [`implicit_gaussian_curvature`] /
//!   [`implicit_mean_curvature`] shorthands. These work on the surface's second
//!   fundamental form, evaluated with Lagrange differentiation from
//!   [`crate::lagrange_polynomials`].
//! - Resampling polylines by arclength: [`sample_discrete_curve`],
//!   [`sample_discrete_bitangents`].

use euclidean::orthogonal_vector;
use linear_isomorphic::*;
use num_traits::float::TotalOrder;

use crate::lagrange_polynomials::partial_derivatives;

pub fn numerical_derivative<T, F, S>(x: S, f: &F, h: S) -> T
where
    T: InnerSpace<S>,
    S: linear_isomorphic::RealField,
    F: Fn(S) -> T,
{
    // Consider computing epsilons as:
    //https://scicomp.stackexchange.com/questions/42489/computing-numerical-derivatives/42491#42491

    // Convoluted but rustc fails to understand 2.0 * h.
    let denom = h + h;
    let num = f(x + h) - f(x - h);

    // Convoluted but rustc fails to understand num / denom.
    num * (S::from(1.0).unwrap() / denom)
}

/// Compute the tangent vector for some vector function (scalar functions are
/// vector functions) at parameter x. Result is *not* normalized.
pub fn tangent<T, F, S>(x: S, f: &F, h: S) -> T
where
    T: InnerSpace<S>,
    S: linear_isomorphic::RealField,
    F: Fn(S) -> T,
{
    numerical_derivative(x, f, h)
}

pub fn numerical_gradient<T, F, S>(point: &T, f: F, dim: usize, epsilon: S) -> T
where
    T: InnerSpace<S>,
    S: linear_isomorphic::RealField + core::ops::Mul<T, Output = T>,
    F: Fn(&T) -> S,
{
    let mut res = point.clone();
    let mut sample = point.clone();
    for i in 0..dim
    {
        sample[i] += epsilon;
        let v1 = f(&sample);
        sample[i] -= S::from(2.0).unwrap() * epsilon;
        let v2 = f(&sample);
        sample[i] = point[i];

        res[i] = (v1 - v2) / S::from(2.0).unwrap() * epsilon;
    }

    res
}

/// Returns the minimum curvature, maximum curvature and a local reference frame
/// with [direction of minimum curvature, direction of maximum curvature,
/// gradient].
pub fn implicit_curvature_features<V, S, F>(
    point: &V,
    fun: F,
    epsilon: S,
) -> (S, S, [V; 3])
where
    V: InnerSpace<S>,
    S: RealField
        + std::iter::Sum
        + std::iter::Product
        + TotalOrder
        + core::ops::Mul<V, Output = V>,
    F: Fn(&V) -> S + Clone,
{
    let n: V = numerical_gradient::<V, _, S>(point, &fun, 3, epsilon).normalized();
    let u: V = orthogonal_vector::<S, V>(&n);
    let v = n.cross(&u);

    let d1 = u.clone();
    let d2 = v.clone();
    let d3 = n.clone();

    // Express the implicit function in the point's tangent plane coordinate frame.
    let fun_as_tangent = move |x, y, z| {
        let mut transformed_point = V::default();
        transformed_point += u.clone() * x;
        transformed_point += v.clone() * y;
        transformed_point += n.clone() * z;

        transformed_point += point.clone();

        fun(&transformed_point)
    };

    let phi_uu = partial_derivatives(
        [S::from(0.).unwrap(); 3],
        &fun_as_tangent,
        [2, 0, 0],
        epsilon,
        3,
    );
    let phi_vv = partial_derivatives(
        [S::from(0.).unwrap(); 3],
        &fun_as_tangent,
        [0, 2, 0],
        epsilon,
        3,
    );
    let phi_uv = partial_derivatives(
        [S::from(0.).unwrap(); 3],
        &fun_as_tangent,
        [1, 1, 0],
        epsilon,
        3,
    );
    let phi_n = partial_derivatives(
        [S::from(0.).unwrap(); 3],
        &fun_as_tangent,
        [0, 0, 1],
        epsilon,
        3,
    );

    let gaussian_curvature = (phi_uu * phi_vv - (phi_uv * phi_uv)) / (phi_n * phi_n);
    let mean_curvature = (phi_uu + phi_vv) / (S::from(2.).unwrap() * phi_n.abs());

    let coeff = (mean_curvature * mean_curvature - gaussian_curvature)
        .abs()
        .sqrt();
    let kmin = mean_curvature - coeff;
    let kmax = mean_curvature + coeff;

    let t1 = d1.clone() * phi_uv + d2.clone() * (kmin * phi_n - phi_uu);
    let t2 = d1.clone() * (kmax * phi_n - phi_vv) + d2.clone() * phi_uv;

    (kmin, kmax, [t1, t2, d3])
}

/// Gaussian curvature of the level set of `fun` that passes through `point`.
///
/// This is the product of the two principal curvatures, `k1 * k2`, so a sphere
/// of radius `r` gives `1 / r^2`, a cylinder or a plane gives `0` because one
/// principal curvature vanishes, and the saddle `z = x^2 - y^2` gives `-4` at
/// the origin. The sign classifies the point: positive is elliptic, a local dome
/// or bowl; zero is parabolic or flat; negative is hyperbolic, a saddle.
///
/// `point` does not have to sit on the zero set of `fun`; whichever level set
/// runs through it is the one being measured. Unlike
/// [`implicit_mean_curvature`], the result does not depend on which way the
/// field increases — it divides by the squared normal derivative, so negating
/// `fun` leaves it unchanged. At a critical point of `fun` the gradient
/// vanishes, there is no tangent plane, and the result is `NaN`.
///
/// `epsilon` is the step for every finite difference taken: the gradient that
/// fixes the tangent frame, and the second derivatives of `fun` within it.
///
/// # Choosing `epsilon`
///
/// Smaller is not better, and this is less forgiving than
/// [`implicit_mean_curvature`] because the answer is a difference of products of
/// second derivatives, so each one's round-off enters twice. The second
/// derivatives divide by `epsilon^2`, and what that amplifies is the
/// representation error of the *value* of `fun` near `point`, not of its
/// curvature. In `f32`, with `fun` of order 10 near `point`, `epsilon` from
/// `0.02` to `0.2` holds five or more digits, while `0.01` returns `0.0549`
/// where the answer is `0.0625`. `f64` has digits to spare and holds up far
/// below that.
///
/// # Example
///
/// ```
/// use nalgebra::Vector3;
/// use numerical_analysis::differential_geometry::implicit_gaussian_curvature;
///
/// // The level set through (4, 0, 0) is the sphere of radius 4.
/// let field = |p: &Vector3<f64>| p.dot(p) - 4.0;
/// let curvature =
///     implicit_gaussian_curvature(&Vector3::new(4.0, 0.0, 0.0), field, 0.05);
///
/// assert!((curvature - 1.0 / 16.0).abs() < 1e-9);
/// ```
///
/// See [`implicit_mean_curvature`] for `(k1 + k2) / 2`, and
/// [`implicit_curvature_features`] for the principal curvatures themselves
/// together with the directions they act along.
pub fn implicit_gaussian_curvature<V, S, F>(point: &V, fun: F, epsilon: S) -> S
where
    V: InnerSpace<S>,
    S: RealField
        + std::iter::Sum
        + std::iter::Product
        + TotalOrder
        + core::ops::Mul<V, Output = V>,
    F: Fn(&V) -> S + Clone,
{
    let n = numerical_gradient(point, &fun, 3, epsilon).normalized();
    let u = orthogonal_vector(&n);
    let v = n.cross(&u);

    // Express the implicit function in the point's tangent plane coordinate frame.
    let fun_as_tangent = move |x, y, z| {
        let mut transformed_point = V::default();
        transformed_point += u.clone() * x;
        transformed_point += v.clone() * y;
        transformed_point += n.clone() * z;

        transformed_point += point.clone();

        fun(&transformed_point)
    };

    let phi_uu = partial_derivatives(
        [S::from(0.).unwrap(); 3],
        &fun_as_tangent,
        [2, 0, 0],
        epsilon,
        3,
    );
    let phi_vv = partial_derivatives(
        [S::from(0.).unwrap(); 3],
        &fun_as_tangent,
        [0, 2, 0],
        epsilon,
        3,
    );
    let phi_uv = partial_derivatives(
        [S::from(0.).unwrap(); 3],
        &fun_as_tangent,
        [1, 1, 0],
        epsilon,
        3,
    );
    let phi_n = partial_derivatives(
        [S::from(0.).unwrap(); 3],
        &fun_as_tangent,
        [0, 0, 1],
        epsilon,
        3,
    );

    (phi_uu * phi_vv - (phi_uv * phi_uv)) / (phi_n * phi_n)
}

/// Mean curvature of the level set of `fun` that passes through `point`.
///
/// Mean here is the average of the two principal curvatures, `(k1 + k2) / 2`,
/// so a sphere of radius `r` gives `1 / r`, a cylinder of radius `r` gives
/// `1 / (2 r)`, and a plane or a minimal surface gives `0`.
///
/// `point` does not have to sit on the zero set of `fun`; whichever level set
/// runs through it is the one being measured. The surface normal is taken to be
/// the gradient of `fun`, so the sign is tied to the direction the field
/// increases in, and negating `fun` negates the result. At a critical point of
/// `fun` the gradient vanishes, there is no tangent plane, and the result is
/// `NaN`.
///
/// `epsilon` is the step for every finite difference taken: the gradient that
/// fixes the tangent frame, and the second derivatives of `fun` within it.
///
/// # Choosing `epsilon`
///
/// Smaller is not better. The second derivatives divide by `epsilon^2`, so
/// round-off in `fun` is amplified by `1 / epsilon^2`, and what gets amplified
/// is the representation error of the *value* of `fun` near `point`, not of its
/// curvature. In `f32`, with `fun` of order 10 near `point`, `epsilon = 0.01`
/// loses about five of the seven digits available and returns `0.234` where the
/// answer is `0.25`, while `0.05` is accurate to seven digits. Prefer `0.02` to
/// `0.1` there. `f64` has digits to spare and holds up far below that.
///
/// # Example
///
/// ```
/// use nalgebra::Vector3;
/// use numerical_analysis::differential_geometry::implicit_mean_curvature;
///
/// // The level set through (4, 0, 0) is the sphere of radius 4.
/// let field = |p: &Vector3<f64>| p.dot(p) - 4.0;
/// let curvature =
///     implicit_mean_curvature(&Vector3::new(4.0, 0.0, 0.0), field, 0.05);
///
/// assert!((curvature - 1.0 / 4.0).abs() < 1e-9);
/// ```
///
/// See [`implicit_gaussian_curvature`] for `k1 * k2` and
/// [`implicit_curvature_features`] for the principal curvatures together with
/// the directions they act along.
pub fn implicit_mean_curvature<V, S, F>(point: &V, fun: F, epsilon: S) -> S
where
    V: InnerSpace<S>,
    S: RealField
        + std::iter::Sum
        + std::iter::Product
        + TotalOrder
        + core::ops::Mul<V, Output = V>,
    F: Fn(&V) -> S + Clone,
{
    let n = numerical_gradient(point, &fun, 3, epsilon).normalized();
    let u = orthogonal_vector(&n);
    let v = n.cross(&u);

    // Express the implicit function in the point's tangent plane coordinate frame.
    let fun_as_tangent = move |x, y, z| {
        let mut transformed_point = V::default();
        transformed_point += u.clone() * x;
        transformed_point += v.clone() * y;
        transformed_point += n.clone() * z;

        transformed_point += point.clone();

        fun(&transformed_point)
    };

    let phi_uu = partial_derivatives(
        [S::from(0.).unwrap(); 3],
        &fun_as_tangent,
        [2, 0, 0],
        epsilon,
        3,
    );
    let phi_vv = partial_derivatives(
        [S::from(0.).unwrap(); 3],
        &fun_as_tangent,
        [0, 2, 0],
        epsilon,
        3,
    );
    let phi_n = partial_derivatives(
        [S::from(0.).unwrap(); 3],
        &fun_as_tangent,
        [0, 0, 1],
        epsilon,
        3,
    );

    (phi_uu + phi_vv) / (S::from(2.).unwrap() * phi_n.abs())
}

/// Computes the Laplace operator for a graph.
/// https://en.wikipedia.org/wiki/Discrete_Laplace_operator#Graph_Laplacians
pub fn sample_discrete_curve<V, S>(t: S, curve: &[V]) -> V
where
    V: linear_isomorphic::InnerSpace<S>,
    S: linear_isomorphic::RealField,
{
    debug_assert!(t.is_finite());
    t.clamp(S::from(0.).unwrap(), S::from(1.).unwrap());

    let t = t.clamp(S::from(0.).unwrap(), S::from(1.).unwrap());

    let arclength = (0..curve.len() - 1)
        .map(|i| (curve[i + 1].clone() - curve[i].clone()).norm())
        .fold(S::from(0.).unwrap(), |acc, x| acc + x);

    let target = t * arclength;

    let mut current = S::from(0.).unwrap();
    let mut i = 0;
    loop
    {
        let l = (curve[i].clone() - curve[i + 1].clone()).norm();

        if current + l >= target || i >= curve.len() - 1
        {
            break;
        }
        current += l;
        i += 1;
    }

    let v = target - current;
    let extent = (curve[i].clone() - curve[i + 1].clone()).norm();

    let t = v / extent;
    curve[i].clone() * (S::from(1.).unwrap() - t) + curve[i + 1].clone() * t
}

pub fn sample_discrete_bitangents<V, S>(t: S, curve: &[V], bitangents: &[V]) -> V
where
    V: linear_isomorphic::InnerSpace<S>,
    S: linear_isomorphic::RealField,
{
    debug_assert!(t.is_finite());
    t.clamp(S::from(0.).unwrap(), S::from(1.).unwrap());

    let t = t.clamp(S::from(0.).unwrap(), S::from(1.).unwrap());

    let arclength = (0..curve.len() - 1)
        .map(|i| (curve[i + 1].clone() - curve[i].clone()).norm())
        .fold(S::from(0.).unwrap(), |acc, x| acc + x);

    let target = t * arclength;

    let mut current = S::from(0.).unwrap();
    let mut i = 0;
    loop
    {
        let l = (curve[i].clone() - curve[i + 1].clone()).norm();

        if current + l >= target || i >= curve.len() - 1
        {
            break;
        }
        current += l;
        i += 1;
    }

    let v = target - current;
    let extent = (curve[i].clone() - curve[i + 1].clone()).norm();

    let t = v / extent;
    bitangents[i].clone() * (S::from(1.).unwrap() - t) + bitangents[i + 1].clone() * t
}

// +| Tests |+ =======================================================
#[cfg(test)]
mod tests
{
    use super::*;

    // `nalgebra` rather than the internal `algebra` crate, which this published
    // crate cannot depend on.
    type Vec3 = nalgebra::Vector3<f64>;

    #[test]
    fn test_implicit_curvature()
    {
        let sphere_sdf = |p: &Vec3| p.dot(&p) - 4.;

        let curvature =
            implicit_mean_curvature(&Vec3::new(4.0, 0., 0.), sphere_sdf, 0.01);
        assert!((curvature - 0.25).abs() < 0.01, "{curvature}");

        let sphere_sdf = |p: &Vec3| p.dot(&p) - 4.;
        let curvature =
            implicit_gaussian_curvature(&Vec3::new(4.0, 0., 0.), sphere_sdf, 0.1);
        assert!((curvature - 1. / 16.).abs() < 0.01, "{}", curvature);
    }
}