diff --git a/.gitignore b/.gitignore index b80da91..041d818 100644 --- a/.gitignore +++ b/.gitignore @@ -5,6 +5,10 @@ target/ # files .DS_Store **/*.rs.bk +*.aux +*.log +*.gz +*.pdf #/target /Cargo.lock diff --git a/Cargo.toml b/Cargo.toml index 9097d0c..b173722 100644 --- a/Cargo.toml +++ b/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 ", ] diff --git a/README.md b/README.md index 2bc3579..beece4e 100644 --- a/README.md +++ b/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) diff --git a/doc-images/equations/_preamble.tex b/doc-images/equations/_preamble.tex new file mode 100644 index 0000000..dbfb0b6 --- /dev/null +++ b/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} +} diff --git a/doc-images/equations/basis-function-zero.svg b/doc-images/equations/basis-function-zero.svg new file mode 100644 index 0000000..7a28935 --- /dev/null +++ b/doc-images/equations/basis-function-zero.svg @@ -0,0 +1,262 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/equations/basis-function-zero.tex b/doc-images/equations/basis-function-zero.tex new file mode 100644 index 0000000..bde9539 --- /dev/null +++ b/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} diff --git a/doc-images/equations/basis-function.svg b/doc-images/equations/basis-function.svg new file mode 100644 index 0000000..04fdac5 --- /dev/null +++ b/doc-images/equations/basis-function.svg @@ -0,0 +1,309 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/equations/basis-function.tex b/doc-images/equations/basis-function.tex new file mode 100644 index 0000000..d27943a --- /dev/null +++ b/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} diff --git a/doc-images/equations/basis-prefactor.svg b/doc-images/equations/basis-prefactor.svg new file mode 100644 index 0000000..cd2caf1 --- /dev/null +++ b/doc-images/equations/basis-prefactor.svg @@ -0,0 +1,282 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/equations/basis-prefactor.tex b/doc-images/equations/basis-prefactor.tex new file mode 100644 index 0000000..0b93230 --- /dev/null +++ b/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} diff --git a/doc-images/equations/control-points.svg b/doc-images/equations/control-points.svg new file mode 100644 index 0000000..5843c4c --- /dev/null +++ b/doc-images/equations/control-points.svg @@ -0,0 +1,388 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/equations/control-points.tex b/doc-images/equations/control-points.tex new file mode 100644 index 0000000..c95df8d --- /dev/null +++ b/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} diff --git a/doc-images/equations/curve-deriv.tex b/doc-images/equations/curve-deriv.tex deleted file mode 100644 index fc5e42f..0000000 --- a/doc-images/equations/curve-deriv.tex +++ /dev/null @@ -1,13 +0,0 @@ -\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}^{(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}] -\end{equation} -\end{document} diff --git a/doc-images/equations/curve-deriv.tex.svg b/doc-images/equations/curve-deriv.tex.svg deleted file mode 100644 index 332b428..0000000 --- a/doc-images/equations/curve-deriv.tex.svg +++ /dev/null @@ -1,268 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/doc-images/equations/curve.svg b/doc-images/equations/curve.svg new file mode 100644 index 0000000..5c13fa4 --- /dev/null +++ b/doc-images/equations/curve.svg @@ -0,0 +1,314 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/equations/curve.tex b/doc-images/equations/curve.tex index 8153ff8..393502a 100644 --- a/doc-images/equations/curve.tex +++ b/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} diff --git a/doc-images/equations/curve.tex.svg b/doc-images/equations/curve.tex.svg deleted file mode 100644 index 613d504..0000000 --- a/doc-images/equations/curve.tex.svg +++ /dev/null @@ -1,152 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/doc-images/equations/generate-equations-images.sh b/doc-images/equations/generate-equations-images.sh new file mode 100755 index 0000000..e31f89c --- /dev/null +++ b/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 diff --git a/doc-images/equations/knots.svg b/doc-images/equations/knots.svg new file mode 100644 index 0000000..57c7812 --- /dev/null +++ b/doc-images/equations/knots.svg @@ -0,0 +1,474 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/equations/knots.tex b/doc-images/equations/knots.tex new file mode 100644 index 0000000..11fa279 --- /dev/null +++ b/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} diff --git a/src/curve/basis/mod.rs b/src/curve/basis/mod.rs new file mode 100644 index 0000000..1b8be03 --- /dev/null +++ b/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); + } +} diff --git a/src/curve/knots/mod.rs b/src/curve/knots/mod.rs index 3c153a8..03abdbd 100644 --- a/src/curve/knots/mod.rs +++ b/src/curve/knots/mod.rs @@ -1,13 +1,24 @@ -//! Knot vectors can be generated with different methods: +#![cfg_attr(feature = "doc-images", +cfg_attr(all(), +doc = ::embed_doc_image::embed_image!("eq-knots", "doc-images/equations/knots.svg")))] +//! Implements the knot vector defining the [spline basis functions][basis]. //! -//! - Uniform knots -//! - Knot-averaging -//! - De-Boor's method +//! The knot vector parametrizing the `k`-th degree curve is composed of `n+p+2 - 2k` scalar values +//! in ascending order, called 'knots'. +//! +//! ![The knot vector][eq-knots] +//! +//! The head and tail contains of `p-k+1` knots of value `0` and `1`, respectively. +//! This leaves `n-p` internal knots in the center. +//! The interval from index `i = p-k,..., n+1-k` is called 'domain'. +//! +//! Different [knot vector generation methods][methods] are available. use std::ops::MulAssign; use thiserror::Error; +use crate::curve::basis; use crate::{ curve::{parameters, parameters::Parameters, CurveError}, types::{KnotVectorDerivatives, VecD, VecDView, VecHelpers}, @@ -64,7 +75,7 @@ impl Knots { let mut Uk: Vec = Vec::with_capacity(degree + 1); Uk.push(knots); - // TODO Add multipliity check + // TODO Add multiplicity check let mut knots = Knots { Uk, p: degree, k_max: 0 }; knots.derive(); @@ -168,7 +179,6 @@ impl Knots { for k in 1..=p { // obtain the `k`-th derivative knot vector from the previous `k-1`-th derivative knot vector segment // by dropping the first and last segment - let segment_of_previous_order_knot_vector = self.Uk[k - 1].segment(1, self.len(k - 1) - 2).clone_owned(); self.Uk.push(segment_of_previous_order_knot_vector); @@ -261,43 +271,10 @@ impl Knots { let n = self.segments(); let pk = p - k; - evaluate(i, pk, n, Uk, u, k) + basis::basis(Uk, i, pk, k, n, u) } } -pub(crate) fn evaluate(i: usize, pk: usize, n: usize, Uk: &VecD, u: f64, k: usize) -> f64 { - // Err(ParameterOutOfBounds) - - if pk == 0 { - // the conditional (i == n - k) && (u == Uk(n+1)) closes the last interval [u_n,u_n+1] - if (Uk[i] <= u && u < Uk[i + 1]) || (i == n - k && u == Uk[n + 1]) { - return 1.0; - } - return 0.0; - } - - let summand1 = if Uk[i + pk] == Uk[i] { - 0.0 - } else { - let g = i; - let h = pk - 1; - (u - Uk[g]) / (Uk[g + h + 1] - Uk[g]) * evaluate(i, h, n, Uk, u, k) /* self.evaluate(k, i, h, u) */ - }; - - let summand2 = if Uk[i + 1 + pk] == Uk[i + 1] { - 0.0 - } else { - let g = i + 1; - let h = pk - 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 + pk] - u) / (Uk[g + h + 1] - Uk[g]) * evaluate(g, h, n, Uk, u, k) - }; - - summand1 + summand2 -} - fn is_valid(knots: Knots) -> bool { knots.segments() == knots.len(0) - (knots.degree() + 2) } @@ -399,7 +376,6 @@ fn rescale(knots: &mut VecD, old_lim: (f64, f64), new_lim: (f64, f64)) { #[cfg(test)] mod tests { - use approx::assert_relative_eq; use nalgebra::dvector; use rstest::rstest; @@ -461,130 +437,6 @@ mod tests { assert_eq!(knots_example(3).domain(), dvector![0.0, 0.5, 1.0]); } - #[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); - } - #[test] fn multiplicity() { let knots = Knots::new(2, dvector![0., 0., 0., 0.25, 0.5, 0.5, 0.75, 1., 1., 1.]); diff --git a/src/curve/mod.rs b/src/curve/mod.rs index badd1b4..8f4bd9b 100644 --- a/src/curve/mod.rs +++ b/src/curve/mod.rs @@ -1,3 +1,23 @@ +//#![warn(missing_docs)] +//#![warn(missing_doc_code_examples)] +#![cfg_attr(feature = "doc-images", +cfg_attr(all(), +doc = ::embed_doc_image::embed_image!("eq-curve", "doc-images/equations/curve.svg")))] +//! Implements the B-spline curve. +//! +//! A B-spline curve can be defined by +//! +//! ![B-spline curve][eq-curve] +//! +//! with the +//! - parameter `u ∈ [0,1]` defining a point on the curve, +//! - derivative order `k`, +//! - number of polynomial spline segments `n`, +//! - spline degree `p`, +//! - `k`-th derivative [knot vector][knots] `U`, +//! - `n+1-k` [spline basis function][basis] `N` of degree `p` defined by the [knot vector][knots] `U`, and +//! - `n+1-k`, `N`-dimensional [control points][points] `P`. + use embed_doc_image::embed_doc_image; use thiserror::Error; @@ -14,6 +34,7 @@ use crate::{ types::VecD, }; +pub mod basis; pub mod generation; pub mod knots; pub mod parameters; diff --git a/src/curve/parameters/mod.rs b/src/curve/parameters/mod.rs index 8bae19d..fb795c1 100644 --- a/src/curve/parameters/mod.rs +++ b/src/curve/parameters/mod.rs @@ -1,4 +1,4 @@ -//! Curve parameters can be generated with different methods: +//! Implements different parameter generation methods. //! //! - Equally spaced parameters //! - Centripetal method diff --git a/src/curve/points/mod.rs b/src/curve/points/mod.rs index 1995f67..dfcacb5 100644 --- a/src/curve/points/mod.rs +++ b/src/curve/points/mod.rs @@ -1,3 +1,16 @@ +#![cfg_attr(feature = "doc-images", +cfg_attr(all(), +doc = ::embed_doc_image::embed_image!("eq-control-points", "doc-images/equations/control-points.svg")))] +//! Implements the control points constituting the control polygon of the curve. +//! +//! The control points of the `k`-th derivative B-spline curve can be derived +//! from the zeroth order control points +//! +//! ![The control points][eq-control-points] +//! +//! Control points are generated and manipulated as part of different curve generation +//! and curve manipulation methods. + use std::ops::MulAssign; use crate::{ @@ -142,8 +155,8 @@ impl ControlPoints { return VecD::zeros(self.dimension()); } - (p - k_max + 1) as f64 / (U0[i + p + 1] - U0[i + k_max]) * - (self.derive_single_point(i + 1, k_max - 1, knots) - self.derive_single_point(i, k_max - 1, knots)) + (p - k_max + 1) as f64 / (U0[i + p + 1] - U0[i + k_max]) + * (self.derive_single_point(i + 1, k_max - 1, knots) - self.derive_single_point(i, k_max - 1, knots)) } pub fn reverse(&mut self) -> &mut Self { diff --git a/src/lib.rs b/src/lib.rs index be21eec..eca474d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2,9 +2,11 @@ //#![warn(missing_doc_code_examples)] #![cfg_attr(feature = "doc-images", cfg_attr(all(), -doc = ::embed_doc_image::embed_image!("eq-curve", "doc-images/equations/curve-deriv.tex.svg"), +doc = ::embed_doc_image::embed_image!("eq-curve", "doc-images/equations/curve.svg"), +doc = ::embed_doc_image::embed_image!("eq-knots", "doc-images/equations/knots.svg"), +doc = ::embed_doc_image::embed_image!("eq-control-points", "doc-images/equations/control-points.svg"), doc = ::embed_doc_image::embed_image!("img-curve", "doc-images/plots/manipulation/insert-before.svg")))] -//! `bsplines` is a Rust library for vectorized, N-dimensional B-spline curves and their derivatives based on +//! **bsplines** is a library for vectorized, N-dimensional B-spline curves and their derivatives based on //! [nalgebra]. //! //! ## Features @@ -21,7 +23,7 @@ doc = ::embed_doc_image::embed_image!("img-curve", "doc-images/plots/manipulatio //! - [splitting][manipulation::split] //! - [merging][manipulation::split] //! -//! ## Mathematical Definition +//! ## What are B-Splines? //! //! B-splines are parametric functions composed of piecewise polynomials with a polynomial degree `p > 0`. //! These piecewise polynomials are joined so that the parametric function is `p-1` times continuously @@ -33,19 +35,6 @@ doc = ::embed_doc_image::embed_image!("img-curve", "doc-images/plots/manipulatio //! complex-shaped and high-dimensional data, while maintaining a low polynomial degree. Because of the polynomial //! nature, all possible derivatives are accessible. //! -//! A B-spline curve `C` can be defined by -//! -//! ![The mathematical definition of a B-spline curve.][eq-curve] -//! -//! with -//! the derivative order `k`, -//! the parameter `u ∈ [0,1]` defining a point on the curve, -//! the number of polynomial spline segments `n`, -//! the spline degree `p`, -//! the `k`-th derivative knot vector `U^(k)`, -//! the `n+1` spline basis function `N_{i,p}` of degree `p` defined by the knot vector `U^(k)`, and -//! the `n+1` control points `P_i` that can be of arbitrary dimension. -//! //! ![A 2D B-Spline curve.][img-curve] //! //! Still, evaluations or spatial manipulations can be executed fast because only local polynomial segments must be diff --git a/src/manipulation/merge.rs b/src/manipulation/merge.rs index 3a6979d..575038b 100644 --- a/src/manipulation/merge.rs +++ b/src/manipulation/merge.rs @@ -21,10 +21,11 @@ use std::ops::{AddAssign, DivAssign, SubAssign}; use nalgebra::SVD; use thiserror::Error; +use crate::curve::basis::basis; use crate::{ curve, curve::{ - knots::{evaluate, is_clamped, is_normed, reversed, Knots}, + knots::{is_clamped, is_normed, reversed, Knots}, points::{ControlPoints, Points}, Curve, CurveError, }, @@ -185,7 +186,7 @@ fn calculateKV(a: &Curve) -> MatD { let mut sum = 0.; for a in m - p..=m - k { - sum += prefactor(p, a, i, k, S, &Vk[0]) * evaluate(a, p - k, m - k, &Vk[k], Vk[0][m + 1], 0); + sum += prefactor(p, a, i, k, S, &Vk[0]) * basis(&Vk[k], a, p - k, 0, m - k, Vk[0][m + 1]); } KV[(k, i - (m + 1 - p))] = sum; } @@ -206,7 +207,7 @@ fn calculateKW(b: &Curve) -> MatD { let mut sum = 0.; for b in 0..=p - k { - sum += prefactor(p, b, j, k, T, &Wk[0]) * evaluate(b, p - k, o - k, &Wk[k], Wk[0][p], 0); + sum += prefactor(p, b, j, k, T, &Wk[0]) * basis(&Wk[k], b, p - k, 0, o - k, Wk[0][p]); } KW[(k, j)] = sum; } @@ -229,7 +230,7 @@ fn calculateIV(a: &Curve) -> MatD { let mut sum = 0.; for a in m - p..=m - k { - sum += prefactor(p, a, i, k, S, &Vk[0]) * evaluate(a, p - k, m - k, &Vk[k], Vk[0][m + 1], 0); + sum += prefactor(p, a, i, k, S, &Vk[0]) * basis(&Vk[k], a, p - k, 0, m - k, Vk[0][m + 1]); } IV[(i - (m + 1 - p), k)] = sum; } @@ -252,7 +253,7 @@ fn calculateJW(b: &Curve) -> MatD { let mut sum = 0.; for b in 0..=p - k { - sum += prefactor(p, b, j, k, T, &Wk[0]) * evaluate(b, p - k, o - k, &Wk[k], Wk[0][p], 0); + sum += prefactor(p, b, j, k, T, &Wk[0]) * basis(&Wk[k], b, p - k, 0, o - k, Wk[0][p]); } JW[(j, k)] = sum; } @@ -273,7 +274,7 @@ fn calculateGV(a: &ConstrainedCurve) -> MatD { for g in 0..=mg { for i in m - p + 1..=m { - GV[(g, i - (m - p + 1))] = evaluate(i, p, m, Vk0, a.constraints.params[g], 0); + GV[(g, i - (m - p + 1))] = basis(Vk0, i, p, 0, m, a.constraints.params[g]); } } GV @@ -289,7 +290,7 @@ fn calculateHW(b: &ConstrainedCurve) -> MatD { for h in 0..=oh { for i in 0..=p - 1 { - HW[(h, i)] = evaluate(i, p, o, Wk0, b.constraints.params[h], 0); + HW[(h, i)] = basis(Wk0, i, p, 0, o, b.constraints.params[h]); } } @@ -315,7 +316,7 @@ fn calculateIpV(a: &ConstrainedCurve) -> MatD { for i in m - p + 1..=m { for g in 0..=mg { - IpV[(i - (m + 1 - p), g)] = evaluate(i, p, m, Vk0, a.constraints.params[g], 0); + IpV[(i - (m + 1 - p), g)] = basis(Vk0, i, p, 0, m, a.constraints.params[g]); } } IpV *= -0.5; @@ -333,7 +334,7 @@ fn calculateJppW(b: &ConstrainedCurve) -> MatD { for j in 0..=p_ - 1 { for h in 0..=oh_ { - JppW[(j, h)] = evaluate(j, p_, o_, Wk0, b.constraints.params[h], 0); + JppW[(j, h)] = basis(Wk0, j, p_, 0, o_, b.constraints.params[h]); } } JppW *= -0.5; @@ -362,8 +363,7 @@ fn calculateKconst(a: &Curve, b: &Curve) -> MatD { for i in m - p..=m { for a in m - p..=m - k { // TODO use point getter instead of .column(i) - sum += - prefactor(p, a, i, k, S, &Vk[0]) * evaluate(a, p - k, m - k, &Vk[k], Vk[0][m + 1], 0) * S.column(i); + sum += prefactor(p, a, i, k, S, &Vk[0]) * basis(&Vk[k], a, p - k, 0, m - k, Vk[0][m + 1]) * S.column(i); } } @@ -371,7 +371,7 @@ fn calculateKconst(a: &Curve, b: &Curve) -> MatD { for b in 0..=p - k { //sumB += -prefactor(p, b, j, k, &T, &Wk[0]) * evaluate(b, p - k, o - k, &Wk[k], Wk[0][p]) * T.row(j); // TODO use point getter instead of .column(j) - sum -= prefactor(p, b, j, k, T, &Wk[0]) * evaluate(b, p - k, o - k, &Wk[k], Wk[0][p], 0) * T.column(j); + sum -= prefactor(p, b, j, k, T, &Wk[0]) * basis(&Wk[k], b, p - k, 0, o - k, Wk[0][p]) * T.column(j); } } Kconst.column_mut(k).sub_assign(&sum); @@ -581,8 +581,8 @@ fn prefactor(p: usize, i: usize, i0: usize, k: usize, P0: &MatD, U0: &VecD) -> f } else if U0[i + p + 1] == U0[i + k] { 0. } else { - (p + 1 - k) as f64 / (U0[i + p + 1] - U0[i + k]) * - (prefactor(p, i + 1, i0, k - 1, P0, U0) - prefactor(p, i, i0, k - 1, P0, U0)) + (p + 1 - k) as f64 / (U0[i + p + 1] - U0[i + k]) + * (prefactor(p, i + 1, i0, k - 1, P0, U0) - prefactor(p, i, i0, k - 1, P0, U0)) } } else { 0. diff --git a/src/manipulation/mod.rs b/src/manipulation/mod.rs index 9bac688..91ba75a 100644 --- a/src/manipulation/mod.rs +++ b/src/manipulation/mod.rs @@ -1,3 +1,5 @@ +//! Implements different curve manipulation methods. + pub mod insert; pub mod merge; pub mod reverse;