Merge branch 'feature/docs'

.gitignore

@@ -5,6 +5,10 @@ target/
 # files
 .DS_Store
 **/*.rs.bk
+*.aux
+*.log
+*.gz
+*.pdf
 
 #/target
 /Cargo.lock

Cargo.toml

@@ -1,7 +1,7 @@
 [package]
 name = "bsplines"
 license = "Apache-2.0"
-version = "0.0.1-alpha.5"
+version = "0.0.1-alpha.6"
 authors = [
     "Michael A. Heuer <[email protected]>",
 ]

README.md

@@ -1,4 +1,4 @@
-# Bsplines
+# `bsplines` Rust Library
 
 [![Crates.io](https://img.shields.io/crates/v/bsplines)](https://crates.io/crates/bsplines)
 [![Docs.rs](https://docs.rs/bsplines/badge.svg)](https://docs.rs/bsplines)

doc-images/equations/_preamble.tex

@@ -0,0 +1,17 @@
+\documentclass[tikz]{standalone}
+\usepackage{amssymb}
+\usepackage{amsmath}
+\usepackage{oubraces}
+\usepackage[customcolors]{hf-tikz}
+\definecolor{docsrscolor}{HTML}{f5f5f5}
+
+\newcommand{\myeqs}[1]{
+    \Huge
+    \begin{tikzpicture}
+        \node[
+            fill=docsrscolor,
+            rounded corners=8mm,
+            inner sep=8mm
+        ]{$#1$};
+    \end{tikzpicture}
+}

doc-images/equations/basis-function-zero.tex

@@ -0,0 +1,11 @@
+\input{_preamble}
+\begin{document}
+    \myeqs{
+        \mathcal{N}_{i,0}^{\boldsymbol{U}^{(k)}}(u)
+        =
+        \begin{cases}
+            1, & u\in\left[u_i^{(k)},u_{i+1}^{(k)}\right)\,\lor\, (i=n-k \land u=1)\\
+            0, & \text{else}
+        \end{cases}
+    }
+\end{document}

doc-images/equations/basis-function.tex

@@ -0,0 +1,10 @@
+\input{_preamble}
+\begin{document}
+    \myeqs{
+        \myeqs{
+            \mathcal{N}_{i,p-k}^{\boldsymbol{U}^{(k)}}(u)
+            = \mu_{i,p-k-1}^{(k)}(u)\, \mathcal{N}_{i,p-k-1}^{\boldsymbol{U}^{(k)}}(u)
+            + \left[1-\mu_{i+1,p-k-1}^{(k)}(u)\right]\, \mathcal{N}_{i+1,p-k-1}^{\boldsymbol{U}^{(k)}}(u)
+        }
+    }
+\end{document}

doc-images/equations/basis-prefactor.tex

@@ -0,0 +1,11 @@
+\input{_preamble}
+\begin{document}
+    \myeqs{
+        \mu_{g,h}^{(k)}(u)
+        =
+        \begin{cases}
+            0 & \text{if}\quad u_{g+h+1}^{(k)} = u_{g}^{(k)}\\[2mm]
+            \frac{u-u_g^{(k)}}{u_{g+h+1}^{(k)}-u_g^{(k)}} & \text{else}
+        \end{cases}
+    }
+\end{document}

doc-images/equations/control-points.tex

@@ -0,0 +1,13 @@
+\input{_preamble}
+\begin{document}
+    \myeqs{
+        \boldsymbol{P}_i^{(k)}=
+        \begin{cases}
+            \boldsymbol{P}_{i}^{(0)} & k=0 \\[2mm]
+            \left.\begin{cases}
+                      \boldsymbol{0} & u_{i+p+1}^{{(0)}} = u_{i+k}^{{(0)}}\\[2mm]
+                      \frac{p-k+1}{u_{i+p+1}^{{(0)}} - u_{i+k}^{{(0)}}} \left( \boldsymbol{P}_{i+1}^{(k-1)} - \boldsymbol{P}_{i}^{(k-1)} \right) & \text{else}%u_{i+p+1} \neq u_{i+k}\\[2mm]
+            \end{cases}\right|& k>0 \\
+        \end{cases}
+    }
+\end{document}

doc-images/equations/curve.tex

@@ -1,13 +1,8 @@
-\documentclass[10pt]{article}
-\usepackage[usenames]{color}
-\usepackage{amssymb}
-\usepackage{amsmath}
-\usepackage{nicefrac}
-\definecolor{mygreen}{rgb}{0.454,0.824,0.208}
-\definecolor{myred}{rgb}{0.8,0.173,0.137}
-
-\usepackage[utf8]{inputenc}
-\begin{equation}\nonumber
-\mathcal{C}(u) = \sum_{i=0}^{n} \mathcal{N}_{i,p}^{\boldsymbol{U}} (u)\, \boldsymbol{P}_i,\quad u \in [u_p,u_{n+1}]
-\end{equation}
+\input{_preamble}
+\begin{document}
+    \myeqs{
+        \mathcal{C}^{(k)}(u) = \frac{\partial^k\mathcal{C}^{(0)}(u)}{\partial u^k} = \sum_{i=0}^{n-k} \mathcal{N}_{i,p-k}^{\boldsymbol{U^{(k)}}} (u)\, \boldsymbol{P}^{(k)}_i,\quad
+        u \in [u_{p-k},u_{n+1-k}],\quad
+        n > p
+    }
 \end{document}

doc-images/equations/generate-equations-images.sh

@@ -0,0 +1,10 @@
+for file in *.tex; do
+    if [ "$file" != "_preamble.tex" ]; then
+        echo "Processing file: $file"
+        pdflatex -shell-escape -synctex=1 "$file" | grep '^!.*' -A200
+        pdf2svg "${file%.tex}.pdf" "${file%.tex}.svg"
+        rm -f "${file%.tex}.aux" "${file%.tex}.log" "${file%.tex}.pdf" "${file%.tex}.pdf" "${file%.tex}.synctex.gz"
+    else
+        echo "Excluded file: $file"
+    fi
+done

doc-images/equations/knots.tex

@@ -0,0 +1,16 @@
+\input{_preamble}
+\begin{document}
+    \myeqs{
+        \boldsymbol{U}^{(k)}
+        =\left\{u_i^{(k)}\right\}
+        \equiv\overunderbraces{
+            &&\br{5}{\text{domain}\vphantom{p}}
+        }{
+            \left\{\vphantom{u^{(k)}}\right.
+            &u_{0}^{(k)},\dots,&u_{p-k}^{(k)}&,
+            &u_{p-k+1}^{(k)}, \dots,u_{n-k}^{(k)}
+            &,&u_{n+1-k}^{(k)}&,\dots,u_{n+p+1-2k}^{(k)}& \left.\vphantom{u^{(k)}}\right\}}
+        {&\br{2}{p-k+1}&&\br{1}{n-p}&&\br{2}{1+p-k}},\quad
+        u_{i}^{(k)}\leq u_{i+1}^{(k)},
+    }
+\end{document}

src/curve/basis/mod.rs

@@ -0,0 +1,191 @@
+#![cfg_attr(feature = "doc-images",
+cfg_attr(all(),
+doc = ::embed_doc_image::embed_image!("eq-basis-function", "doc-images/equations/basis-function.svg"),
+doc = ::embed_doc_image::embed_image!("eq-basis-prefactor", "doc-images/equations/basis-prefactor.svg"),
+doc = ::embed_doc_image::embed_image!("eq-basis-function-zero", "doc-images/equations/basis-function-zero.svg")))]
+//! Evaluates the basis spline functions using the Cox-de Boor-Mansfield recurrence relation
+//!
+//! ![The Cox-de Boor-Mansfield recurrence relation][eq-basis-function]
+//!
+//! with the basis functions of degree `p = 0`
+//!
+//! ![Basis function of degree zero][eq-basis-function-zero]
+//!
+//! where the conditional `⋁ (i = n - k ⋀ u = U_{n+1-k)` closes the last interval
+//! and the pre-factors
+//!
+//! ![Pre-factors][eq-basis-prefactor]
+
+use crate::types::VecD;
+
+/// Evaluates the `i`-th basis spline function of degree `p`
+///
+/// ## Arguments
+///
+/// - `i` the index with `i ∈ {0, 1, ..., n}`
+/// - `p` the spline degree
+/// - `k` the derivative order
+/// - `U` the knot vector
+pub fn basis(Uk: &VecD, i: usize, p: usize, k: usize, n: usize, u: f64) -> f64 {
+    if p == 0 {
+        if (Uk[i] <= u && u < Uk[i + 1]) || (i == n - k && u == Uk[n + 1 - k]) {
+            return 1.0;
+        }
+        return 0.0;
+    }
+
+    let summand1 = if Uk[i + p] == Uk[i] {
+        0.0
+    } else {
+        let g = i;
+        let h = p - 1;
+        (u - Uk[g]) / (Uk[g + h + 1] - Uk[g]) * basis(Uk, i, h, k, n, u)
+    };
+
+    let summand2 = if Uk[i + 1 + p] == Uk[i + 1] {
+        0.0
+    } else {
+        let g = i + 1;
+        let h = p - 1;
+
+        // The following equation is numerically more stable than
+        // `(1.0 - ((u - Uk[g]) / (Uk[g + h + 1] - Uk[g]))) * self.evaluate(k, g, h, u)`
+        (Uk[g + p] - u) / (Uk[g + h + 1] - Uk[g]) * basis(Uk, g, h, k, n, u)
+    };
+
+    summand1 + summand2
+}
+
+#[cfg(test)]
+mod tests {
+    use approx::assert_relative_eq;
+    use nalgebra::dvector;
+
+    use crate::curve::knots::Knots;
+
+    const SEGMENTS: usize = 4;
+
+    #[test]
+    fn basis_func_degree3() {
+        let k = 0;
+        let p = 3;
+        let knots = Knots::new(p, dvector![0., 0., 0., 0., 1. / 3., 2. / 3., 1., 1., 1., 1.]);
+
+        // Basis function i = 0
+        let mut i = 0;
+        assert_eq!(knots.evaluate(k, i, p, 0.0), 1.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 6.), 1. / 8.);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 3.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 2.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 2. / 3.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 5. / 6.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1.), 0.0);
+
+        i = 1;
+        assert_eq!(knots.evaluate(k, i, p, 0.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 6.), 19. / 32.);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 3.), 1. / 4.);
+        assert_relative_eq!(knots.evaluate(k, i, p, 1. / 2.), 1. / 32., epsilon = f64::EPSILON.sqrt());
+        assert_eq!(knots.evaluate(k, i, p, 2. / 3.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 5. / 6.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1.), 0.0);
+
+        i = 2;
+        assert_eq!(knots.evaluate(k, i, p, 0.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 6.), 25. / 96.);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 3.), 7. / 12.);
+        assert_relative_eq!(knots.evaluate(k, i, p, 1. / 2.), 15. / 32., epsilon = f64::EPSILON.sqrt());
+        assert_relative_eq!(knots.evaluate(k, i, p, 2. / 3.), 1. / 6., epsilon = f64::EPSILON.sqrt());
+        assert_relative_eq!(knots.evaluate(k, i, p, 5. / 6.), 1. / 48., epsilon = f64::EPSILON.sqrt());
+        assert_eq!(knots.evaluate(k, i, p, 1.0), 0.0);
+
+        i = 3;
+        assert_eq!(knots.evaluate(k, i, p, 0.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 6.), 1. / 48.);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 3.), 1. / 6.);
+        assert_relative_eq!(knots.evaluate(k, i, p, 1. / 2.), 15. / 32., epsilon = f64::EPSILON.sqrt());
+        assert_relative_eq!(knots.evaluate(k, i, p, 2. / 3.), 7. / 12., epsilon = f64::EPSILON.sqrt());
+        assert_relative_eq!(knots.evaluate(k, i, p, 5. / 6.), 25. / 96., epsilon = f64::EPSILON.sqrt());
+        assert_eq!(knots.evaluate(k, i, p, 1.0), 0.0);
+
+        i = 4;
+        assert_eq!(knots.evaluate(k, i, p, 0.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 6.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 3.), 0.0);
+        assert_relative_eq!(knots.evaluate(k, i, p, 1. / 2.), 1. / 32., epsilon = f64::EPSILON.sqrt());
+        assert_relative_eq!(knots.evaluate(k, i, p, 2. / 3.), 1. / 4., epsilon = f64::EPSILON.sqrt());
+        assert_relative_eq!(knots.evaluate(k, i, p, 5. / 6.), 19. / 32., epsilon = f64::EPSILON.sqrt());
+        assert_eq!(knots.evaluate(k, i, p, 1.0), 0.0);
+
+        i = 5;
+        assert_eq!(knots.evaluate(k, i, p, 0.0), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 6.), 0.);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 3.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 2.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 2. / 3.), 0.0);
+        assert_relative_eq!(knots.evaluate(k, i, p, 5. / 6.), 1. / 8., epsilon = f64::EPSILON.sqrt());
+        assert_eq!(knots.evaluate(k, i, p, 1.), 1.0);
+    }
+
+    #[test]
+    fn basis_func_degree4() {
+        let k = 1;
+        let p = 4;
+        let knots = Knots::new(p, dvector![0., 0., 0., 0., 0., 1. / 3., 2. / 3., 1., 1., 1., 1., 1.]);
+
+        // Basis function i = 0
+        let mut i = 0;
+        assert_eq!(knots.evaluate(k, i, p, 0.0), 1.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 6.), 1. / 8.);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 3.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 2.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 2. / 3.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 5. / 6.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1.), 0.0);
+
+        i = 1;
+        assert_eq!(knots.evaluate(k, i, p, 0.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 6.), 19. / 32.);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 3.), 1. / 4.);
+        assert_relative_eq!(knots.evaluate(k, i, p, 1. / 2.), 1. / 32., epsilon = f64::EPSILON.sqrt());
+        assert_eq!(knots.evaluate(k, i, p, 2. / 3.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 5. / 6.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1.), 0.0);
+
+        i = 2;
+        assert_eq!(knots.evaluate(k, i, p, 0.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 6.), 25. / 96.);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 3.), 7. / 12.);
+        assert_relative_eq!(knots.evaluate(k, i, p, 1. / 2.), 15. / 32., epsilon = f64::EPSILON.sqrt());
+        assert_relative_eq!(knots.evaluate(k, i, p, 2. / 3.), 1. / 6., epsilon = f64::EPSILON.sqrt());
+        assert_relative_eq!(knots.evaluate(k, i, p, 5. / 6.), 1. / 48., epsilon = f64::EPSILON.sqrt());
+        assert_eq!(knots.evaluate(k, i, p, 1.0), 0.0);
+
+        i = 3;
+        assert_eq!(knots.evaluate(k, i, p, 0.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 6.), 1. / 48.);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 3.), 1. / 6.);
+        assert_relative_eq!(knots.evaluate(k, i, p, 1. / 2.), 15. / 32., epsilon = f64::EPSILON.sqrt());
+        assert_relative_eq!(knots.evaluate(k, i, p, 2. / 3.), 7. / 12., epsilon = f64::EPSILON.sqrt());
+        assert_relative_eq!(knots.evaluate(k, i, p, 5. / 6.), 25. / 96., epsilon = f64::EPSILON.sqrt());
+        assert_eq!(knots.evaluate(k, i, p, 1.0), 0.0);
+
+        i = 4;
+        assert_eq!(knots.evaluate(k, i, p, 0.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 6.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 3.), 0.0);
+        assert_relative_eq!(knots.evaluate(k, i, p, 1. / 2.), 1. / 32., epsilon = f64::EPSILON.sqrt());
+        assert_relative_eq!(knots.evaluate(k, i, p, 2. / 3.), 1. / 4., epsilon = f64::EPSILON.sqrt());
+        assert_relative_eq!(knots.evaluate(k, i, p, 5. / 6.), 19. / 32., epsilon = f64::EPSILON.sqrt());
+        assert_eq!(knots.evaluate(k, i, p, 1.0), 0.0);
+
+        i = 5;
+        assert_eq!(knots.evaluate(k, i, p, 0.0), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 6.), 0.);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 3.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 1. / 2.), 0.0);
+        assert_eq!(knots.evaluate(k, i, p, 2. / 3.), 0.0);
+        assert_relative_eq!(knots.evaluate(k, i, p, 5. / 6.), 1. / 8., epsilon = f64::EPSILON.sqrt());
+        assert_eq!(knots.evaluate(k, i, p, 1.), 1.0);
+    }
+}