From 1b08d957413ae8e6ab0e9dc53f3ad30eb49f1a63 Mon Sep 17 00:00:00 2001 From: Michael Heuer Date: Sun, 4 Feb 2024 23:13:22 +0100 Subject: [PATCH] initial version --- .cargo/config.toml | 2 + .editorconfig | 16 + .github/workflows/ci.yml | 45 + .gitignore | 10 + Cargo.toml | 57 + LICENSE | 201 ++++ README.md | 39 + doc-images/plots/derivatives-reversed.svg | 1066 +++++++++++++++++ doc-images/plots/derivatives.svg | 1066 +++++++++++++++++ doc-images/plots/generation/fit-loose-all.svg | 254 ++++ .../generation/fit-loose-half-penalized.svg | 232 ++++ .../plots/generation/fit-loose-half.svg | 232 ++++ doc-images/plots/generation/interpolation.svg | 254 ++++ doc-images/plots/generation/manual.svg | 254 ++++ doc-images/plots/generation/points.svg | 215 ++++ .../plots/manipulation/insert-after.svg | 214 ++++ .../plots/manipulation/insert-before.svg | 212 ++++ .../merge-after-left-end-constrained.svg | 212 ++++ .../merge-after-right-start-constrained.svg | 212 ++++ doc-images/plots/manipulation/merge-after.svg | 212 ++++ .../plots/manipulation/merge-before.svg | 219 ++++ .../plots/manipulation/reverse-after.svg | 212 ++++ .../plots/manipulation/reverse-before.svg | 212 ++++ doc-images/plots/manipulation/split-after.svg | 219 ++++ .../plots/manipulation/split-before.svg | 212 ++++ doc-images/src/main.rs | 233 ++++ doc-images/src/visualization.rs | 185 +++ rustfmt.toml | 10 + src/curve/generation/mod.rs | 151 +++ src/curve/knots/methods.rs | 331 +++++ src/curve/knots/mod.rs | 686 +++++++++++ src/curve/mod.rs | 477 ++++++++ src/curve/parameters/methods.rs | 150 +++ src/curve/parameters/mod.rs | 43 + src/curve/points/methods/fit/fixed.rs | 227 ++++ src/curve/points/methods/fit/loose.rs | 167 +++ src/curve/points/methods/fit/mod.rs | 130 ++ src/curve/points/methods/interpolation.rs | 71 ++ src/curve/points/methods/mod.rs | 2 + src/curve/points/mod.rs | 222 ++++ src/docs-header.html | 83 ++ src/lib.rs | 62 + src/manipulation/insert.rs | 170 +++ src/manipulation/merge.rs | 691 +++++++++++ src/manipulation/mod.rs | 4 + src/manipulation/reverse.rs | 24 + src/manipulation/split.rs | 218 ++++ src/types.rs | 77 ++ 48 files changed, 10493 insertions(+) create mode 100644 .cargo/config.toml create mode 100644 .editorconfig create mode 100644 .github/workflows/ci.yml create mode 100644 .gitignore create mode 100644 Cargo.toml create mode 100644 LICENSE create mode 100644 README.md create mode 100644 doc-images/plots/derivatives-reversed.svg create mode 100644 doc-images/plots/derivatives.svg create mode 100644 doc-images/plots/generation/fit-loose-all.svg create mode 100644 doc-images/plots/generation/fit-loose-half-penalized.svg create mode 100644 doc-images/plots/generation/fit-loose-half.svg create mode 100644 doc-images/plots/generation/interpolation.svg create mode 100644 doc-images/plots/generation/manual.svg create mode 100644 doc-images/plots/generation/points.svg create mode 100644 doc-images/plots/manipulation/insert-after.svg create mode 100644 doc-images/plots/manipulation/insert-before.svg create mode 100644 doc-images/plots/manipulation/merge-after-left-end-constrained.svg create mode 100644 doc-images/plots/manipulation/merge-after-right-start-constrained.svg create mode 100644 doc-images/plots/manipulation/merge-after.svg create mode 100644 doc-images/plots/manipulation/merge-before.svg create mode 100644 doc-images/plots/manipulation/reverse-after.svg create mode 100644 doc-images/plots/manipulation/reverse-before.svg create mode 100644 doc-images/plots/manipulation/split-after.svg create mode 100644 doc-images/plots/manipulation/split-before.svg create mode 100644 doc-images/src/main.rs create mode 100644 doc-images/src/visualization.rs create mode 100644 rustfmt.toml create mode 100644 src/curve/generation/mod.rs create mode 100644 src/curve/knots/methods.rs create mode 100644 src/curve/knots/mod.rs create mode 100644 src/curve/mod.rs create mode 100644 src/curve/parameters/methods.rs create mode 100644 src/curve/parameters/mod.rs create mode 100644 src/curve/points/methods/fit/fixed.rs create mode 100644 src/curve/points/methods/fit/loose.rs create mode 100644 src/curve/points/methods/fit/mod.rs create mode 100644 src/curve/points/methods/interpolation.rs create mode 100644 src/curve/points/methods/mod.rs create mode 100644 src/curve/points/mod.rs create mode 100644 src/docs-header.html create mode 100644 src/lib.rs create mode 100644 src/manipulation/insert.rs create mode 100644 src/manipulation/merge.rs create mode 100644 src/manipulation/mod.rs create mode 100644 src/manipulation/reverse.rs create mode 100644 src/manipulation/split.rs create mode 100644 src/types.rs diff --git a/.cargo/config.toml b/.cargo/config.toml new file mode 100644 index 0000000..9d53892 --- /dev/null +++ b/.cargo/config.toml @@ -0,0 +1,2 @@ +[build] +rustdocflags = ["--html-in-header", "./src/docs-header.html"] diff --git a/.editorconfig b/.editorconfig new file mode 100644 index 0000000..e1c2589 --- /dev/null +++ b/.editorconfig @@ -0,0 +1,16 @@ +# EditorConfig http://EditorConfig.org + +# top-most EditorConfig file +root = true + +# All files +[*] +charset = utf-8 +end_of_line = lf +indent_size = 4 +indent_style = space +insert_final_newline = true +trim_trailing_whitespace = true + +[*.yml] +indent_size = 2 diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..2b91cf2 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,45 @@ +name: "CI" + +on: + workflow_dispatch: + pull_request: + push: + branches: + - "main" + +jobs: + check: + name: "Check" + runs-on: "ubuntu-latest" + steps: + - uses: "actions/checkout@v3" + - uses: "dtolnay/rust-toolchain@stable" + - run: "cargo check" + + test: + name: "Test Suite" + runs-on: "ubuntu-latest" + steps: + - uses: "actions/checkout@v3" + - uses: "dtolnay/rust-toolchain@stable" + - run: "cargo test" + + fmt: + name: "Rustfmt" + runs-on: "ubuntu-latest" + steps: + - uses: "actions/checkout@v3" + - uses: "dtolnay/rust-toolchain@nightly" + with: + components: "rustfmt" + - run: "cargo +nightly fmt --all -- --check" + + clippy: + name: "Clippy" + runs-on: "ubuntu-latest" + steps: + - uses: "actions/checkout@v3" + - uses: "dtolnay/rust-toolchain@nightly" + with: + components: "clippy" + - run: "cargo clippy -- -D warnings" diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..b80da91 --- /dev/null +++ b/.gitignore @@ -0,0 +1,10 @@ +# directories +target/ +.idea + +# files +.DS_Store +**/*.rs.bk + +#/target +/Cargo.lock diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..e8d0d33 --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,57 @@ +[package] +name = "bsplines" +version = "0.0.1-alpha.0" +authors = [ + "Michael A. Heuer ", +] +description = "A B-Spline curve library built on top of nalgebra" +keywords = ["b-spline", "spline", "curve"] +categories = [ + "mathematics", + "algorithms", + "graphics", + "no-std", +] +repository = "https://github.com/michael-a-heuer/bsplines" +documentation = "https://michael-a-heuer.github.io/bsplines/bsplines/" +edition = "2021" +license = "MIT" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html +[lib] +name = "bsplines" +path = "src/lib.rs" + +[[example]] +name = "doc-images" +path = "doc-images/src/main.rs" + +[lints.rust] +unsafe_code = "forbid" +non_snake_case = "allow" +dead_code = "allow" + +[lints.clippy] +enum_glob_use = "deny" +type_complexity = "allow" + +[dependencies] +embed-doc-image = "0.1.4" +nalgebra = "0.32.3" +thiserror = "1.0.56" + + +[dev-dependencies] +rstest = "0.18.2" +approx = "0.5.1" +plotters = "0.3.5" +plotters-arrows = "0.1.0" +num-traits = "0.2.17" + +[features] +doc-images = [] + +[package.metadata.docs.rs] +rustdoc-args = [ + "--html-in-header", "./src/docs-header.html", +] diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..261eeb9 --- /dev/null +++ b/LICENSE @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "[]" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright [yyyy] [name of copyright owner] + + 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. diff --git a/README.md b/README.md new file mode 100644 index 0000000..04bd54d --- /dev/null +++ b/README.md @@ -0,0 +1,39 @@ +# Bsplines + +[![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) +[![License](https://img.shields.io/crates/l/bsplines)](https://www.apache.org/licenses/LICENSE-2.0) + +Rust library for vectorized, N-dimensional B-spline curves and their derivatives based +on [nalgebra](https://docs.rs/nalgebra/latest/nalgebra/). + +## 🚧 This Library is Under Construction 🚧 + +-[ ] Use iterators and simplify loops +-[ ] Use `thiserror` +-[ ] Refactor visibility and folder structure +-[ ] Refactor method selection and settings structs + +## Docs + +To build the docs locally, use + +```sh +cargo doc \ + --package bsplines \ + --no-deps \ + --features doc-images \ + --document-private-items \ + --open +``` + +or + +```sh +cargo doc \ + --package bsplines \ + --no-deps \ + --features doc-images \ + --open +``` + diff --git a/doc-images/plots/derivatives-reversed.svg b/doc-images/plots/derivatives-reversed.svg new file mode 100644 index 0000000..b496669 --- /dev/null +++ b/doc-images/plots/derivatives-reversed.svg @@ -0,0 +1,1066 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-4.0 + + + +-2.0 + + + +0.0 + + + +2.0 + + + +4.0 + + + +6.0 + + + +8.0 + + + + +0.0 + + + +0.1 + + + +0.2 + + + +0.3 + + + +0.4 + + + +0.5 + + + +0.6 + + + +0.7 + + + +0.8 + + + +0.9 + + + +1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/derivatives.svg b/doc-images/plots/derivatives.svg new file mode 100644 index 0000000..5c48189 --- /dev/null +++ b/doc-images/plots/derivatives.svg @@ -0,0 +1,1066 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-4.0 + + + +-2.0 + + + +0.0 + + + +2.0 + + + +4.0 + + + +6.0 + + + +8.0 + + + + +0.0 + + + +0.1 + + + +0.2 + + + +0.3 + + + +0.4 + + + +0.5 + + + +0.6 + + + +0.7 + + + +0.8 + + + +0.9 + + + +1.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/generation/fit-loose-all.svg b/doc-images/plots/generation/fit-loose-all.svg new file mode 100644 index 0000000..f01cb39 --- /dev/null +++ b/doc-images/plots/generation/fit-loose-all.svg @@ -0,0 +1,254 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/generation/fit-loose-half-penalized.svg b/doc-images/plots/generation/fit-loose-half-penalized.svg new file mode 100644 index 0000000..e64e897 --- /dev/null +++ b/doc-images/plots/generation/fit-loose-half-penalized.svg @@ -0,0 +1,232 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/generation/fit-loose-half.svg b/doc-images/plots/generation/fit-loose-half.svg new file mode 100644 index 0000000..8e82574 --- /dev/null +++ b/doc-images/plots/generation/fit-loose-half.svg @@ -0,0 +1,232 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/generation/interpolation.svg b/doc-images/plots/generation/interpolation.svg new file mode 100644 index 0000000..f01cb39 --- /dev/null +++ b/doc-images/plots/generation/interpolation.svg @@ -0,0 +1,254 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/generation/manual.svg b/doc-images/plots/generation/manual.svg new file mode 100644 index 0000000..07c6fb1 --- /dev/null +++ b/doc-images/plots/generation/manual.svg @@ -0,0 +1,254 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/generation/points.svg b/doc-images/plots/generation/points.svg new file mode 100644 index 0000000..dcad6b5 --- /dev/null +++ b/doc-images/plots/generation/points.svg @@ -0,0 +1,215 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/manipulation/insert-after.svg b/doc-images/plots/manipulation/insert-after.svg new file mode 100644 index 0000000..96ed282 --- /dev/null +++ b/doc-images/plots/manipulation/insert-after.svg @@ -0,0 +1,214 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/manipulation/insert-before.svg b/doc-images/plots/manipulation/insert-before.svg new file mode 100644 index 0000000..b5251d1 --- /dev/null +++ b/doc-images/plots/manipulation/insert-before.svg @@ -0,0 +1,212 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/manipulation/merge-after-left-end-constrained.svg b/doc-images/plots/manipulation/merge-after-left-end-constrained.svg new file mode 100644 index 0000000..d68a93d --- /dev/null +++ b/doc-images/plots/manipulation/merge-after-left-end-constrained.svg @@ -0,0 +1,212 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/manipulation/merge-after-right-start-constrained.svg b/doc-images/plots/manipulation/merge-after-right-start-constrained.svg new file mode 100644 index 0000000..f5f5189 --- /dev/null +++ b/doc-images/plots/manipulation/merge-after-right-start-constrained.svg @@ -0,0 +1,212 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/manipulation/merge-after.svg b/doc-images/plots/manipulation/merge-after.svg new file mode 100644 index 0000000..b0909a9 --- /dev/null +++ b/doc-images/plots/manipulation/merge-after.svg @@ -0,0 +1,212 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/manipulation/merge-before.svg b/doc-images/plots/manipulation/merge-before.svg new file mode 100644 index 0000000..49d770d --- /dev/null +++ b/doc-images/plots/manipulation/merge-before.svg @@ -0,0 +1,219 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/manipulation/reverse-after.svg b/doc-images/plots/manipulation/reverse-after.svg new file mode 100644 index 0000000..aefdcdc --- /dev/null +++ b/doc-images/plots/manipulation/reverse-after.svg @@ -0,0 +1,212 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/manipulation/reverse-before.svg b/doc-images/plots/manipulation/reverse-before.svg new file mode 100644 index 0000000..b5251d1 --- /dev/null +++ b/doc-images/plots/manipulation/reverse-before.svg @@ -0,0 +1,212 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/manipulation/split-after.svg b/doc-images/plots/manipulation/split-after.svg new file mode 100644 index 0000000..d021d91 --- /dev/null +++ b/doc-images/plots/manipulation/split-after.svg @@ -0,0 +1,219 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc-images/plots/manipulation/split-before.svg b/doc-images/plots/manipulation/split-before.svg new file mode 100644 index 0000000..16d222d --- /dev/null +++ b/doc-images/plots/manipulation/split-before.svg @@ -0,0 +1,212 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + +-3.0 + + + +-2.0 + + + +-1.0 + + + +0.0 + + + +1.0 + + + +2.0 + + + +3.0 + + + + + + + + + + + + + + + + + + diff --git a/doc-images/src/main.rs b/doc-images/src/main.rs new file mode 100644 index 0000000..e2bc766 --- /dev/null +++ b/doc-images/src/main.rs @@ -0,0 +1,233 @@ +use std::ops::AddAssign; + +use nalgebra::{dmatrix, dvector}; +use plotters::{prelude::*, style::full_palette::TEAL}; + +use bsplines::{ + curve::{ + generation::{ + generate, + Generation::{Interpolation, LeastSquaresFit, Manual}, + }, + knots::Generation::Uniform, + points::{ + methods::fit::{Method::LooseEnds, Penalization}, + ControlPoints, DataPoints, Points, + }, + Curve, + }, + manipulation::{ + merge::{merge, merge_from, merge_to}, + split::split_and_normalize, + }, +}; + +use crate::visualization::Limits; + +mod visualization; + +const PLOTS_DIR: &str = "doc-images/plots/"; + +const RED_50: RGBAColor = RGBAColor(255, 0, 0, 0.5); +const RED_100: RGBAColor = RGBAColor(255, 0, 0, 1.0); +const BLUE_50: RGBAColor = RGBAColor(0, 0, 255, 0.5); +const BLUE_100: RGBAColor = RGBAColor(0, 0, 255, 1.0); + +const PURPLE_50: RGBAColor = RGBAColor(200, 0, 200, 0.5); +const PURPLE_100: RGBAColor = RGBAColor(200, 0, 200, 1.0); + +fn limits() -> Limits { + Limits { min: vec![-3.0, -3.0], max: vec![3.0, 3.0] } +} + +fn scatteredDataPoints() -> DataPoints { + let dp = DataPoints::new(dmatrix![ + -2.50,-2.45,-2.15,-1.70,-1.50,-1.35,-1.20, 0.05, 0.20, 0.55, 0.65, 1.00, 1.20, 1.50, 1.75, 2.00, 2.15, 2.50; + -2.55,-2.10,-2.45,-2.60,-2.15,-2.15,-1.85,-1.20,-0.70,-0.90,-0.20, 2.00, 0.95, 1.40,-0.70,-1.90,-1.70,-2.15; + ]); + + dp +} + +fn example_spline(p: usize) -> Curve { + generate(Manual { + degree: p, + points: ControlPoints::new(dmatrix![ + -2.5,-1.5,-0.5, 1.0, 2.0, 0.0; + -2.5, 1.0,-1.5,-2.0, 1.0, 2.0; + ]), + knots: Uniform, + }) + .unwrap() +} + +fn fit_plots() { + let p = 2; + let dp = scatteredDataPoints(); + let bs_max = generate(LeastSquaresFit { + degree: p, + points: &dp, + intended_segments: dp.segments(), + method: LooseEnds, + penalization: None, + }) + .unwrap(); + + let bs_half = generate(LeastSquaresFit { + degree: p, + points: &dp, + intended_segments: dp.count() / 3, + method: LooseEnds, + penalization: None, + }) + .unwrap(); + + let bs_half_penalized = generate(LeastSquaresFit { + degree: p, + points: &dp, + intended_segments: dp.count() / 3, + method: LooseEnds, + penalization: Some(Penalization { lambda: 0.5, kappa: 2 }), + }) + .unwrap(); + + visualization::generate_2d_plot("generation/points.svg", vec![], &limits(), Some(&dp)); + visualization::generate_2d_plot("generation/fit-loose-all.svg", vec![(&bs_max, RED_100)], &limits(), Some(&dp)); + visualization::generate_2d_plot("generation/fit-loose-half.svg", vec![(&bs_half, RED_100)], &limits(), Some(&dp)); + visualization::generate_2d_plot( + "generation/fit-loose-half-penalized.svg", + vec![(&bs_half_penalized, RED_100)], + &limits(), + Some(&dp), + ); +} + +fn interpolation_plot() { + let dp = scatteredDataPoints(); + let c = generate(Interpolation { degree: 2, points: &dp }).unwrap(); + + visualization::generate_2d_plot("generation/interpolation.svg", vec![(&c, RED_100)], &limits(), Some(&dp)); +} + +fn manual_plot() { + let dp = scatteredDataPoints(); + let c = generate(Manual { degree: 2, knots: Uniform, points: ControlPoints::new(dp.matrix().clone()) }).unwrap(); + visualization::generate_2d_plot("generation/manual.svg", vec![(&c, RED_100)], &limits(), Some(&dp)); +} + +fn derivatives_plot() { + let mut bs_k0 = generate(Manual { + degree: 3, + points: ControlPoints::new(dmatrix![ + -0.25, -0.05, 0.0, 0.05, 0.25; + ]), + knots: Uniform, + }) + .unwrap(); + + let lim = Limits { min: vec![-4.4], max: vec![9.0] }; + + let bs_k1 = bs_k0.get_derivative_curve(1); + let bs_k2 = bs_k0.get_derivative_curve(2); + let bs_k3 = bs_k0.get_derivative_curve(3); + + visualization::generate_1d_plot( + "derivatives.svg", + vec![(&bs_k0, RED_100), (&bs_k1, PURPLE_100), (&bs_k2, BLUE_100), (&bs_k3, TEAL.to_rgba())], + &lim, + ); + + bs_k0.reverse(); // p = 3 CORRECT + let bs_k1_rev = bs_k0.get_derivative_curve(1); // p = 2 WRONG + let bs_k2_rev = bs_k0.get_derivative_curve(2); // p = 1 CORRECT + let bs_k3_rev = bs_k0.get_derivative_curve(3); // p = 0 WRONG + + visualization::generate_1d_plot( + "derivatives-reversed.svg", + vec![(&bs_k0, RED_100), (&bs_k1_rev, PURPLE_100), (&bs_k2_rev, BLUE_100), (&bs_k3_rev, TEAL.to_rgba())], + &lim, + ); +} + +fn insert_plots() { + let p = 2; + let mut c = example_spline(p); + let lim = Limits { min: vec![-3., -3.], max: vec![3., 3.] }; + let u = 0.8; + visualization::generate_2d_plot("manipulation/insert-before.svg", vec![(&c, RED_100)], &lim, None); + + c.insert_times(u, 1).unwrap(); + + visualization::generate_2d_plot("manipulation/insert-after.svg", vec![(&c, BLUE_100)], &lim, None); +} + +fn split_plots() { + let c = example_spline(2); + let lim = Limits { min: vec![-3., -3.], max: vec![3., 3.] }; + visualization::generate_2d_plot("manipulation/split-before.svg", vec![(&c, PURPLE_100)], &lim, None); + let (a, b) = split_and_normalize(&c, 0.5, (true, true)).unwrap(); + visualization::generate_2d_plot("manipulation/split-after.svg", vec![(&a, RED_100), (&b, BLUE_100)], &lim, None); +} + +fn reverse_plots() { + let mut c = example_spline(2); + let lim = Limits { min: vec![-3., -3.], max: vec![3., 3.] }; + visualization::generate_2d_plot("manipulation/reverse-before.svg", vec![(&c, RED_100)], &lim, None); + c.reverse(); + visualization::generate_2d_plot("manipulation/reverse-after.svg", vec![(&c, PURPLE_100)], &lim, None); +} + +fn merge_plots() { + let c = example_spline(2); + let (l_unshifted, r_unshifted) = split_and_normalize(&c, 0.5, (true, true)).unwrap(); + + let mut P_a = l_unshifted.points.matrix().clone(); + let mut P_b = r_unshifted.points.matrix().clone(); + + // Shift points + P_a.column_mut(P_a.ncols() - 1).add_assign(dvector![0.25, -0.25]); + P_b.column_mut(0).add_assign(dvector![0.25, 0.25]); + + let a = Curve::new(l_unshifted.knots.clone(), ControlPoints::new(P_a)).unwrap(); + let b = Curve::new(r_unshifted.knots.clone(), ControlPoints::new(P_b)).unwrap(); + + let lim = limits(); + + visualization::generate_2d_plot("manipulation/merge-before.svg", vec![(&a, RED_100), (&b, BLUE_100)], &lim, None); + visualization::generate_2d_plot( + "manipulation/merge-after.svg", + vec![(&merge(&a, &b).unwrap(), PURPLE_100)], + &lim, + None, + ); + visualization::generate_2d_plot( + "manipulation/merge-after-left-end-constrained.svg", + vec![(&merge_from(&a, &b).unwrap(), PURPLE_100)], + &lim, + None, + ); + visualization::generate_2d_plot( + "manipulation/merge-after-right-start-constrained.svg", + vec![(&merge_to(&a, &b).unwrap(), PURPLE_100)], + &lim, + None, + ); +} + +fn main() -> Result<(), Box> { + // Derivatives + derivatives_plot(); + + // Generation + manual_plot(); + interpolation_plot(); + fit_plots(); + + // Manipulation + reverse_plots(); + insert_plots(); + split_plots(); + merge_plots(); + + Ok(()) +} diff --git a/doc-images/src/visualization.rs b/doc-images/src/visualization.rs new file mode 100644 index 0000000..a4df797 --- /dev/null +++ b/doc-images/src/visualization.rs @@ -0,0 +1,185 @@ +use plotters::{backend::SVGBackend, chart::ChartContext, coord::types::RangedCoordf64, prelude::*}; +use plotters_arrows::TriangleArrow; + +use bsplines::curve::{ + points::{ControlPoints, DataPoints, Points}, + Curve, +}; + +use crate::PLOTS_DIR; + +const IMG_SIZE: (u32, u32) = (400, 400); +const NUM_POINTS: usize = 200; + +fn linspace(start: f64, end: f64, num_points: usize) -> Vec { + assert!(num_points > 1, "Number of points must be greater than 1"); + + let step = (end - start) / (num_points - 1) as f64; + + (0..num_points).map(|i| start + i as f64 * step).collect() +} + +pub struct Limits { + pub min: Vec, + pub max: Vec, +} + +pub fn draw_parametrized_spline_2d( + chart_context: &mut ChartContext>, + curve: &Curve, + derivative: usize, + color: RGBAColor, +) { + let point_size = 1; + let u_values = linspace(0., 1., NUM_POINTS); + let data = u_values.iter().cloned().map(|u| { + let v = &curve.evaluate_derivative(u, derivative).unwrap(); + let x = v[0]; + let y = v[1]; + (x, y) + }); + chart_context.draw_series(LineSeries::new(data, color.filled().stroke_width(point_size)).point_size(0)).unwrap(); + + let mut knot_values = curve.knots.vector().data.as_vec().clone(); + knot_values.dedup(); + let knot_data = knot_values.iter().cloned().map(|u| { + let v = &curve.evaluate_derivative(u, derivative).unwrap(); + let x = v[0]; + let y = v[1]; + (x, y) + }); + + chart_context.draw_series(knot_data.map(|point| Circle::new(point, point_size * 2, color))).unwrap(); + + if curve.degree() > 1 { + draw_arrow(chart_context, curve, derivative, color.filled(), 0.9); + } +} + +fn draw_arrow( + chart_context: &mut ChartContext>, + curve: &Curve, + derivative: usize, + style: ShapeStyle, + u: f64, +) { + let v = &curve.evaluate_derivative(u, derivative).unwrap(); + let dv = &curve.evaluate_derivative(u, derivative + 1).unwrap(); + + let x = v[0]; + let y = v[1]; + let dx = dv[0]; + let dy = dv[1]; + + let norm = (dx * dx + dy * dy).sqrt(); + let dx_n = dx / norm * 0.15; + let dy_n = dy / norm * 0.15; + + let arrow = TriangleArrow::new((x - dx_n, y - dy_n), (x, y), style).width(10).head(10); + chart_context.plotting_area().draw(&arrow).unwrap(); +} + +pub fn draw_spline_1d( + chart_context: &mut ChartContext>, + curve: &Curve, + derivative: usize, + style: ShapeStyle, + dim: usize, +) { + let u_values = linspace(0., 1., NUM_POINTS); + + let data = u_values.iter().cloned().map(|u| { + let v = &curve.evaluate_derivative(u, derivative).unwrap(); + (u, v[dim]) + }); + chart_context.draw_series(LineSeries::new(data, style).point_size(1)).unwrap(); +} + +pub fn draw_control_polygon_2d( + chart_context: &mut ChartContext>, + points: &ControlPoints, + color: RGBAColor, +) { + let point_size = 3; + chart_context + .draw_series( + points.matrix().column_iter().map(|point| Circle::new((point[0], point[1]), point_size, color.filled())), + ) + .unwrap(); + + chart_context.draw_series(LineSeries::new(points.matrix().column_iter().map(|v| (v[0], v[1])), color)).unwrap(); +} + +pub fn draw_data_points_2d( + chart_context: &mut ChartContext>, + points: &DataPoints, + connected: bool, +) { + let point_size = 3; + chart_context + .draw_series( + points.matrix().column_iter().map(|point| Circle::new((point[0], point[1]), point_size, BLACK.filled())), + ) + .unwrap(); + + if connected { + chart_context.draw_series(LineSeries::new(points.matrix().column_iter().map(|v| (v[0], v[1])), BLACK)).unwrap(); + } +} + +pub fn generate_2d_plot(filename: &str, splines: Vec<(&Curve, RGBAColor)>, limits: &Limits, data: Option<&DataPoints>) { + for (c, _) in &splines { + assert_eq!(c.dimension(), 2); + } + + let mut path = String::from(PLOTS_DIR); + path.push_str(filename); + let area = SVGBackend::new(&path, IMG_SIZE).into_drawing_area(); + area.fill(&RGBAColor(255, 255, 255, 0.85)).unwrap(); + + let mut chart_builder = ChartBuilder::on(&area); + chart_builder.margin(10).set_left_and_bottom_label_area_size(20); + + let mut chart_context = + chart_builder.build_cartesian_2d(limits.min[0]..limits.max[0], limits.min[1]..limits.max[1]).unwrap(); + chart_context.configure_mesh().draw().unwrap(); + + if let Some(dp) = data { + assert_eq!(dp.dimension(), 2); + + draw_data_points_2d(&mut chart_context, &dp, false) + } + + for (c, color) in splines { + draw_control_polygon_2d(&mut chart_context, &c.points, color); + draw_parametrized_spline_2d(&mut chart_context, c, 0, color); + } + + area.present() + .expect("Unable to write result to file, please make sure 'plotters-doc-data' dir exists under current dir"); + println!("Result has been saved to {}", path); +} + +pub fn generate_1d_plot(filename: &str, splines: Vec<(&Curve, RGBAColor)>, limits: &Limits) { + for (c, _) in &splines { + assert_eq!(c.dimension(), 1); + } + + let mut path = String::from(PLOTS_DIR); + path.push_str(filename); + let area = SVGBackend::new(&path, IMG_SIZE).into_drawing_area(); + + let mut chart_builder = ChartBuilder::on(&area); + chart_builder.margin(10).set_left_and_bottom_label_area_size(20); + + let mut chart_context = chart_builder.build_cartesian_2d(0.0..1.0, limits.min[0]..limits.max[0]).unwrap(); + chart_context.configure_mesh().draw().unwrap(); + + for (bs, color) in splines { + draw_spline_1d(&mut chart_context, bs, 0, color.filled(), 0); + } + + area.present() + .expect("Unable to write result to file, please make sure 'plotters-doc-data' dir exists under current dir"); + println!("Result has been saved to {}", path); +} diff --git a/rustfmt.toml b/rustfmt.toml new file mode 100644 index 0000000..065c114 --- /dev/null +++ b/rustfmt.toml @@ -0,0 +1,10 @@ +binop_separator = "Back" +comment_width = 120 +imports_granularity = "Crate" +max_width = 120 +reorder_imports = true +tab_spaces = 4 +trailing_semicolon = true +use_field_init_shorthand = true +use_small_heuristics = "Max" +wrap_comments = true diff --git a/src/curve/generation/mod.rs b/src/curve/generation/mod.rs new file mode 100644 index 0000000..2cd9c4e --- /dev/null +++ b/src/curve/generation/mod.rs @@ -0,0 +1,151 @@ +#![cfg_attr(feature = "doc-images", +cfg_attr(all(), +doc = ::embed_doc_image::embed_image!("points", "doc-images/plots/generation/points.svg"), +doc = ::embed_doc_image::embed_image!("manual", "doc-images/plots/generation/manual.svg"), +doc = ::embed_doc_image::embed_image!("interpolation", "doc-images/plots/generation/interpolation.svg"), +doc = ::embed_doc_image::embed_image!("fit-loose-all", "doc-images/plots/generation/fit-loose-all.svg"), +doc = ::embed_doc_image::embed_image!("fit-loose-half", "doc-images/plots/generation/fit-loose-half.svg"), +doc = ::embed_doc_image::embed_image!("fit-loose-half-penalized", "doc-images/plots/generation/fit-loose-half-penalized.svg")))] +//! Generates a curve with different methods knot into the curve. +//! +//! ## Methods +//! +//! - Manual control polygon +//! - Interpolation +//! - Least-Squares Fit +//! - fixed start and end points +//! - unpenalized/penalized +//! +//! +//! | Raw Data | Manual Control Polygon | Interpolation | +//! |:------------------------------|:--------------------------|:--------------------------| +//! | ![][points] | ![][manual] | ![][interpolation] | +//! | Scattered, 2-dimensional data points
(`$N=18`).

| Curve of degree `$p=2$` with `$n=N-1$`
segments and control points generated
directly from the data points. | Curve of degree `$p=2$` and `$n = N-1$`
segments interpolating the data points.

| +//! |
| | +//! | Least-Squares Fit | Least-Squares Fit | Penalized Least-Squares Fit | +//! | ![][fit-loose-all] | ![][fit-loose-half] | ![][fit-loose-half-penalized] | +//! | Curve of degree `$p=2$` with `$n=N-1$`
segments approximating the data points.
| Curve of degree `$p=2$` with `$n=N/3$`
segments approximating the data points.
| Curve of degree `$p=2$` with `$n=N/3$`
segments approximating the data
points penalized with `$\lambda=1$`, `$\kappa=2$`. | + +use crate::curve::{ + knots, parameters, + points::{ + methods::{fit, fit::Penalization, interpolation}, + ControlPoints, DataPoints, Points, + }, + Curve, CurveError, +}; + +#[derive()] +pub enum Generation<'a> { + Manual { + degree: usize, + points: ControlPoints, + knots: knots::Generation, + }, + Interpolation { + degree: usize, + points: &'a DataPoints, + }, + LeastSquaresFit { + degree: usize, + points: &'a DataPoints, + intended_segments: usize, + method: fit::Method, + penalization: Option, + }, +} + +/// Returns a B-Spline +/// +/// # Arguments +/// +/// * `degree` - The degree of the spline +/// +/// # Examples +/// ``` +/// use nalgebra::dmatrix; +/// use bsplines::curve::generation::{generate, Generation::Manual}; +/// use bsplines::curve::knots::Generation::Uniform; +/// use bsplines::curve::points::ControlPoints; +/// +/// // Create a coordinate matrix containing with four 2D points. +/// let points = ControlPoints::new(dmatrix![ +/// // 1 2 3 4 5 +/// -2.0,-2.0,-1.0, 0.5, 1.5; // x +/// -1.0, 0.0, 1.0, 1.0, 2.0; // y +/// ]); +/// let degree = 2; +/// let curve = generate(Manual{degree, points, knots: Uniform}).unwrap(); +/// println!("{:?}", curve.evaluate(0.5)); +/// ``` +pub fn generate(generation: Generation) -> Result { + match generation { + Generation::Manual { degree, points, knots: method } => { + // TODO more sanity checks + let n = points.segments(); + let p = degree; + + if points.segments() < degree { + return Err(CurveError::DegreeAndSegmentsMismatch { p, n }); + } + + let data_points = DataPoints::new(points.matrix().clone()); + let knots = match method { + knots::Generation::Uniform => knots::methods::uniform(p, n), + knots::Generation::Manual { knots } => Ok(knots.clone()), + knots::Generation::Method { + // TODO default methods deBoor + chord length + parameter_method, + knot_method, + } => { + let params = parameters::generate(&data_points, parameter_method); + knots::generate(p, n, ¶ms, knot_method) + } + }; + Curve::new(knots?, points) + } + Generation::Interpolation { degree, points } => { + // TODO allow other methods + uniform + let params = parameters::generate(points, parameters::Method::EquallySpaced); // TODO Piegl: ChordLength + DeBoor + let knots = knots::generate(degree, points.segments(), ¶ms, knots::Method::Uniform)?; + let points = ControlPoints::new_with_capacity( + interpolation::interpolate(&knots, points, ¶ms), + knots.degree() + 1, + ); + + Curve::new(knots, points) + } + Generation::LeastSquaresFit { degree, points, intended_segments, method, penalization } => { + // TODO allow other methods + uniform + let params = parameters::generate(points, parameters::Method::EquallySpaced); + let knots = knots::generate(degree, intended_segments, ¶ms, knots::Method::Uniform)?; + + /*let (knots, params) = match penalization { + Some(..) => { + let params = parameters::generate(&points, parameters::Method::Centripetal); + let knots = knots::generate(degree, segments, ¶ms, knots::Method::DeBoor)?; + (knots, params) + } + None => { + let params = parameters::generate(points, parameters::Method::EquallySpaced); + let knots = knots::generate(degree, segments, ¶ms, knots::Method::Uniform)?; + (knots, params) + } + };*/ + + let points = match method { + fit::Method::FixedEnds => ControlPoints::new_with_capacity( + fit::fixed::fit(&knots, points, ¶ms, penalization) + .map_err(|err| CurveError::FitError { err })?, + degree + 1, + ), + fit::Method::LooseEnds => ControlPoints::new_with_capacity( + fit::loose::fit(&knots, points, ¶ms, penalization) + .map_err(|err| CurveError::FitError { err })?, + degree + 1, + ), + }; + Curve::new(knots, points) + } + } +} diff --git a/src/curve/knots/methods.rs b/src/curve/knots/methods.rs new file mode 100644 index 0000000..fd42ee1 --- /dev/null +++ b/src/curve/knots/methods.rs @@ -0,0 +1,331 @@ +use crate::{ + curve::{knots::Knots, parameters::Parameters, CurveError}, + types::VecD, +}; + +fn inputCheck(degree: usize, segments: usize) -> Result<(), CurveError> { + match degree { + 0 => Err(CurveError::DegreeTooLow { p: degree, limit: 0 }), + degree if degree > segments => Err(CurveError::DegreeAndSegmentsMismatch { p: degree, n: segments }), + _ => Ok(()), + } +} + +/// see eq. (9.7) in `Piegl1997` +/// +/// ## Note +/// Use this method only if the control points are evenly distributed. +pub fn uniform(degree: usize, segments: usize) -> Result { + let p = degree; + let n = segments; + + inputCheck(p, n)?; + + let internal_knot_count = n - p; + + let mut U = VecD::zeros(p + n + 2); + + for i in 1..=internal_knot_count { + U[degree + i] = i as f64 / (internal_knot_count + 1) as f64 + } + + // Set the tail clamp to one. + for i in n + 1..U.len() { + U[i] = 1.; + } + + //TODO the equally spaced method should not be used in conjunction with the uniform method." + Ok(Knots::new(degree, U)) +} + +/// see eq. (9.8) in `Piegl1997` +pub fn averaging(degree: usize, segments: usize, parameters: &Parameters) -> Result { + let p = degree; + let n = segments; + + inputCheck(p, n)?; + + let internal_knot_count = n - p; + + let U_bar = parameters.vector(); + + let mut U = VecD::zeros(p + n + 2); + + for j in 1..=internal_knot_count { + let mut parameter_sum = 0.; + + let lim = j + p - 1; + for i in j..=lim { + parameter_sum += U_bar[i]; + } + U[p + j] = parameter_sum / p as f64; + } + + for j in n + 1..U.len() { + U[j] = 1.; + } + + Ok(Knots::new(degree, U)) +} + +/// see eqs. (9.68-9.69) in `Piegl1997` +/// +/// ## Note +/// If this method is used, the parameters must be generated by the chord length method +/// to prevent the resulting system of linear equations from becoming singular. +/// It guarantees that every knot span contains at least one u_bar. +/// According to deBoor, this ensures that the $$N\times N$$ cofficient matrix is positive definite and +/// well-conditioned, which is important for the least-squares fitting and interpolation of data points. +pub fn de_boor(degree: usize, segments: usize, parameters: &Parameters) -> Result { + let p = degree; + let n = segments; + + inputCheck(p, n)?; + + let U_bar = parameters.vector(); + + let internal_knot_count = n - p; + let internal_knot_spans = internal_knot_count + 1; + + let d = (n + 1) as f64 / (internal_knot_spans) as f64; // eq. 9.68 in `Piegl1997` + + let mut U = VecD::zeros(p + n + 2); + + for j in 1..=internal_knot_count { + let jd = j as f64 * d; + + // Floor `j * d` into an non-negative integer + let i = jd as usize; + + // Calculate the remainder + let alpha = jd - i as f64; + + U[degree + j] = (1f64 - alpha) * U_bar[i - 1] + alpha * U_bar[i]; + } + + // Set the tail clamp to one. + for j in n + 1..U.len() { + U[j] = 1.; + } + + Ok(Knots::new(degree, U)) +} + +#[cfg(test)] +mod tests { + use nalgebra::dvector; + + use crate::curve::parameters::methods::equally_spaced; + + use super::*; + + mod uniform { + use rstest::rstest; + + use super::*; + + #[test] + fn degree_0_errors() { + let degree = 0; + let segments = degree + 1; // data points and control points counts are chosen to be identical here + + assert!(uniform(degree, segments).is_err()); + } + + #[test] + fn degree_1() { + assert_eq!(uniform(1, 4).unwrap().Uk[0], dvector![0., 0., 0.25, 0.5, 0.75, 1., 1.]); + } + + #[test] + fn degree_2() { + assert_eq!(uniform(2, 4).unwrap().Uk[0], dvector![0., 0., 0., 1. / 3., 2. / 3., 1., 1., 1.]); + } + + #[test] + fn degree_3() { + assert_eq!(uniform(3, 4).unwrap().Uk[0], dvector![0., 0., 0., 0., 0.5, 1., 1., 1., 1.]); + } + + #[test] + fn it_creates_knot_vectors_delta0() { + // segment number and degree are equal + for i in 1..3 { + assert_eq!(uniform(i, i).unwrap().Uk[0], VecD::from_vec([vec![0.; i + 1], vec![1.; i + 1]].concat())); + } + } + + #[test] + fn it_creates_knot_vectors_delta1() { + // segment number are greater than degree by one + for i in 1..3 { + assert_eq!( + uniform(i, i + 1).unwrap().Uk[0], + VecD::from_vec([vec![0f64; i + 1], vec![0.5f64], vec![1f64; i + 1]].concat()) + ); + } + } + + #[test] + fn it_creates_knot_vectors_delta2() { + // segment number are greater than degree by two + for i in 1..3 { + assert_eq!( + uniform(i, i + 2).unwrap().Uk[0], + VecD::from_vec([vec![0f64; i + 1], vec![1f64 / 3f64, 2f64 / 3f64], vec![1f64; i + 1]].concat()) + ); + } + } + + #[rstest(degree, case(1), case(2), case(3))] + fn segments_eq_degree(degree: usize) { + let segments = degree; // data points and control points counts are chosen to be identical here + + let head = vec![0.0; degree + 1]; + let tail = vec![1.0; degree + 1]; + + assert_eq!(uniform(degree, segments).unwrap().Uk[0], VecD::from_vec([head, tail].concat())); + } + + #[rstest(degree, case(1), case(2), case(3))] + fn segments_eq_degree_plus_1(degree: usize) { + let segments = degree + 1; // data points and control points counts are chosen to be identical here + + let head = vec![0.0; degree + 1]; + let internal = vec![0.5]; + let tail = vec![1.0; degree + 1]; + + assert_eq!(uniform(degree, segments).unwrap().Uk[0], VecD::from_vec([head, internal, tail].concat())); + } + } + + mod averaging { + use rstest::rstest; + + use crate::curve::knots::methods::averaging; + + use super::*; + + #[test] + fn degree_0_errors() { + let degree = 0; + let segments = degree + 1; // data points and control points counts are chosen to be identical here + let params = equally_spaced(segments); + + assert!(averaging(degree, segments, ¶ms).is_err()); + } + + #[test] + fn degree_1() { + let segments = 4; // data points and control points counts are chosen to be identical here + let params = equally_spaced(segments); + assert_eq!(averaging(1, segments, ¶ms).unwrap().Uk[0], dvector![0., 0., 0.25, 0.5, 0.75, 1., 1.]); + } + + #[test] + fn degree_2() { + let segments = 4; // data points and control points counts are chosen to be identical here + let params = equally_spaced(segments); + assert_eq!(averaging(2, segments, ¶ms).unwrap().Uk[0], dvector![0., 0., 0., 0.375, 0.625, 1., 1., 1.]); + } + + #[test] + fn degree_3() { + let segments = 4; // data points and control points counts are chosen to be identical here + let params = equally_spaced(4); + assert_eq!(averaging(3, segments, ¶ms).unwrap().Uk[0], dvector![0., 0., 0., 0., 0.5, 1., 1., 1., 1.]); + } + + #[rstest(degree, case(1), case(2), case(3))] + fn segments_eq_degree(degree: usize) { + let segments = degree; // data points and control points counts are chosen to be identical here + let params = equally_spaced(4); + + let head = vec![0.0; degree + 1]; + let tail = vec![1.0; degree + 1]; + assert_eq!(averaging(degree, segments, ¶ms).unwrap().Uk[0], VecD::from_vec([head, tail].concat())); + } + + #[rstest(degree, case(1), case(2), case(3))] + fn segments_eq_degree_plus_1(degree: usize) { + let segments = degree + 1; // data points and control points counts are chosen to be identical here + let params = equally_spaced(segments); + + let head = vec![0.0; degree + 1]; + let internal = vec![0.5]; + let tail = vec![1.0; degree + 1]; + + assert_eq!( + averaging(degree, segments, ¶ms).unwrap().Uk[0], + VecD::from_vec([head, internal, tail].concat()) + ); + } + } + + mod de_boor { + use approx::assert_relative_eq; + use nalgebra::dvector; + use rstest::rstest; + + use super::*; + + #[test] + fn degree_0_errors() { + let degree = 0; + let segments = degree + 1; // data points and control points counts are chosen to be identical here + let params = equally_spaced(segments); + + assert!(de_boor(degree, segments, ¶ms).is_err()); + } + + #[test] + fn degree_1() { + let segments = 4; // data points and control points counts are chosen to be identical here + let params = equally_spaced(segments); + assert_eq!(de_boor(1, segments, ¶ms).unwrap().Uk[0], dvector![0., 0., 0.0625, 0.375, 0.6875, 1., 1.]); + } + + #[test] + fn degree_2() { + let segments = 4; // data points and control points counts are chosen to be identical here + let params = equally_spaced(segments); + assert_relative_eq!( + de_boor(2, segments, ¶ms).unwrap().Uk[0], + dvector![0., 0., 0., 1. / 6., 1. / 3. + 0.25, 1., 1., 1.], + epsilon = f64::EPSILON.sqrt() + ); + } + + #[test] + fn degree_3() { + let segments = 4; // data points and control points counts are chosen to be identical here + let params = equally_spaced(4); + assert_eq!(de_boor(3, segments, ¶ms).unwrap().Uk[0], dvector![0., 0., 0., 0., 0.375, 1., 1., 1., 1.]); + } + + #[rstest(degree, case(1), case(2), case(3))] + fn segments_eq_degree(degree: usize) { + let segments = degree; // data points and control points counts are chosen to be identical here + let params = equally_spaced(4); + + let head = vec![0.0; degree + 1]; + let tail = vec![1.0; degree + 1]; + assert_eq!(de_boor(degree, segments, ¶ms).unwrap().Uk[0], VecD::from_vec([head, tail].concat())); + } + + #[rstest(degree, internal_knot, case(1, 1. / 4.), case(2, 2. / 6.), case(3, 3. / 8.))] + fn segments_eq_degree_plus_1(degree: usize, internal_knot: f64) { + let segments = degree + 1; // data points and control points counts are chosen to be identical here + let params = equally_spaced(segments); + + let head = vec![0.0; degree + 1]; + let internal = vec![internal_knot]; + let tail = vec![1.0; degree + 1]; + + assert_eq!( + de_boor(degree, segments, ¶ms).unwrap().Uk[0], + VecD::from_vec([head, internal, tail].concat()) + ); + } + } +} diff --git a/src/curve/knots/mod.rs b/src/curve/knots/mod.rs new file mode 100644 index 0000000..b9b62be --- /dev/null +++ b/src/curve/knots/mod.rs @@ -0,0 +1,686 @@ +//! Knot vectors can be generated with different methods: +//! +//! - Uniform knots +//! - Knot-averaging +//! - De-Boor's method + +use std::ops::MulAssign; + +use thiserror::Error; + +use crate::{ + curve::{parameters, parameters::Parameters, CurveError}, + types::{KnotVectorDerivatives, VecD, VecDView, VecHelpers}, +}; + +pub mod methods; + +#[derive(Debug, Clone)] +pub struct Knots { + pub(crate) Uk: KnotVectorDerivatives, + pub(crate) p: usize, + pub(crate) k_max: usize, +} + +pub enum DomainKnotComparatorType { + Left, + LeftOrEqual, + RightOrEqual, + Right, +} + +pub enum Generation { + Uniform, + Manual { knots: Knots }, + Method { parameter_method: parameters::Method, knot_method: Method }, +} + +pub enum Method { + Uniform, + DeBoor, + Averaging, +} + +#[derive(Error, Debug, PartialEq)] +pub enum KnotError { + #[error("Parameter `u = {u}` lies outside the interval `[{lower_bound}, {upper_bound}]`.")] + ParameterOutOfBounds { u: f64, lower_bound: f64, upper_bound: f64 }, + // TODO needed? + //#[error("Knot `u = {u}` has a muliplicity of {multiplicity}.")] + //InvalidMultiplicity { u: f64, multiplicity: usize }, +} + +pub fn generate(degree: usize, segments: usize, params: &Parameters, method: Method) -> Result { + match method { + Method::Uniform => methods::uniform(degree, segments), + Method::DeBoor => methods::de_boor(degree, segments, params), + Method::Averaging => methods::averaging(degree, segments, params), + } +} + +impl Knots { + // generates params automatically, if none are provided. + pub fn new(degree: usize, knots: VecD) -> Self { + let mut Uk: Vec = Vec::with_capacity(degree + 1); + Uk.push(knots); + + // TODO Add multipliity check + + let mut knots = Knots { Uk, p: degree, k_max: 0 }; + knots.derive(); + knots + } + + pub fn vector(&self) -> &VecD { + &self.Uk[0] + } + + pub fn vector_mut(&mut self) -> &mut VecD { + &mut self.Uk[0] + } + + /// # Arguments + /// * `k` - The `$k$`-th derivative knot vector. + pub fn vector_derivative(&self, derivative: usize) -> &VecD { + &self.Uk[derivative] + } + + pub fn vector_derivative_mut(&mut self, derivative: usize) -> &mut VecD { + &mut self.Uk[derivative] + } + + pub fn degree(&self) -> usize { + self.p + } + + pub fn segments(&self) -> usize { + self.segments_derivative(0) + } + + pub fn segments_derivative(&self, k: usize) -> usize { + self.Uk[k].len() - (self.p + 2 - 2 * k) + } + + pub fn len(&self, k: usize) -> usize { + self.Uk[k].len() + } + + pub fn internal_count(&self) -> usize { + self.segments() - self.p + } + + pub fn internal(&self) -> VecDView { + self.Uk[0].segment(self.p + 1, self.internal_count()) + } + + pub fn internal_knot(&self, i: usize) -> f64 { + self.internal()[i] + } + + pub fn domain_count(&self) -> usize { + self.segments() - self.p + 2 + } + + pub fn domain(&self) -> VecDView { + self.domain_derivative(0) + } + + pub fn domain_derivative(&self, k: usize) -> VecDView { + self.Uk[k].segment(self.p - k, self.domain_count()) + } + + pub fn domain_knot(&self, i: usize) -> f64 { + self.domain()[i] + } + + pub fn multiplicity(&self, u: f64) -> usize { + self.domain().iter().filter(|&x| *x == u).count() + } + + pub(crate) fn reverse(&mut self) -> &mut Self { + for knots in self.Uk.iter_mut() { + reverse(knots); + } + self + } + + pub fn normalize(&mut self) -> &mut Self { + for knots in self.Uk.iter_mut() { + normalize(knots); + } + self + } + + pub fn rescale(&mut self, old_lim: (f64, f64), new_lim: (f64, f64)) { + for knots in self.Uk.iter_mut() { + rescale(knots, old_lim, new_lim); + } + } + + pub fn max_derivative(&self) -> usize { + self.k_max + } + + pub fn derive(&mut self) { + let p = self.p; + + self.Uk.truncate(1); + 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); + } + self.k_max = p; + } + + /// Returns the index `i` of the knot on the domain interval + /// `$[u_{p-k}^{(k)},u_{n+1-k}^{(k)})$`, + /// being less or equal than `u`. + fn find_index(&self, u: f64, k: usize) -> Option { + let Uk = self.vector_derivative(k); + let pk = self.degree() - k; + + let first = pk; + + if u < Uk[first] { + return None; + } + + let last = Uk.len() - pk - 1; + let mut low = first; + let mut high = last; + + while low < high { + let mid = high - (high - low) / 2; + + if Uk[mid] > u { + high = mid - 1; + } else { + low = mid; + } + } + Some(high) + } + + /// Returns the index `i` of the knot on the domain interval + /// `$[u_{p-k}^{(k)},u_{n+1-k}^{(k)}]$`, + /// being lower, equal, or higher than `u`. + pub fn find_idx(&self, u: f64, k: usize, comparator: DomainKnotComparatorType) -> usize { + let Uk = self.vector_derivative(k); + let pk = self.degree() - k; + match comparator { + DomainKnotComparatorType::Left => { + let lim = self.segments() + 1 - k; + let mut i = pk; + + while u > Uk[i + 1] && i + 1 < lim { + i += 1; + } + i + } + DomainKnotComparatorType::LeftOrEqual => { + let lim = self.segments() + 1 - k; + let mut i = pk; + + while u >= Uk[i + 1] && i + 1 < lim { + i += 1; + if Uk[i + 1] == Uk[i] { + break; + } + } + i + } + DomainKnotComparatorType::RightOrEqual => { + let mut i = Uk.len() - 1 - pk; + + while u <= Uk[i - 1] && i > pk { + i -= 1; + if Uk[i - 1] == Uk[i] { + break; + } + } + i + } + DomainKnotComparatorType::Right => { + let mut i = Uk.len() - 1 - pk; + + while u < Uk[i - 1] && i - 1 > pk { + i -= 1; + } + i + } + } + } + + /// `p` the degree of this basis function of the kth degree spline - not of the 0th degree spline + pub fn evaluate(&self, k: usize, i: usize, p: usize, u: f64) -> f64 { + let Uk = &self.Uk[k]; + let n = self.segments(); + let pk = p - k; + + evaluate(i, pk, n, Uk, u, k) + } +} + +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) +} + +pub fn is_clamped(knots: &Knots) -> bool { + let U0 = knots.vector(); + let clamp_size = knots.p + 1; + + let is_head_clamped = U0.iter().take(clamp_size).all(|&u| u == 0.0); + let is_tail_clamped = U0.iter().rev().take(clamp_size).all(|&u| u == 1.0); + + is_head_clamped && is_tail_clamped +} + +pub fn is_normed(knots: &Knots) -> bool { + let U0 = knots.vector(); + + let is_min_zero = U0.iter().min_by(|a, b| a.partial_cmp(b).unwrap()) == Some(&0.0); + let is_max_unity = U0.iter().max_by(|a, b| a.partial_cmp(b).unwrap()) == Some(&1.0); + + is_min_zero && is_max_unity +} + +pub fn is_sorted(knots: &Knots) -> bool { + let mut it = knots.Uk[0].iter(); + match it.next() { + None => true, + Some(first) => it + .scan(first, |state, next| { + let cmp = *state <= next; + *state = next; + Some(cmp) + }) + .all(|b| b), + } +} + +pub fn is_uniform(knots: &Knots) -> Result { + let U0 = knots.vector(); + + let expected = methods::uniform(knots.degree(), knots.segments())?; + + Ok(U0.eq(expected.vector())) + // U0.relative_eq(expected.vector(None), NUMERICAL_PRECISION, 0.0) // TODO test +} + +pub(crate) fn reverse(knots: &mut VecD) { + let nrows = knots.nrows(); + let half_nrows = knots.len() / 2; + + for i in 0..half_nrows { + knots.swap_rows(i, nrows - 1 - i); + } + + knots.add_scalar_mut(-1.0); + knots.mul_assign(-1.0); +} + +pub(crate) fn reversed(knots: &VecD) -> VecD { + let mut copy = knots.clone(); + reverse(&mut copy); + copy +} + +pub fn normalize(knots: &mut VecD) { + let old_lim = (knots.min(), knots.max()); + + rescale(knots, old_lim, (0.0, 1.0)) +} + +pub fn normalized(knots: &mut VecD) -> VecD { + let mut copy = knots.clone(); + normalize(&mut copy); + copy +} + +fn rescaled_knot(mut knot: f64, old_lim: (f64, f64), new_lim: (f64, f64)) -> f64 { + knot -= old_lim.0; + knot /= old_lim.1 - old_lim.0; + knot *= new_lim.1 - new_lim.0; + knot += new_lim.0; + + knot +} + +fn rescaled(knots: &VecD, old_lim: (f64, f64), new_lim: (f64, f64)) -> VecD { + let mut copy = knots.clone(); + rescale(&mut copy, old_lim, new_lim); + copy +} + +fn rescale(knots: &mut VecD, old_lim: (f64, f64), new_lim: (f64, f64)) { + let n = knots.len(); + *knots -= VecD::repeat(n, old_lim.0); + *knots /= old_lim.1 - old_lim.0; + *knots *= new_lim.1 - new_lim.0; + *knots += VecD::repeat(n, new_lim.0); +} + +#[cfg(test)] +mod tests { + use approx::assert_relative_eq; + use nalgebra::dvector; + use rstest::rstest; + + use super::*; + + const SEGMENTS: usize = 4; + + fn knots_example(degree: usize) -> Knots { + methods::uniform(degree, SEGMENTS).unwrap() + } + + #[rstest(degree, case(1), case(2), case(3))] + fn segments(degree: usize) { + assert_eq!(knots_example(degree).segments(), SEGMENTS); + } + + #[rstest(degree, expected, case(1, 3), case(2, 2), case(3, 1))] + fn internal_count(degree: usize, expected: usize) { + assert_eq!(knots_example(degree).internal_count(), expected); + } + + #[test] + fn internal_degree_1() { + assert_eq!(knots_example(1).internal(), dvector![0.25, 0.5, 0.75]); + } + + #[test] + fn internal_degree_2() { + assert_eq!(knots_example(2).internal(), dvector![1.0 / 3.0, 2.0 / 3.0]); + } + + #[test] + fn internal_degree_3() { + assert_eq!(knots_example(3).internal(), dvector![0.5]); + } + + #[rstest(degree, expected, case(1, 5), case(2, 4), case(3, 3))] + fn domain_count(degree: usize, expected: usize) { + assert_eq!(knots_example(degree).domain_count(), expected); + } + + #[test] + fn domain_degree_1() { + assert_eq!(knots_example(1).domain(), dvector![0.0, 0.25, 0.5, 0.75, 1.0]); + } + + #[test] + fn domain_degree_2() { + assert_eq!(knots_example(2).domain(), dvector![0.0, 1.0 / 3.0, 2.0 / 3.0, 1.0]); + } + + #[test] + fn domain_degree_2_derivative() { + assert_eq!(knots_example(2).domain_derivative(1), dvector![0.0, 1.0 / 3.0, 2.0 / 3.0, 1.0]); + } + + #[test] + fn domain_degree_3() { + 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.]); + + assert_eq!(knots.multiplicity(0.2), 0); + assert_eq!(knots.multiplicity(0.25), 1); + assert_eq!(knots.multiplicity(0.5), 2); + assert_eq!(knots.multiplicity(0.), 1); + assert_eq!(knots.multiplicity(1.), 1); + } + + #[test] + fn knot_derivation() { + let degree = 3; + let knots = methods::uniform(degree, SEGMENTS).unwrap(); + + let U0 = dvector![0.0, 0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0, 1.0]; + let U1 = dvector![0.0, 0.0, 0.0, 0.5, 1.0, 1.0, 1.0]; + let U2 = dvector![0.0, 0.0, 0.5, 1.0, 1.0]; + let U3 = dvector![0.0, 0.5, 1.0]; + + assert_eq!(knots.vector_derivative(0), &U0); + assert_eq!(knots.vector_derivative(1), &U1); + assert_eq!(knots.vector_derivative(2), &U2); + assert_eq!(knots.vector_derivative(3), &U3); + } + + #[test] + fn normalize() { + let mut knots = Knots::new(1, dvector![1.0, 1.0, 1.5, 2.0, 2.0]); + knots.normalize(); + assert_eq!(knots.vector(), &dvector![0.0, 0.0, 0.5, 1.0, 1.0]); + } + + #[test] + fn reverse() { + let mut knots = Knots::new(1, dvector![0.0, 0.0, 0.6, 1.0, 1.0]); + knots.reverse(); + assert_eq!(knots.vector(), &dvector![0.0, 0.0, 0.4, 1.0, 1.0]); + } + + #[test] + fn is_sorted_test() { + assert!(is_sorted(&Knots::new(1, dvector![0.0, 0.0, 0.5, 1.0, 1.0]))); + assert!(!is_sorted(&Knots::new(1, dvector![0.0, 1.0, 0.5, 1.0, 1.0]))); + } + + #[test] + fn is_clamped_test() { + assert!(is_clamped(&Knots::new(1, dvector![0.0, 0.0, 0.5, 1.0, 1.0]))); + assert!(!is_uniform(&Knots::new(1, dvector![0.0, 1.0, 0.5, 1.0, 1.0])).unwrap()); + //TODO assert_eq!(is_clamped(&Knots::new(1, dvector![2.0, 2.0, 0.5, 3.0, 3.0])), true); + } + + #[test] + fn is_normed_test() { + assert!(is_normed(&Knots::new(1, dvector![0.0, 0.0, 0.5, 1.0, 1.0]))); + assert!(!is_normed(&Knots::new(1, dvector![0.0, 0.0, 1.5, 1.0, 1.0]))); + } + + #[test] + fn is_uniform_test() { + assert!(is_uniform(&Knots::new(1, dvector![0.0, 0.0, 0.25, 0.5, 0.75, 1.0, 1.0])).unwrap()); + assert!(!is_uniform(&Knots::new(1, dvector![0.0, 0.0, 0.25, 0.75, 1.0, 1.0])).unwrap()); + } + + #[rstest(u, expected, case(0.24, 1), case(0.25, 1), case(0.26, 2), case(0.74, 3), case(0.75, 3), case(0.76, 4))] + fn test_find_idx_of_left_domain_knot(u: f64, expected: usize) { + assert_eq!(knots_example(1).find_idx(u, 0, DomainKnotComparatorType::Left), expected); + } + + #[rstest(u, expected, case(0.24, 1), case(0.25, 2), case(0.26, 2), case(0.74, 3), case(0.75, 4), case(0.76, 4))] + fn test_find_idx_of_left_or_equal_domain_knot(u: f64, expected: usize) { + assert_eq!(knots_example(1).find_idx(u, 0, DomainKnotComparatorType::LeftOrEqual), expected); + } + + #[rstest(u, expected, case(0.24, 1), case(0.25, 2), case(0.26, 2), case(0.74, 3), case(0.75, 4), case(0.76, 4))] + fn test_find_idx_of_left_or_equal_domain_knot_bisection(u: f64, expected: usize) { + assert_eq!(knots_example(1).find_index(u, 0).unwrap(), expected); + } + + #[test] + fn test_find_idx_of_left_or_equal_domain_knot_bisection_limits() { + assert_eq!(knots_example(1).find_index(-0.1, 0), None); + assert_eq!(knots_example(1).find_index(0.0, 0), Some(1)); + assert_eq!(knots_example(1).find_index(1.0, 0), Some(5)); + assert_eq!(knots_example(1).find_index(1.2, 0), Some(5)); + } + + #[rstest(u, expected, case(0.24, 2), case(0.25, 2), case(0.26, 3), case(0.74, 4), case(0.75, 4), case(0.76, 5))] + fn test_find_idx_of_right_or_equal_domain_knot(u: f64, expected: usize) { + assert_eq!(knots_example(1).find_idx(u, 0, DomainKnotComparatorType::RightOrEqual), expected); + } + + #[rstest(u, expected, case(0.24, 2), case(0.25, 3), case(0.26, 3), case(0.74, 4), case(0.75, 5), case(0.76, 5))] + fn test_find_idx_of_right_domain_knot(u: f64, expected: usize) { + assert_eq!(knots_example(1).find_idx(u, 0, DomainKnotComparatorType::Right), expected); + } +} diff --git a/src/curve/mod.rs b/src/curve/mod.rs new file mode 100644 index 0000000..954df82 --- /dev/null +++ b/src/curve/mod.rs @@ -0,0 +1,477 @@ +use embed_doc_image::embed_doc_image; +use thiserror::Error; + +use crate::{ + curve::{ + knots::Knots, + points::{methods::fit::FitError, ControlPoints, Points}, + }, + manipulation::{ + insert::{insert, InsertError}, + merge::{merge, merge_with_constraints, ConstrainedCurve, Constraints}, + split::{split, SplitError}, + }, + types::VecD, +}; + +pub mod generation; +pub mod knots; +pub mod parameters; +pub mod points; + +#[embed_doc_image("spline", "doc-images/plots/derivatives.svg")] +#[derive(Debug, Clone)] +pub struct Curve { + pub knots: Knots, + pub points: ControlPoints, +} + +#[derive(Error, Debug, PartialEq)] +pub enum CurveError { + #[error("Parameter `u = {u}` lies outside the interval `[{lower_bound}, {upper_bound}]`.")] + ParameterOutOfBounds { u: f64, lower_bound: f64, upper_bound: f64 }, + #[error( + "The number of polynomial segments `n = {n}` of the curve must be \ + greater than or equal to its polynomial degree `p = {p}." + )] + DegreeAndSegmentsMismatch { p: usize, n: usize }, + + #[error("The derivative order `k = {k}` cannot be greater than curve degree `p = {p}.")] + DegreeAndDerivativeOrderMismatch { p: usize, k: usize }, + + #[error("The maximal derivative order k_max = {k_max} cannot be greater than the spline degree p = {p}.")] + DerivativeNotAvailable { k_max: usize, p: usize }, + + #[error("Degree `p = {p}` is too low and must be greater than `{limit}`")] + DegreeTooLow { p: usize, limit: usize }, + + #[error("Curve generation failed with error {err}.")] + FitError { err: FitError }, + // TODO + //#[error("Knot generation failed with error {err}.")] + //KnotError { err: KnotError }, +} + +impl Curve { + /// Returns a B-Spline + /// + /// # Arguments + /// + /// * `degree` - The degree of the spline + /// + /// # Examples + /// ``` + /// use nalgebra::{dmatrix, dvector}; + /// use bsplines::curve::Curve; + /// use bsplines::curve::generation::{generate, Generation::Manual}; + /// use bsplines::curve::knots; + /// use bsplines::curve::points::ControlPoints; + /// + /// // Create a coordinate matrix containing with four 2D points. + /// let points = ControlPoints::new(dmatrix![ + /// // 1 2 3 4 5 + /// -2.0,-2.0,-1.0, 0.5, 1.5; // x + /// -1.0, 0.0, 1.0, 1.0, 2.0; // y + /// ]); + /// let degree = 2; + /// let knots = knots::methods::uniform(degree, points.segments()).unwrap(); + /// let curve = Curve::new(knots, points).unwrap(); + /// println!("{:?}", curve.evaluate(0.5)); + /// ``` + pub fn new(knots: Knots, points: ControlPoints) -> Result { + // TODO more sanity checks + + match (knots.degree(), points.segments()) { + (p, n) if n < p => Err(CurveError::DegreeAndSegmentsMismatch { p, n }), + _ => { + let mut c = Self { knots, points }; + c.calculate_derivatives(); + Ok(c) + } + } + } + + pub fn degree(&self) -> usize { + self.knots.degree() + } + + pub fn segments(&self) -> usize { + self.points.segments() + } + + /// Returns the dimension of the curve. + pub fn dimension(&self) -> usize { + self.points.dimension() + } + + pub fn evaluate(&self, u: f64) -> Result { + self.evaluate_derivative(u, 0) + } + + pub fn evaluate_derivative(&self, u: f64, k: usize) -> Result { + if !(0.0..=1.0).contains(&u) { + return Err(CurveError::ParameterOutOfBounds { u, lower_bound: 0.0, upper_bound: 1.0 }); + } + + let p = self.degree(); + + let mut value = VecD::zeros(self.points.dimension()); + + if k <= p { + let n = self.segments(); + let l = self.knots.find_idx(u, k, knots::DomainKnotComparatorType::LeftOrEqual); + + for i in l - (p - k)..=n - k { + value += self.knots.evaluate(k, i, p, u) * self.points.matrix_derivative(k).column(i); + } + } + Ok(value) + } + + pub fn max_derivative(&self) -> usize { + let kmax_knots = self.knots.max_derivative(); + let kmax_points = self.points.max_derivative(); + + assert_eq!( + kmax_knots, kmax_points, + "The available derivatives of the knots and control points do not match up {} != {}.", + kmax_knots, kmax_points + ); + + kmax_knots + } + + /// Reverses the curve. + pub fn reverse(&mut self) -> &mut Self { + self.knots.reverse(); + self.points.reverse(); + self + } + + /// Prepends another curve. + /// + /// The end of another curve is attached to the beginning of this curve, + /// while maintaining continuity of all derivatives + /// This affects the first and last `$p$` control points of the two curves, respectively, + /// and removes `$p$` control points in total. + /// + /// # Examples + /// + /// ``` + /// use approx::relative_eq; + /// use nalgebra::dmatrix; + /// use bsplines::curve::Curve; + /// use bsplines::curve::generation::generate; + /// use bsplines::curve::generation::Generation::Manual; + /// use bsplines::curve::knots::Generation::Uniform; + /// use bsplines::curve::points::{ControlPoints, Points}; + /// + /// let mut a = generate(Manual { + /// degree: 2, + /// points: ControlPoints::new(dmatrix![ 1.0, 2.0, 3.0;]), + /// knots: Uniform, + /// }).unwrap(); + /// let b = generate(Manual { + /// degree: 2, + /// points: ControlPoints::new(dmatrix![-3.0,-2.0,-1.0;]), + /// knots: Uniform, + /// }).unwrap(); + /// let c = a.prepend(&b); + /// + /// relative_eq!(c.points.matrix(), &dmatrix![-3.0,-2.0, 2.0, 3.0;], epsilon = f64::EPSILON); + /// ``` + pub fn prepend(&mut self, other: &Self) -> &mut Self { + let c = merge(other, self).unwrap(); + self.knots = c.knots; + self.points = c.points; + self + } + + /// Prepends another curve with maximally `$p-1$` constraints. + pub fn prepend_constrained( + &mut self, + constraints_self: Constraints, + other: &Self, + constraints_other: Constraints, + ) -> &mut Self { + let c = merge_with_constraints( + &ConstrainedCurve { curve: other, constraints: constraints_self }, + &ConstrainedCurve { curve: self, constraints: constraints_other }, + ) + .unwrap(); + self.knots = c.knots; + self.points = c.points; + self + } + + /// Appends another curve. + /// + /// The end of this curve is attached to the beginning of the other curve, + /// while maintaining continuity of all derivatives + /// This affects the first and last `$p$` control points of the two curves, respectively, + /// and removes `$p$` control points in total. + /// + /// # Examples + /// + /// ``` + /// use approx::relative_eq; + /// use nalgebra::dmatrix; + /// use bsplines::curve::Curve; + /// use bsplines::curve::generation::generate; + /// use bsplines::curve::generation::Generation::Manual; + /// use bsplines::curve::knots::Generation::Uniform; + /// use bsplines::curve::points::{ControlPoints, Points}; + /// + /// let mut a = generate(Manual { + /// degree: 2, + /// points: ControlPoints::new(dmatrix![-3.0,-2.0,-1.0;]), + /// knots: Uniform, + /// }).unwrap(); + /// let b = generate(Manual { + /// degree: 2, + /// points: ControlPoints::new(dmatrix![ 1.0, 2.0, 3.0;]), + /// knots: Uniform, + /// }).unwrap(); + /// + /// let c = a.append(&b); + /// + /// relative_eq!(c.points.matrix(), &dmatrix![-3.0,-2.0, 2.0, 3.0;], epsilon = f64::EPSILON); + /// ``` + pub fn append(&mut self, other: &Self) -> &mut Self { + let c = merge(self, other).unwrap(); + self.knots = c.knots; + self.points = c.points; + self + } + + /// Appends another curve with maximally `$p-1$` constraints. + pub fn append_constrained( + &mut self, + constraints_self: Constraints, + other: &Self, + constraints_other: Constraints, + ) -> &mut Self { + let c = merge_with_constraints( + &ConstrainedCurve { curve: self, constraints: constraints_self }, + &ConstrainedCurve { curve: other, constraints: constraints_other }, + ) + .unwrap(); + self.knots = c.knots; + self.points = c.points; + self + } + + /// Splits a curve into two at parameter `$u$`. + /// + /// # Arguments + /// * `u` - The parameter`$u$` that must lie in the interval `$(0,1)$`. + pub fn split(&self, u: f64) -> Result<(Self, Self), SplitError> { + split(self, u) + } + + /// Inserts a knot into the curve at parameter `$u$`. + /// + /// # Arguments + /// * `u` - The parameter`$u$` that must lie in the interval `$(0,1)$`. + pub fn insert(&mut self, u: f64) -> Result<&mut Self, InsertError> { + self.insert_times(u, 1)?; + Ok(self) + } + + /// Inserts a knot `x` times into the curve at parameter `$u$`. + /// + /// # Arguments + /// * `u` - The parameter`$u$` that must lie in the interval `$(0,1)$`. + /// * `x` - The number of insertions of parameter `$u$`. + pub fn insert_times(&mut self, u: f64, x: usize) -> Result<&mut Self, InsertError> { + for _ in 0..x { + insert(self, u)?; + } + Ok(self) + } + + pub(crate) fn calculate_derivatives(&mut self) { + self.knots.derive(); + self.points.derive(&self.knots); + } + + /// Returns the curve describing the `$k$`-th derivative of the current one. + /// # Arguments + /// + /// * `k` - The derivative + pub fn get_derivative_curve(&self, k: usize) -> Self { + let knots = Knots::new(self.degree() - k, self.knots.vector_derivative(k).clone()); + let points = ControlPoints::new(self.points.matrix_derivative(k).clone()); + Curve { knots, points } + } +} + +#[cfg(test)] +mod tests { + use approx::assert_relative_eq; + use nalgebra::{dmatrix, dvector}; + use rstest::fixture; + + use points::DataPoints; + + use crate::curve::{ + generation::{ + generate, + Generation::{Interpolation, Manual}, + }, + knots::Generation::Uniform, + }; + + use super::*; + + #[fixture] + /// A one-dimensional, linear test curve with default degree two. + fn c(#[default(2)] degree: usize) -> Curve { + let c = generate(Manual { + degree, + points: ControlPoints::new(dmatrix![ + 1., 3., 5.; + 2., 4., 6.; + ]), + knots: Uniform, + }) + .unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 1., 1., 1.]); + c + } + + mod evaluate { + use rstest::rstest; + + use super::*; + + #[rstest] + fn non_existing_derivative(c: Curve) { + let k = 3; + assert_eq!(c.evaluate_derivative(0.5, k), Ok(dvector![0., 0.])); + } + + #[rstest] //TODO test for orders > 0 + fn outside_lower_bound(c: Curve) { + let u = -0.1; + assert_eq!( + c.evaluate_derivative(u, 0), + Err(CurveError::ParameterOutOfBounds { u, lower_bound: 0.0, upper_bound: 1.0 }) + ); + } + + #[rstest] + fn outside_upper_bound(c: Curve) { + let u = 1.1; + assert_eq!( + c.evaluate_derivative(u, 0), + Err(CurveError::ParameterOutOfBounds { u, lower_bound: 0.0, upper_bound: 1.0 }) + ); + } + + #[test] + fn derivative_out_of_bounds_error() { + let p = 2; + let points = dmatrix![1., 1., 1., 1.;]; + let mut c = generate(Manual { degree: p, points: ControlPoints::new(points), knots: Uniform }).unwrap(); + + c.insert(0.5).unwrap(); + + assert_eq!(c.evaluate_derivative(0.9, 1).unwrap(), dvector![0.]); + } + + #[test] + fn evaluate_p_repeated_knots() { + let p = 3; + let points = dmatrix![-1., -0.5, 0.5, 1.;]; + let mut c = generate(Manual { degree: p, points: ControlPoints::new(points), knots: Uniform }).unwrap(); + let u = 0.5; + let expectedEvaluationResult = dvector![0.0]; + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 0., 1., 1., 1., 1.]); + assert_eq!(c.evaluate(0.0).unwrap(), dvector![-1.]); + assert_eq!(c.evaluate(1.0).unwrap(), dvector![1.]); + + insert(&mut c, u).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 0., u, 1., 1., 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1., -0.75, 0.0, 0.75, 1.;]); + assert_eq!(c.evaluate(u).unwrap(), expectedEvaluationResult); + assert_eq!(c.evaluate(0.0).unwrap(), dvector![-1.]); + assert_eq!(c.evaluate(1.0).unwrap(), dvector![1.]); + + insert(&mut c, u).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 0., u, u, 1., 1., 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1., -0.75, -0.375, 0.375, 0.75, 1.;]); + //assert_eq!(c.evaluate_naive(u).unwrap(), expectedEvaluationResult); + assert_eq!(c.evaluate(u).unwrap(), expectedEvaluationResult); + assert_eq!(c.evaluate(0.0).unwrap(), dvector![-1.]); + assert_eq!(c.evaluate(1.0).unwrap(), dvector![1.]); + + insert(&mut c, u).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 0., u, u, u, 1., 1., 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1., -0.75, -0.375, 0.0, 0.375, 0.75, 1.;]); + //assert_eq!(c.evaluate_naive(u).unwrap(), expectedEvaluationResult); + assert_eq!(c.evaluate(u).unwrap(), expectedEvaluationResult); + assert_eq!(c.evaluate(0.0).unwrap(), dvector![-1.]); + assert_eq!(c.evaluate(1.0).unwrap(), dvector![1.]); + } + + #[rstest] + fn start(c: Curve) { + assert_eq!(c.evaluate_derivative(0., 0).unwrap(), dvector![1., 2.]) + } + + #[rstest] + fn middle(c: Curve) { + assert_eq!(c.evaluate_derivative(0.5, 0).unwrap(), dvector![3., 4.]) + } + + #[rstest] + fn end(c: Curve) { + assert_eq!(c.evaluate_derivative(1., 0).unwrap(), dvector![5., 6.]) + } + } + + #[test] + fn reverse() { + let mut c = generate(Manual { + degree: 2, + points: ControlPoints::new(dmatrix![ + 1., 3., 5.; + 2., 4., 6.; + ]), + knots: Uniform, + }) + .unwrap(); + + let knots_before = c.knots.vector().clone(); + let points_before = c.points.matrix().clone(); + c.reverse(); + + let points_after = dmatrix![ + 5., 3., 1.; + 6., 4., 2.; + ]; + assert_eq!(c.knots.vector(), &knots_before); + assert_eq!(c.points.matrix(), &points_after); + + c.reverse(); + + assert_eq!(c.knots.vector(), &knots_before); + assert_eq!(c.points.matrix(), &points_before); + } + + #[test] + fn interpolate_linear() { + let points = DataPoints::new(dmatrix![ + 1., 2., 3., 4.; + 1., 2., 3., 4.; + ]); + + let c = generate(Interpolation { degree: 1, points: &points }).unwrap(); + + assert_eq!(c.evaluate(0.0).unwrap(), dvector![1., 1.]); + assert_relative_eq!(c.evaluate(1. / 3.).unwrap(), dvector![2., 2.], epsilon = f64::EPSILON.sqrt()); + assert_eq!(c.evaluate(2. / 3.).unwrap(), dvector![3., 3.]); + assert_eq!(c.evaluate(1.0).unwrap(), dvector![4., 4.]); + } +} diff --git a/src/curve/parameters/methods.rs b/src/curve/parameters/methods.rs new file mode 100644 index 0000000..68bacf9 --- /dev/null +++ b/src/curve/parameters/methods.rs @@ -0,0 +1,150 @@ +use crate::{ + curve::{ + parameters::Parameters, + points::{DataPoints, Points}, + }, + types::VecD, +}; + +/// Creates a parameters for every data point and distributes them equally ranging from 0 to 1. +/// This method is not recommended, as it can produce erratic shapes (such as loops) when the data is unevenly spaced. +/// see eq. (9.3) in `Piegl1997` +pub fn equally_spaced(segments: usize) -> Parameters { + let m = segments; //data_points.nrows() - 1; + let mut u_bar = VecD::zeros(m + 1); + + for g in 1..m { + u_bar[g] = g as f64 / m as f64; + } + u_bar[m] = 1f64; + + Parameters { vector: u_bar, segments } +} + +/// see eqs. (9.4) and (9.5) in `Piegl1997` +pub fn centripetal(points: &DataPoints) -> Parameters { + let m = points.segments(); //data_points.nrows() - 1; + + let mut sum = 0.0; + + for g in 1..=m { + let diff = points.get(g) - points.get(g - 1); + sum += diff.norm().sqrt() + } + + debug_assert!(sum >= f64::EPSILON.sqrt(), "sum {} is too small. Use the equally spaced methods instead", sum); + + let mut u_bar = VecD::zeros(m + 1); + + for g in 1..m { + let diff = points.get(g) - points.get(g - 1); + u_bar[g] = u_bar[g - 1] + diff.norm().sqrt() / sum; + } + + u_bar[m] = 1.0; + + Parameters { vector: u_bar, segments: points.segments() } +} + +pub fn chord_length(points: &DataPoints) -> Parameters { + let m = points.segments(); + + let mut sum = 0f64; + + for g in 1..=m { + let diff = points.get(g) - points.get(g - 1); + sum += diff.norm(); + } + + let mut u_bar = VecD::zeros(m + 1); + for g in 1..m { + let diff = points.get(g) - points.get(g - 1); + u_bar[g] = u_bar[g - 1] + diff.norm() / sum; + } + + u_bar[m] = 1f64; + + Parameters { vector: u_bar, segments: points.segments() } +} + +#[cfg(test)] +mod tests { + use nalgebra::{dmatrix, dvector}; + + use super::*; + + mod equally_spaced { + use super::*; + + #[test] + fn test() { + let points = DataPoints::new(dmatrix![1.0, 2.0, 3.0, 4.0, 5.0;]); + let knots = equally_spaced(points.segments()); + assert_eq!(knots.vector, dvector![0., 0.25, 0.5, 0.75, 1.]); + } + } + + mod chord_length { + use super::*; + + #[test] + fn linear_1() { + let points = DataPoints::new(dmatrix![1.0, 2.0, 3.0, 4.0, 5.0;]); + let knots = chord_length(&points); + assert_eq!(knots.vector, dvector![0., 0.25, 0.5, 0.75, 1.]); + } + + #[test] + fn linear_2() { + let points = DataPoints::new(dmatrix![1.0, 3.0, 5.0;]); + let knots = chord_length(&points); + assert_eq!(knots.vector, dvector![0., 0.5, 1.]); + } + + #[test] + fn non_linear_1() { + let points = DataPoints::new(dmatrix![1.0, 2.0, 5.0;]); + let knots = chord_length(&points); + assert_eq!(knots.vector, dvector![0., 0.25, 1.]); + } + + #[test] + fn non_linear_2() { + let points = DataPoints::new(dmatrix![1.0, 4.0, 5.0;]); + let knots = chord_length(&points); + assert_eq!(knots.vector, dvector![0., 0.75, 1.]); + } + } + + mod centripetal { + use super::*; + + #[test] + fn linear_1() { + let points = DataPoints::new(dmatrix![1.0, 2.0, 3.0, 4.0, 5.0;]); + let knots = centripetal(&points); + assert_eq!(knots.vector, dvector![0., 0.25, 0.5, 0.75, 1.]); + } + + #[test] + fn linear_2() { + let points = DataPoints::new(dmatrix![1.0, 3.0, 5.0;]); + let knots = centripetal(&points); + assert_eq!(knots.vector, dvector![0., 0.5, 1.]); + } + + #[test] + fn non_linear_1() { + let points = DataPoints::new(dmatrix![1.0, 2.0, 11.0;]); + let knots = centripetal(&points); + assert_eq!(knots.vector, dvector![0., 0.25, 1.]); + } + + #[test] + fn non_linear_2() { + let points = DataPoints::new(dmatrix![1.0, 10.0, 11.0;]); + let knots = centripetal(&points); + assert_eq!(knots.vector, dvector![0., 0.75, 1.]); + } + } +} diff --git a/src/curve/parameters/mod.rs b/src/curve/parameters/mod.rs new file mode 100644 index 0000000..8bae19d --- /dev/null +++ b/src/curve/parameters/mod.rs @@ -0,0 +1,43 @@ +//! Curve parameters can be generated with different methods: +//! +//! - Equally spaced parameters +//! - Centripetal method +//! - Chord-length method + +use crate::{curve::points::DataPoints, types::VecD}; + +pub mod methods; + +#[derive(Debug, Clone)] +pub struct Parameters { + vector: VecD, + segments: usize, +} + +impl Parameters { + pub fn new(vector: VecD, segments: usize) -> Self { + Parameters { vector, segments } + } + + pub fn vector(&self) -> &VecD { + &self.vector + } + + pub fn segments(&self) -> usize { + self.segments + } +} + +pub enum Method { + EquallySpaced, + Centripetal, + ChordLength, +} + +pub fn generate(points: &DataPoints, method: Method) -> Parameters { + match method { + Method::EquallySpaced => methods::equally_spaced(points.segments()), + Method::ChordLength => methods::chord_length(points), + Method::Centripetal => methods::centripetal(points), + } +} diff --git a/src/curve/points/methods/fit/fixed.rs b/src/curve/points/methods/fit/fixed.rs new file mode 100644 index 0000000..0124b5c --- /dev/null +++ b/src/curve/points/methods/fit/fixed.rs @@ -0,0 +1,227 @@ +use std::ops::SubAssign; + +use crate::{ + curve::{ + knots::Knots, + parameters::Parameters, + points::{ + methods::fit::{compute_svd, difference_operator, input_checks, FitError, Penalization}, + DataPoints, Points, + }, + }, + types::{MatD, VecD}, +}; + +pub fn fit( + knots: &Knots, + points: &DataPoints, + params: &Parameters, + penalization: Option, +) -> Result { + input_checks(knots, points, params, &penalization)?; + + let Q = generate_qvectors(knots, points, params); + let Qmat = calculate_constant_terms_matrix(knots, points, params, &Q); + let Nmat = calculate_coefficient_matrix(knots, points, params); + + let svd = compute_svd(knots, &Nmat, &penalization, Box::new(calculate_finite_difference_matrix))?; + let internal_control_points = svd.solve(&Qmat.transpose(), f64::EPSILON.sqrt()).unwrap().transpose(); + + let n = knots.segments(); + let m = points.segments(); + let mut control_points = MatD::zeros(points.dimension(), n + 1); + + // Set the first and last control point to the original data points + control_points.column_mut(0).copy_from(&points.get(0)); + control_points.column_mut(n).copy_from(&points.get(m)); + + // TODO copy matrix view instead of individual columns + for i in 1..=n - 1 { + control_points.column_mut(i).copy_from(&internal_control_points.column(i - 1)); + } + + Ok(control_points) +} + +fn generate_qvectors(knots: &Knots, points: &DataPoints, params: &Parameters) -> MatD { + let p = knots.degree(); + let n = knots.segments(); + let m = points.segments(); + let dim = points.dimension(); + + let mut Q = MatD::zeros(dim, m + 1); + + let U_bar = params.vector(); + + for g in 1..=m - 1 { + Q.column_mut(g).copy_from(&points.get(g)); + let u = U_bar[g]; + + Q.column_mut(g).sub_assign(knots.evaluate(0, 0, p, u) * points.get(0)); + Q.column_mut(g).sub_assign(knots.evaluate(0, n, p, u) * points.get(m)); + } + + Q +} + +fn calculate_constant_terms_matrix(knots: &Knots, points: &DataPoints, params: &Parameters, Q: &MatD) -> MatD { + let p = knots.degree(); + let n = knots.segments(); //settings.intended_control_points_count - 1; + let m = points.segments(); + let dim = points.dimension(); + + let U_bar = params.vector(); + + let mut Qmat = MatD::zeros(dim, n - 1); + + let mut accum = VecD::zeros(dim); + for i in 1..=n - 1 { + accum *= 0.0; + + for g in 1..=m - 1 { + let u = U_bar[g]; + accum += knots.evaluate(0, i, p, u) * Q.column(g); + } + Qmat.column_mut(i - 1).copy_from(&accum); + } + + Qmat +} + +fn calculate_coefficient_matrix(knots: &Knots, points: &DataPoints, params: &Parameters) -> MatD { + let p = knots.degree(); + let n = knots.segments(); + let m = points.segments(); + + let U_bar = params.vector(); + + let mut Nmat = MatD::zeros(m - 1, n - 1); + for g in 1..=m - 1 { + let u = U_bar[g]; + for i in 1..=n - 1 { + Nmat[(g - 1, i - 1)] = knots.evaluate(0, i, p, u); + } + } + Nmat +} + +fn calculate_finite_difference_matrix(kappa: usize, knots: &Knots) -> MatD { + let n = knots.segments(); //settings.intended_control_points_count - 1; + assert!( + kappa <= n - 2, + "The parameter kappa = {} must be greater than the number of spline segments n - 2 = {}", + kappa, + n - 2 + ); + + let mut delta_mat = MatD::zeros(n - 1 - kappa, n - 1); + + for i in 0..=n - kappa - 2 { + for j in 0..=n - 2 { + delta_mat[(i, j)] = difference_operator(i, j, kappa) as f64; + } + } + + delta_mat +} + +#[cfg(test)] +mod tests { + use approx::assert_relative_eq; + use nalgebra::{dmatrix, dvector}; + + use crate::curve::{ + knots, + knots::Method::{Averaging, Uniform}, + parameters, + parameters::Method::{ChordLength, EquallySpaced}, + points::{methods::fit::test_data_points, ControlPoints}, + Curve, + }; + + use super::*; + + #[test] + fn calculate_finite_difference_matrix_kappa1_test() { + let knots = Knots::new(1, dvector![0.0, 0.0, 0.25, 0.5, 0.75, 1.0, 1.0]); + let mat = calculate_finite_difference_matrix(1, &knots); + let expected = dmatrix![ + -1.0, 1.0, 0.0; + 0.0,-1.0, 1.0; + ]; + assert_eq!(mat, expected); + // TODO check + } + + #[test] + fn calculate_finite_difference_matrix_kappa2_test() { + let knots = Knots::new(1, dvector![0.0, 0.0, 0.25, 0.5, 0.75, 1.0, 1.0]); + let mat = calculate_finite_difference_matrix(2, &knots); + let expected = dmatrix![ + 1.0,-2.0, 1.0; + ]; + assert_eq!(mat, expected); + // TODO check + } + + // TODO reuse tests for fixed and loose + #[test] + fn unpenalized_linear() { + let points = DataPoints::new(dmatrix![ + 1., 2., 3., 4., 5.; + 1., 2., 3., 4., 5.; + ]); + + let params = parameters::generate(&points, EquallySpaced); + let knots = knots::generate(1, points.segments(), ¶ms, Uniform).unwrap(); + assert_eq!( + crate::curve::points::methods::fit::loose::fit(&knots, &points, ¶ms, None).unwrap(), + *points.matrix() + ); + } + + #[test] + fn penalized_linear() { + let points = DataPoints::new(dmatrix![ + 1., 2., 3., 4., 5.; + 1., 2., 3., 4., 5.; + ]); + + let params = parameters::generate(&points, ChordLength); + let knots = knots::generate(1, points.segments(), ¶ms, Averaging).unwrap(); + assert_relative_eq!( + fit(&knots, &points, ¶ms, Some(Penalization { lambda: 0.5, kappa: 2 })).unwrap(), + points.matrix(), + epsilon = f64::EPSILON.sqrt() + ); + } + + #[test] + fn unpenalized_nonlinear() { + let p = 1; + let data_points = test_data_points(10); + + let n = data_points.segments(); + let params = parameters::generate(&data_points, EquallySpaced); + let knots = knots::generate(p, n, ¶ms, Uniform).unwrap(); + + assert_relative_eq!( + fit(&knots, &data_points, ¶ms, None).unwrap(), + data_points.matrix(), + epsilon = f64::EPSILON.sqrt() + ); + } + + #[test] + fn penalized_nonlinear() { + let p = 2; + let data_points = test_data_points(10); + + let params = parameters::generate(&data_points, EquallySpaced); + let knots = knots::generate(p, data_points.segments(), ¶ms, Uniform).unwrap(); + let points = fit(&knots, &data_points, ¶ms, Some(Penalization { lambda: 1.0, kappa: 2 })).unwrap(); + let c = Curve::new(knots, ControlPoints::new(points)).unwrap(); + + assert_relative_eq!(c.evaluate(0.5).unwrap(), dvector![0.0, 0.0], epsilon = f64::EPSILON.sqrt()); + } +} diff --git a/src/curve/points/methods/fit/loose.rs b/src/curve/points/methods/fit/loose.rs new file mode 100644 index 0000000..2194e91 --- /dev/null +++ b/src/curve/points/methods/fit/loose.rs @@ -0,0 +1,167 @@ +use crate::{ + curve::{ + knots::Knots, + parameters::Parameters, + points::{ + methods::fit::{compute_svd, difference_operator, input_checks, FitError, Penalization}, + DataPoints, Points, + }, + }, + types::MatD, +}; + +pub fn fit( + knots: &Knots, + points: &DataPoints, + params: &Parameters, + penalization: Option, +) -> Result { + input_checks(knots, points, params, &penalization)?; + + let Nmat = calculate_coefficient_matrix(knots, points, params); + + let svd = compute_svd(knots, &Nmat, &penalization, Box::new(calculate_finite_difference_matrix))?; + let mat = Nmat.transpose() * points.matrix().transpose(); + let control_points = svd.solve(&mat, f64::EPSILON.sqrt()).unwrap().transpose(); + + Ok(control_points) +} + +fn calculate_coefficient_matrix(knots: &Knots, points: &DataPoints, params: &Parameters) -> MatD { + let p = knots.degree(); + let n = knots.segments(); + let m = points.segments(); + + let U_bar = params.vector(); + + let mut Nmat = MatD::zeros(m + 1, n + 1); + for g in 0..=m { + let u = U_bar[g]; + for i in 0..=n { + Nmat[(g, i)] = knots.evaluate(0, i, p, u); + } + } + Nmat +} + +fn calculate_finite_difference_matrix(kappa: usize, knots: &Knots) -> MatD { + let n = knots.segments(); + assert!(kappa <= n, "The parameter kappa = {} must be greater than the number of spline segments n = {}", kappa, n); + + let mut delta_mat = MatD::zeros(n + 1 - kappa, n + 1); + + for i in 0..=n - kappa { + for j in 0..=n { + delta_mat[(i, j)] = difference_operator(i, j, kappa) as f64; + } + } + + delta_mat +} + +#[cfg(test)] +mod tests { + use approx::assert_relative_eq; + use nalgebra::{dmatrix, dvector}; + + use crate::curve::{ + knots, + knots::Method::{Averaging, Uniform}, + parameters, + parameters::Method::{ChordLength, EquallySpaced}, + points::{methods::fit::test_data_points, ControlPoints}, + Curve, + }; + + use super::*; + + #[test] + fn calculate_finite_difference_matrix_kappa1_test() { + let knots = Knots::new(1, dvector![0.0, 0.0, 0.25, 0.5, 0.75, 1.0, 1.0]); + let mat = calculate_finite_difference_matrix(1, &knots); + let expected = dmatrix![ + -1.0, 1.0, 0.0, 0.0, 0.0; + 0.0,-1.0, 1.0, 0.0, 0.0; + 0.0, 0.0,-1.0, 1.0, 0.0; + 0.0, 0.0, 0.0,-1.0, 1.0; + ]; + assert_eq!(mat, expected); + // TODO check + } + + #[test] + fn calculate_finite_difference_matrix_kappa2_test() { + let knots = Knots::new(1, dvector![0.0, 0.0, 0.25, 0.5, 0.75, 1.0, 1.0]); + let mat = calculate_finite_difference_matrix(2, &knots); + let expected = dmatrix![ + 1.0,-2.0, 1.0, 0.0, 0.0; + 0.0, 1.0,-2.0, 1.0, 0.0; + 0.0, 0.0, 1.0,-2.0, 1.0; + ]; + assert_eq!(mat, expected); + // TODO check + } + + #[test] + fn unpenalized_linear() { + let points = DataPoints::new(dmatrix![ + 1., 2., 3., 4., 5.; + 1., 2., 3., 4., 5.; + ]); + + let params = parameters::generate(&points, EquallySpaced); + let knots = knots::generate(1, points.segments(), ¶ms, Uniform).unwrap(); + assert_eq!(fit(&knots, &points, ¶ms, None).unwrap(), *points.matrix()); + } + + #[test] + fn penalized_linear() { + let points = DataPoints::new(dmatrix![ + 1., 2., 3., 4., 5.; + 1., 2., 3., 4., 5.; + ]); + + let params = parameters::generate(&points, ChordLength); + let knots = knots::generate(1, points.segments(), ¶ms, Averaging).unwrap(); + assert_relative_eq!( + fit(&knots, &points, ¶ms, Some(Penalization { lambda: 0.5, kappa: 2 })).unwrap(), + points.matrix(), + epsilon = f64::EPSILON.sqrt() + ); + } + + #[test] + fn unpenalized_nonlinear() { + let p = 1; + let data_points = test_data_points(10); + + let n = data_points.segments(); + let params = parameters::generate(&data_points, EquallySpaced); + let knots = knots::generate(p, n, ¶ms, Uniform).unwrap(); + + assert_relative_eq!( + fit(&knots, &data_points, ¶ms, None).unwrap(), + data_points.matrix(), + epsilon = f64::EPSILON.sqrt() + ); + } + + #[test] + fn penalized_nonlinear() { + let p = 2; + let data_points = test_data_points(10); + + let params = parameters::generate(&data_points, EquallySpaced); + let knots = knots::generate(p, data_points.segments(), ¶ms, Uniform).unwrap(); + let points = crate::curve::points::methods::fit::fixed::fit( + &knots, + &data_points, + ¶ms, + Some(Penalization { lambda: 1.0, kappa: 2 }), + ) + .unwrap(); + let c = Curve::new(knots, ControlPoints::new(points)).unwrap(); + + assert_relative_eq!(c.evaluate(0.5).unwrap(), dvector![0.0, 0.0], epsilon = f64::EPSILON.sqrt()); + } +} diff --git a/src/curve/points/methods/fit/mod.rs b/src/curve/points/methods/fit/mod.rs new file mode 100644 index 0000000..ab60832 --- /dev/null +++ b/src/curve/points/methods/fit/mod.rs @@ -0,0 +1,130 @@ +use nalgebra::{DMatrix, Dyn, SVD}; +use thiserror::Error; + +use crate::{ + curve::{ + knots::{is_uniform, Knots}, + parameters::Parameters, + points::DataPoints, + }, + types::MatD, +}; + +pub mod fixed; +pub mod loose; + +pub enum Method { + FixedEnds, + LooseEnds, +} + +#[derive(Error, Debug, PartialEq)] +pub enum FitError { + #[error("The penalization parameter `lambda = {lambda}` cannot be negative.")] + NegativeLambda { lambda: f64 }, + + #[error( + "The number of data point segments m = {m} must be greater than the requested polynomial segments n = {n}" + )] + RequestedPolynomialSegmentAndDataSegementMismatch { n: usize, m: usize }, + + #[error( + "The requested number of polynomial segments n = {n} must be equal or greater than the spline degree p = {p}" + )] + RequestedPolynomialSegmentAndSplineDegreeMismatch { n: usize, p: usize }, + + #[error("The requested number of polynomial segments n = {n} must be larger than penalization offset kappa kappa = {kappa}")] + RequestedPolynomialSegmentAndPenalizationKappaMismatch { n: usize, kappa: usize }, + + #[error("The number of data point segments m = {m} must be equal to the number of parameter segments mp = {mp}.")] + DataSegmentsAndParameterSegmentsMismatch { m: usize, mp: usize }, + + #[error("Penalization requires equidistant/uniform knots.")] // See `Eilers1996`. + NonUniformKnots, +} + +pub struct Penalization { + pub lambda: f64, + pub kappa: usize, + // TODO add `new` method and assertions according to below + // A B-spline curve C(u) of degree p can be generated via least-squares minimization and results in + // an approximation of the (m + 1) data points with dimension N . The number of control points (n + 1) + // can be specified but must be smaller then the number of data points and greater than the spline degree (m > n ≥ + // p). +} + +fn input_checks( + knots: &Knots, + points: &DataPoints, + params: &Parameters, + penalization: &Option, +) -> Result<(), FitError> { + match (knots.segments(), points.segments(), params.segments(), knots.degree(), penalization) { + (n, m, _, _, _) if n > m => Err(FitError::RequestedPolynomialSegmentAndDataSegementMismatch { n, m }), + (_, m, mp, _, _) if m != mp => Err(FitError::DataSegmentsAndParameterSegmentsMismatch { m, mp }), + (n, _, _, p, _) if n < p => Err(FitError::RequestedPolynomialSegmentAndSplineDegreeMismatch { n, p }), + (n, _, _, _, Some(pen)) if n - 1 < pen.kappa => { + Err(FitError::RequestedPolynomialSegmentAndPenalizationKappaMismatch { n, kappa: pen.kappa }) + } + _ => Ok(()), + } +} + +pub fn compute_svd( + knots: &Knots, + Nmat: &MatD, + penalization: &Option, + calculate_finite_difference_matrix: Box MatD>, +) -> Result, FitError> { + let mut mat = Nmat.transpose() * Nmat; + + match penalization { + Some(penalization) => match penalization.lambda { + l if l < 0.0 => return Err(FitError::NegativeLambda { lambda: l }), + l if l > 0.0 && !is_uniform(knots).unwrap() /*TODO refactor errors and eliminate unwrap*/ => return Err(FitError::NonUniformKnots), + l => { + let delta_mat = calculate_finite_difference_matrix(penalization.kappa, knots); + mat += l * (delta_mat.transpose() * delta_mat); + } + }, + None => {} + } + + Ok(SVD::new(mat, true, true)) +} + +// finite difference matrix penalizing BSplines +fn difference_operator(i: usize, j: usize, kappa: usize) -> isize { + match kappa { + 1 => { + if i == j { + return -1; + } + if i + 1 == j { + return 1; + } + 0 + } + k if k > 1 => difference_operator(i + 1, j, k - 1) - difference_operator(i, j, k - 1), + _ => 0, + } +} + +// TODO remove and use explicit input +fn test_data_points(npoints: usize) -> DataPoints { + let inc = 5.0 / npoints as f64; + let shift = (npoints - 1) as f64 * inc / 2.0; + DataPoints::new(DMatrix::from_fn(2, npoints, |r, c| { + if r == 0 { + // x coords + c as f64 * inc - shift + } else { + // y coord + if c % 2 == 0 { + 0.5 + } else { + -0.5 + } + } + })) +} diff --git a/src/curve/points/methods/interpolation.rs b/src/curve/points/methods/interpolation.rs new file mode 100644 index 0000000..a864bc0 --- /dev/null +++ b/src/curve/points/methods/interpolation.rs @@ -0,0 +1,71 @@ +use nalgebra::SVD; + +use crate::{ + curve::{ + knots::Knots, + parameters::Parameters, + points::{DataPoints, Points}, + }, + types::MatD, +}; + +pub fn interpolate(knots: &Knots, points: &DataPoints, params: &Parameters) -> MatD { + let p = knots.degree(); + let m = points.segments(); + let n = m; + + let Ubar = params.vector(); + + let mut Nmat = MatD::zeros(points.count(), points.count()); + for i in 0..=n { + for g in 0..=m { + Nmat[(g, i)] = knots.evaluate(0, i, p, Ubar[g]); + } + } + + let svd = SVD::new(Nmat, true, true); + let mat = points.matrix().transpose(); + svd.solve(&mat, f64::EPSILON.sqrt()).unwrap().transpose() +} + +#[cfg(test)] +mod tests { + use approx::assert_relative_eq; + use nalgebra::dmatrix; + + use crate::curve::{knots, knots::Method::Averaging, parameters, parameters::Method::ChordLength}; + + use super::*; + + #[test] + fn linear() { + let points = DataPoints::new(dmatrix![ + 1., 2., 3., 4.; + 1., 2., 3., 4.; + ]); + + let params = parameters::generate(&points, ChordLength); + let knots = knots::generate(1, points.segments(), ¶ms, Averaging).unwrap(); + + assert_eq!(interpolate(&knots, &points, ¶ms), points.matrix); + } + + #[test] + fn quadratic() { + let points = DataPoints::new(dmatrix![ + 1., 2., 3., 4.; + 1., 2., 3., 4.; + ]); + let params = parameters::generate(&points, ChordLength); + let knots = knots::generate(2, points.segments(), ¶ms, Averaging).unwrap(); + + assert_relative_eq!( + interpolate(&knots, &points, ¶ms), + dmatrix![ + 1., 1.75, 3.25, 4.; + 1., 1.75, 3.25, 4.; + ], + epsilon = f64::EPSILON.sqrt() + ); + } +} diff --git a/src/curve/points/methods/mod.rs b/src/curve/points/methods/mod.rs new file mode 100644 index 0000000..2475cbd --- /dev/null +++ b/src/curve/points/methods/mod.rs @@ -0,0 +1,2 @@ +pub mod fit; +pub mod interpolation; diff --git a/src/curve/points/mod.rs b/src/curve/points/mod.rs new file mode 100644 index 0000000..1995f67 --- /dev/null +++ b/src/curve/points/mod.rs @@ -0,0 +1,222 @@ +use std::ops::MulAssign; + +use crate::{ + curve::knots::Knots, + types::{ControlPointDerivatves, MatD, VecD, VecDView, VecDViewMut}, +}; + +pub mod methods; + +#[derive(PartialEq, Debug, Clone)] +pub struct ControlPoints { + pub(crate) Pk: ControlPointDerivatves, + k_max: usize, +} + +#[derive(PartialEq, Debug, Clone)] +pub struct DataPoints { + pub(self) matrix: MatD, +} + +pub trait Points { + fn matrix(&self) -> &MatD; + fn matrix_mut(&mut self) -> &mut MatD; + + fn get(&self, i: usize) -> VecDView { + self.matrix().column(i) + } + + fn get_mut(&mut self, i: usize) -> VecDViewMut { + self.matrix_mut().column_mut(i) + } + + fn dimension(&self) -> usize { + self.matrix().nrows() + } + + fn count(&self) -> usize { + self.matrix().ncols() + } + + fn is_empty(&self) -> bool { + self.matrix().is_empty() + } +} + +impl Points for DataPoints { + fn matrix(&self) -> &MatD { + &self.matrix + } + + fn matrix_mut(&mut self) -> &mut MatD { + &mut self.matrix + } +} + +impl DataPoints { + pub fn new(matrix: MatD) -> Self { + DataPoints { matrix } + } + + pub fn reverse(&mut self) -> &mut Self { + reverse(self.matrix_mut()); + self + } + pub fn segments(&self) -> usize { + self.count() - 1 + } +} + +impl Points for ControlPoints { + fn matrix(&self) -> &MatD { + &self.Pk[0] + } + + fn matrix_mut(&mut self) -> &mut MatD { + &mut self.Pk[0] + } +} + +impl ControlPoints { + pub fn new(points: MatD) -> Self { + ControlPoints { Pk: vec![points], k_max: 0 } + } + + pub fn new_with_capacity(points: MatD, capacity: usize) -> ControlPoints { + let mut Pk: Vec = Vec::with_capacity(capacity); + Pk.push(points); + + ControlPoints { Pk, k_max: 0 } + } + + pub fn segments(&self) -> usize { + self.count() - 1 + } + + pub fn matrix_derivative(&self, derivative: usize) -> &MatD { + assert!(derivative <= self.k_max, "Derivative {} is not calculated", derivative); + &self.Pk[derivative] + } + + pub fn matrix_derivative_mut(&mut self, derivative: usize) -> &mut MatD { + assert!(derivative <= self.k_max, "Derivative {} is not calculated", derivative); + &mut self.Pk[derivative] + } + + pub fn count(&self) -> usize { + self.count_derivative(0) + } + pub fn count_derivative(&self, k: usize) -> usize { + self.Pk[k].ncols() + } + + pub fn max_derivative(&self) -> usize { + self.k_max + } + + pub fn derive(&mut self, knots: &Knots) { + let p = knots.degree(); + let n = self.segments(); + + self.Pk.truncate(1); + for k in 1..=p { + let mut Pnew = MatD::zeros(self.dimension(), n - k + 1); + // TODO iter over points instead + for (i, mut col) in Pnew.column_iter_mut().enumerate() { + col.copy_from(&self.derive_single_point(i, k, knots)); + } + self.Pk.push(Pnew); + } + self.k_max = p; + } + + fn derive_single_point(&self, i: usize, k_max: usize, knots: &Knots) -> VecD { + let p = knots.degree(); + + if k_max == 0 { + return self.Pk[0].column(i).clone_owned(); + } + + let U0 = knots.vector(); + if U0[i + p + 1] == U0[i + k_max] { + 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)) + } + + pub fn reverse(&mut self) -> &mut Self { + for k in 0..=self.k_max { + let Pk = self.matrix_derivative_mut(k); + reverse(Pk); + + // Odd derivatives change their sign upon reversal + if k % 2 == 1 { + Pk.mul_assign(-1.0); + } + } + self + } +} + +pub(crate) fn reverse(points: &mut MatD) { + let ncols = points.ncols(); + let half_ncols = points.ncols() / 2; + + for i in 0..half_ncols { + points.swap_columns(i, ncols - 1 - i); + } +} + +pub(crate) fn reversed(points: &MatD) -> MatD { + let mut copy = points.clone(); + reverse(&mut copy); + copy +} + +#[cfg(test)] +mod tests { + use nalgebra::dmatrix; + + use super::*; + + fn control_points_example() -> ControlPoints { + ControlPoints::new(dmatrix![ + 1., 3., 5., 7.; + 2., 4., 6., 8.; + ]) + } + + #[test] + fn dimension() { + assert_eq!(control_points_example().dimension(), 2); + } + + #[test] + fn count() { + assert_eq!(control_points_example().count(), 4); + } + + #[test] + fn segments() { + assert_eq!(control_points_example().segments(), 3); + } + + #[test] + fn reverse() { + assert_eq!( + control_points_example().reverse().Pk, + vec![dmatrix![ + 7., 5., 3., 1.; + 8., 6., 4., 2.; + ]] + ); + } + + /* TODO + #[test] + fn reverse_derivatives() { + todo!() + }*/ +} diff --git a/src/docs-header.html b/src/docs-header.html new file mode 100644 index 0000000..8df31a0 --- /dev/null +++ b/src/docs-header.html @@ -0,0 +1,83 @@ + + + + diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..6cfe901 --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,62 @@ +//#![warn(missing_docs)] +//#![warn(missing_doc_code_examples)] +#![cfg_attr(feature = "doc-images", +cfg_attr(all(), +doc = ::embed_doc_image::embed_image!("split-after", "doc-images/plots/manipulation/split-after.svg")))] +//! `bsplines` is a Rust library for vectorized, N-dimensional B-spline curves and their derivatives based on +//! [nalgebra]. +//! +//! ## Features +//! - Create `$x$`-dimensional (`$x = 1, 2, 3,\dots, N$`) curves of arbitrary polynomial degree `$p$`. +//! - Efficient curve evaluation for all available derivatives `$k = 0, 1, \dots, p$`. +//! - Built with [nalgebra](https://crates.io/crates/nalgebra) to store point data in contiguous arrays +//! - Multiple methods for +//! - [curve generation][curve::generation] +//! - [curve parametrization][curve::parameters] +//! - [knot generation][curve::knots] +//! - [curve manipulation][manipulation] +//! - [knot insertion][manipulation::insert] +//! - [reversing][manipulation::reverse] +//! - [splitting][manipulation::split] +//! - [merging][manipulation::split] +//! +//! ## Mathematical Definition +//! +//! 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 +//! differentiable. The overall functions are parametrized over finite domains with the co-domain being an +//! `$N$`-dimensional vector space. They can describe curves, but also surfaces. +//! +//! A B-spline curve can be defined by +//! ```math +//! \mathcal{C}(u) = \sum_{i=0}^{n} \mathcal{N}_{i,p}^{\boldsymbol{U}} (u)\, \boldsymbol{P}_i +//! ``` +//! with +//! `$u$`, a parameter on the curve, usually limited to `$u\in \left[0,1\right]$`. +//! `$n$` the number of polynomial spline segments +//! `$p$`, the spline degree +//! `$\boldsymbol{U}$`, the knot vector +//! `$\mathcal{N}$`, the `$n+1$` spline basis functions of degree `p` +//! `$\boldsymbol{P}$`, the `$n+1$` control points that can be of arbitrary dimension +//! +//! These characteristics lead to many desirable properties. +//! The piecewise definition makes B-spline functions versatile allowing to interpolate or approximate +//! complex-shaped and high-dimensional data, while maintaining a low polynomial degree. Because of the polynomial +//! nature, all possible derivatives are accessible. +//! +//! ![A 2D B-Spline curve.][split-after] +//! +//! Still, evaluations or spatial manipulations can be executed fast because only local polynomial segments must be +//! considered and the associated numerical procedures are stable. Lastly, polynomials represent a memory-efficient way +//! of storing spatial information as few polynomial coefficients suffice to describe complex shapes. + +//! ## Literature: +//! | | | +//! |-----------:|:-------------------------------------------------------------------------------------------------------------------------------------------------------------------| +//! | Piegl1997 | Piegl, L., Tiller, W. The NURBS Book. Monographs in Visual Communication. Springer, Berlin, Heidelberg, 2nd ed., 1997. | +//! | Eilers1996 | Eilers, P. H. C., Marx, B. D., Flexible smoothing with B -splines and penalties, Stat. Sci., 11(2) (1996) 89–121. | +//! | Tai2003 | Tai, C.-L., Hu, S.-M., Huang, Q.-X., Approximate merging of B-spline curves via knot adjustment and constrained optimization, Comput. Des., 35(10) (2003) 893–899. | + +pub mod curve; +pub mod manipulation; +pub mod types; diff --git a/src/manipulation/insert.rs b/src/manipulation/insert.rs new file mode 100644 index 0000000..d38d150 --- /dev/null +++ b/src/manipulation/insert.rs @@ -0,0 +1,170 @@ +#![cfg_attr(feature = "doc-images", +cfg_attr(all(), +doc = ::embed_doc_image::embed_image!("insert-before", "doc-images/plots/manipulation/insert-before.svg"), +doc = ::embed_doc_image::embed_image!("insert-after", "doc-images/plots/manipulation/insert-after.svg")))] +//! Inserts an additional knot into the curve. +//! +//! | A curve before knot insertion. | The curve after knot insertion at `$u=4/5$`. | +//! |:------------------------------:|:--------------------------------------------:| +//! | ![][insert-before] | ![][insert-after] | + +use std::ops::AddAssign; + +use thiserror::Error; + +use crate::{ + curve::{knots::DomainKnotComparatorType, points::Points, Curve}, + types::MatD, +}; + +#[derive(Error, Debug, PartialEq)] +pub enum InsertError { + #[error("Parameter `u = {u}` lies outside the interval `({lower_bound}, {upper_bound})`.")] + OutOfBounds { u: f64, lower_bound: f64, upper_bound: f64 }, + + #[error( + "The knot `u = {u}` has a multiplicity of `m = {m}` already. \ + Therefore, the knot cannot be inserted as this would exceed the maximum \ + multiplicity corresponding to the curve degree with `p = {p}`." + )] + MultiplicityError { u: f64, m: usize, p: usize }, +} + +/// Knot insertion algorithm by Boehm +/// `u` the knot to be inserted. The value must be in `$u\in(0, 1)$` +pub fn insert(c: &mut Curve, u: f64) -> Result<(), InsertError> { + if u <= 0.0 || u >= 1.0 { + return Err(InsertError::OutOfBounds { u, lower_bound: 0.0, upper_bound: 1.0 }); + } + + let p = c.degree(); + + let m = c.knots.multiplicity(u); + if c.knots.multiplicity(u) > p { + return Err(InsertError::MultiplicityError { u, m, p }); + } + + let dim = c.points.dimension(); + + let U_old = c.knots.vector(); + let P_old = c.points.matrix(); + + let l = c.knots.find_idx(u, 0, DomainKnotComparatorType::LeftOrEqual); + + // Insert u into the knot vector + let U_new = U_old.clone().insert_row(l + 1, u); + + // Compute the new control points. + // Only the control points `l-p+1` to `l` change. + let control_point_count = c.points.count(); + + let mut P_new = MatD::zeros(dim, control_point_count + 1); + + let top_cols = l - p + 1; + P_new.columns_mut(0, top_cols).copy_from(&P_old.columns(0, top_cols)); + + let bot_cols = control_point_count - l; + P_new.columns_mut(P_new.ncols() - bot_cols, bot_cols).copy_from(&P_old.columns(P_old.ncols() - bot_cols, bot_cols)); + + let mut alpha: f64; + for i in (l - p + 1)..=l { + alpha = (u - U_old[i]) / (U_old[i + p] - U_old[i]); + + P_new.column_mut(i).add_assign((1. - alpha) * P_old.column(i - 1) + alpha * P_old.column(i)); + } + + c.knots.Uk[0] = U_new; + c.points.Pk[0] = P_new; + c.calculate_derivatives(); + Ok(()) +} + +#[cfg(test)] +mod tests { + use nalgebra::{dmatrix, dvector}; + + use crate::curve::{ + generation::{generate, Generation::Manual}, + knots::Generation::Uniform, + points::ControlPoints, + }; + + use super::*; + + #[test] + fn degree_1() { + let mut c = + generate(Manual { degree: 1, points: ControlPoints::new(dmatrix![-1., 1.;]), knots: Uniform }).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 1., 1.]); + + insert(&mut c, 0.5).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0.5, 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1., 0., 1.;]); + } + + #[test] + fn degree_2() { + let mut c = + generate(Manual { degree: 2, points: ControlPoints::new(dmatrix![-1., 0., 1.;]), knots: Uniform }).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 1., 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1., 0., 1.;]); + + insert(&mut c, 0.5).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 0.5, 1., 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1., -0.5, 0.5, 1.;]); + } + + #[test] + fn degree_2_preexisting_knot() { + let mut c = + generate(Manual { degree: 2, points: ControlPoints::new(dmatrix![-1.5, -0.5, 0.5, 1.5;]), knots: Uniform }) + .unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 0.5, 1., 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1.5, -0.5, 0.5, 1.5;]); + + insert(&mut c, 0.5).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 0.5, 0.5, 1., 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1.5, -0.5, 0.0, 0.5, 1.5;]); + } + + #[test] + fn degree_1_repeated_knot() { + let mut c = + generate(Manual { degree: 1, points: ControlPoints::new(dmatrix![-1., 1.;]), knots: Uniform }).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 1., 1.]); + + insert(&mut c, 0.5).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0.5, 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1., 0., 1.;]); + + insert(&mut c, 0.5).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0.5, 0.5, 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1., 0., 0., 1.;]); + } + + #[test] + fn degree_2_repeated_knots() { + let mut c = + generate(Manual { degree: 2, points: ControlPoints::new(dmatrix![-1., 0., 1.;]), knots: Uniform }).unwrap(); + let u = 0.5; + let expectedEvaluationResult = dvector![0.0]; + + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 1., 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1., 0., 1.;]); + //assert_eq!(c.evaluate(u).unwrap(), expectedEvaluationResult); + + insert(&mut c, u).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., u, 1., 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1., -0.5, 0.5, 1.;]); + //assert_eq!(c.evaluate(u).unwrap(), expectedEvaluationResult); + + insert(&mut c, u).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., u, u, 1., 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1., -0.5, 0.0, 0.5, 1.;]); + assert_eq!(c.evaluate(u).unwrap(), expectedEvaluationResult); + + insert(&mut c, u).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., u, u, u, 1., 1., 1.]); + assert_eq!(c.points.matrix(), &dmatrix![-1., -0.5, 0.0, 0.0, 0.5, 1.;]); + } +} diff --git a/src/manipulation/merge.rs b/src/manipulation/merge.rs new file mode 100644 index 0000000..3a6979d --- /dev/null +++ b/src/manipulation/merge.rs @@ -0,0 +1,691 @@ +#![cfg_attr(feature = "doc-images", +cfg_attr(all(), +doc = ::embed_doc_image::embed_image!("merge-before", "doc-images/plots/manipulation/merge-before.svg"), +doc = ::embed_doc_image::embed_image!("merge-after", "doc-images/plots/manipulation/merge-after.svg"), +doc = ::embed_doc_image::embed_image!("merge-after-left-end-constrained", "doc-images/plots/manipulation/merge-after-left-end-constrained.svg"), +doc = ::embed_doc_image::embed_image!("merge-after-right-start-constrained", "doc-images/plots/manipulation/merge-after-right-start-constrained.svg")))] +//! Combine two independent curves into one. +//! +//! | Two curves before merging. | The resulting curve after merging. | +//! |:--------------------------------------------:|:-----------------------------------------------:| +//! | ![][merge-before] | ![][merge-after] | +//! +//! Constraints can be set so that certain points on the curve stay fixed. +//! +//! | After merging with the left end constrained. | After merging with the right start constrained. | +//! |:--------------------------------------------:|:-----------------------------------------------:| +//! | ![][merge-after-left-end-constrained] | ![][merge-after-right-start-constrained] | + +use std::ops::{AddAssign, DivAssign, SubAssign}; + +use nalgebra::SVD; +use thiserror::Error; + +use crate::{ + curve, + curve::{ + knots::{evaluate, is_clamped, is_normed, reversed, Knots}, + points::{ControlPoints, Points}, + Curve, CurveError, + }, + manipulation::merge::MergeError::CurveGenerationFailure, + types::{MatD, VecD, VecHelpers}, +}; + +pub struct Constraints { + pub params: Vec, +} + +impl Constraints { + pub fn count(&self) -> usize { + self.params.len() + } + + pub fn segments(&self) -> usize { + self.count() - 1 + } +} + +pub struct ConstrainedCurve<'a> { + pub(crate) curve: &'a Curve, + pub(crate) constraints: Constraints, +} + +#[derive(Error, Debug, PartialEq)] +pub enum MergeError { + #[error("The degree of the left curve `p = {left}` differs from the right curve `p = {right}`.")] + DegreeMismatch { left: usize, right: usize }, + + #[error("The dimension of the left curve `dim = {left}` differs from the right curve `dim = {right}`.")] + DimensionMismatch { left: usize, right: usize }, + + #[error("Curves must be clamped")] + UnclampedCurve, + + #[error("Curves must be normed.")] + UnnormedCurve, + + #[error( + "The total number of constrained points `{totalConstraints}` must be \ + less than the curve degree `p = {degree}`. \ + Otherwise no solution for the linear system of equations exists." + )] + TooManyConstraints { totalConstraints: usize, degree: usize }, + + #[error("Curve generation failed with error {err}.")] + CurveGenerationFailure { err: CurveError }, +} + +// Keeps the start of spline 1 fixed +pub fn merge_from(a: &Curve, b: &Curve) -> Result { + merge_with_constraints( + &ConstrainedCurve { curve: a, constraints: Constraints { params: vec![1.] } }, + &ConstrainedCurve { curve: b, constraints: Constraints { params: vec![] } }, + ) +} + +// Keeps the end of spline 2 fixed +pub fn merge_to(a: &Curve, b: &Curve) -> Result { + merge_with_constraints( + &ConstrainedCurve { curve: a, constraints: Constraints { params: vec![] } }, + &ConstrainedCurve { curve: b, constraints: Constraints { params: vec![0.] } }, + ) +} + +pub fn merge(a: &Curve, b: &Curve) -> Result { + merge_with_constraints( + &ConstrainedCurve { curve: a, constraints: Constraints { params: vec![] } }, + &ConstrainedCurve { curve: b, constraints: Constraints { params: vec![] } }, + ) +} + +pub(crate) fn merge_with_constraints(a: &ConstrainedCurve, b: &ConstrainedCurve) -> Result { + let p_a = a.curve.degree(); + let p_b = b.curve.degree(); + + if p_a != p_b { + return Err(MergeError::DegreeMismatch { left: p_a, right: p_b }); + } + + if a.curve.dimension() != b.curve.dimension() { + return Err(MergeError::DimensionMismatch { left: a.curve.dimension(), right: b.curve.dimension() }); + } + + if !is_clamped(&a.curve.knots) || !is_clamped(&b.curve.knots) { + return Err(MergeError::UnclampedCurve); + } + + if !is_normed(&a.curve.knots) || !is_normed(&b.curve.knots) { + return Err(MergeError::UnnormedCurve); + } + + let totalConstraints = a.constraints.count() + b.constraints.count(); + if totalConstraints >= p_a { + return Err(MergeError::TooManyConstraints { totalConstraints, degree: p_a }); + } + + let shifts = solve_linear_equation_system(a, b); + let (Sshifted, Tshifted) = generate_shifted_control_point_vectors_of_spline1and2(a.curve, b.curve, &shifts); + + // adjust and merge knot vectors + let (Vadj, Wrev, Wadj) = adjust_knots_of_both_splines(a.curve, b.curve); + let U_merged = create_merged_knot_vector(a.curve, b.curve, &Vadj, &Wadj); + + // adjust control points + let (Sadj, Tadj) = + adjust_shifted_control_points_of_both_splines(a.curve, &Sshifted, &Tshifted, &Vadj, &Wrev, &Wadj); + + // do actual merge + let P_merged = generate_control_point_vector_of_merged_spline(a.curve, b.curve, &Sadj, &Tadj); + + Curve::new(Knots::new(p_a, U_merged), ControlPoints::new(P_merged)).map_err(|err| CurveGenerationFailure { err }) +} + +fn construct_Nmat(a: &ConstrainedCurve, b: &ConstrainedCurve) -> MatD { + let p = a.curve.degree(); + + let n_constraints_a = a.constraints.count(); + let n_constraints_b = b.constraints.count(); + + let nmat_dim = 3 * p + n_constraints_a + n_constraints_b; + let mut Nmat = MatD::zeros(nmat_dim, nmat_dim); + + Nmat.view_mut((0, 0), (2 * p, 2 * p)).copy_from(&MatD::identity(2 * p, 2 * p)); + + Nmat.view_mut((2 * p, 0), (p, p)).copy_from(&calculateKV(a.curve)); + Nmat.view_mut((2 * p, p), (p, p)).copy_from(&calculateKW(b.curve)); + + Nmat.view_mut((0, 2 * p), (p, p)).copy_from(&calculateIV(a.curve)); + Nmat.view_mut((p, 2 * p), (p, p)).copy_from(&calculateJW(b.curve)); + + // calculate block matrices for the constraints + if n_constraints_a > 0 { + Nmat.view_mut((3 * p, 0), (n_constraints_a, p)).copy_from(&calculateGV(a)); + Nmat.view_mut((0, 3 * p), (p, n_constraints_a)).copy_from(&calculateIpV(a)); + } + if n_constraints_b > 0 { + Nmat.view_mut((3 * p + n_constraints_a, p), (n_constraints_b, p)).copy_from(&calculateHW(b)); + Nmat.view_mut((p, 3 * p + n_constraints_a), (p, n_constraints_b)).copy_from(&calculateJppW(b)); + } + + Nmat +} + +fn calculateKV(a: &Curve) -> MatD { + let p = a.degree(); + let m = a.segments(); + + let Vk = &a.knots.Uk; + let S = a.points.matrix(); + + let mut KV = MatD::zeros(p, p); + + for k in 0..=p - 1 { + for i in m - p + 1..=m { + 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); + } + KV[(k, i - (m + 1 - p))] = sum; + } + } + KV +} + +fn calculateKW(b: &Curve) -> MatD { + let p = b.degree(); + let o = b.segments(); + let Wk = &b.knots.Uk; + let T = b.points.matrix(); + + let mut KW = MatD::zeros(p, p); + + for k in 0..=p - 1 { + for j in 0..=p - 1 { + 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); + } + KW[(k, j)] = sum; + } + } + KW *= -1.0; + + KW +} + +fn calculateIV(a: &Curve) -> MatD { + let p = a.degree(); + let m = a.segments(); + let Vk = &a.knots.Uk; + let S = a.points.matrix(); + + let mut IV = MatD::zeros(p, p); + + for i in m - p + 1..=m { + for k in 0..=p - 1 { + 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); + } + IV[(i - (m + 1 - p), k)] = sum; + } + } + IV *= 0.5; + + IV +} + +fn calculateJW(b: &Curve) -> MatD { + let p = b.degree(); + let o = b.segments(); + let Wk = &b.knots.Uk; + let T = b.points.matrix(); + + let mut JW = MatD::zeros(p, p); + + for j in 0..=p - 1 { + for k in 0..=p - 1 { + 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); + } + JW[(j, k)] = sum; + } + } + JW *= -0.5; + + JW +} + +fn calculateGV(a: &ConstrainedCurve) -> MatD { + let p = a.curve.degree(); + let m = a.curve.segments(); + let Vk0 = a.curve.knots.vector(); + + let mg = a.constraints.segments(); + + let mut GV = MatD::zeros(mg + 1, p); + + 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 +} + +fn calculateHW(b: &ConstrainedCurve) -> MatD { + let p = b.curve.degree(); + let o = b.curve.segments(); + let Wk0 = b.curve.knots.vector(); + + let oh = b.constraints.segments(); + let mut HW = MatD::zeros(oh + 1, p); + + for h in 0..=oh { + for i in 0..=p - 1 { + HW[(h, i)] = evaluate(i, p, o, Wk0, b.constraints.params[h], 0); + } + } + + // TODO use below + /*let mut HW = MatD::zeros(constraints_b.count(), p); + + for (h, u) in constraints_b.params.iter().enumerate() { + for i in 0..=p - 1 { + HW[(h, i)] = evaluate(i, p, o, Wk0, *u); + } + }*/ + HW +} + +fn calculateIpV(a: &ConstrainedCurve) -> MatD { + let p = a.curve.degree(); + let m = a.curve.segments(); + let Vk0 = a.curve.knots.vector(); + + let mg = a.constraints.segments(); + + let mut IpV = MatD::zeros(p, mg + 1); + + 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 *= -0.5; + IpV +} + +fn calculateJppW(b: &ConstrainedCurve) -> MatD { + let p_ = b.curve.degree(); + let o_ = b.curve.segments(); + let Wk0 = b.curve.knots.vector(); + + let oh_ = b.constraints.segments(); + + let mut JppW = MatD::zeros(p_, oh_ + 1); + + 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 *= -0.5; + JppW +} + +fn calculateKconst(a: &Curve, b: &Curve) -> MatD { + let p = a.degree(); + let dim = a.dimension(); + + let mut Kconst = MatD::zeros(dim, p); + let mut sum = VecD::zeros(dim); + + let m = a.segments(); + let o = b.segments(); + + let Vk = &a.knots.Uk; + let Wk = &b.knots.Uk; + + let S = a.points.matrix(); + let T = b.points.matrix(); + + for k in 0..=p - 1 { + sum.fill(0.0); + + 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); + } + } + + for j in 0..=p { + 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); + } + } + Kconst.column_mut(k).sub_assign(&sum); + } + Kconst +} + +fn constructConstMat(a: &Curve, b: &Curve, total_constraints: usize) -> MatD { + let p = a.degree(); + let dim = a.dimension(); + + let mut mat = MatD::zeros(dim, 3 * p + total_constraints); + + let Kconst = calculateKconst(a, b); + mat.view_mut((0, 2 * p), (dim, p)).copy_from(&Kconst); + + mat +} + +fn solve_linear_equation_system(a: &ConstrainedCurve, b: &ConstrainedCurve) -> MatD { + let total_constraints = a.constraints.count() + b.constraints.count(); + + let Nmat = construct_Nmat(a, b); + let const_mat = constructConstMat(a.curve, b.curve, total_constraints); + + SVD::new(Nmat, true, true).solve(&const_mat.transpose(), f64::EPSILON.sqrt()).unwrap().transpose() +} + +fn generate_shifted_control_point_vectors_of_spline1and2(a: &Curve, b: &Curve, shifts: &MatD) -> (MatD, MatD) { + let p = a.degree(); + let mut Sshifted = a.points.matrix().clone(); + let mut Tshifted = b.points.matrix().clone(); + + Sshifted.columns_mut(Sshifted.ncols() - p, p).add_assign(shifts.columns(0, p)); + + Tshifted.columns_mut(0, p).add_assign(shifts.columns(p, p)); + + (Sshifted, Tshifted) +} + +fn adjust_knots_of_both_splines(a: &Curve, b: &Curve) -> (VecD, VecD, VecD) { + let p = a.degree(); + + let m = a.segments(); + let o = b.segments(); + + let V0 = a.knots.vector(); + let W0 = b.knots.vector(); + + // Adjust left knots + let Vadj = adjust_knots(p, V0, m, W0); + + // Adjust right knots + let Vrev = reversed(V0); + let Wrev = reversed(W0); + let Wadjrev = adjust_knots(p, &Wrev, o, &Vrev); + let Wadj = reversed(&Wadjrev).add_scalar(1.); + + (Vadj, Wrev, Wadj) +} + +fn adjust_knots(p: usize, Uleft: &VecD, nLeft: usize, Uright: &VecD) -> VecD { + let mut UleftAdj = VecD::zeros(nLeft + p + 2); + + UleftAdj.head_mut(nLeft + 2).copy_from(&Uleft.head(nLeft + 2)); + + UleftAdj.tail_mut(p).copy_from(&Uright.segment(p + 1, p).add_scalar(1.)); + + UleftAdj +} + +fn create_merged_knot_vector(a: &Curve, b: &Curve, Vadj: &VecD, Wadj: &VecD) -> VecD { + let m = a.segments(); + let o = b.segments(); + + // construct the knot vector of the merged spline + let mut Umerged = VecD::zeros(m + 2 + o + 1); + + Umerged.head_mut(m + 2).copy_from(&Vadj.head(m + 2)); + Umerged.tail_mut(o + 1).copy_from(&Wadj.tail(o + 1)); + + // rescale the knot vector of the merged spline (u in [0,2]) to [0,1] by dividing through the last element + Umerged.div_assign(Umerged[m + o + 2]); + + Umerged +} + +fn generateDerivativeControlPoint( + idx: usize, + k: usize, + Pshifted: &MatD, + U0: &VecD, + p: usize, + n: usize, + dim: usize, +) -> VecD { + let mut controlPoint = VecD::zeros(dim); + + assert!(idx <= n - k, "Index out of bounds: only n-k control points exist"); + for zeroOrderIdx in 0..=n { + controlPoint += prefactor(p, idx, zeroOrderIdx, k, Pshifted, U0) * Pshifted.column(zeroOrderIdx); + } + controlPoint +} + +fn adjust_shifted_control_points( + Pshifted: &MatD, //shifted_control_points: &MatD, + U0: &VecD, //original_knot_vector: &VecD, + Uadj: &VecD, //adjusted_knot_vector: &VecD, + p: usize, + n: usize, + dim: usize, +) -> MatD { + let mut padj = MatD::zeros(dim, n + 1); + let mut qki: Vec> = vec![Vec::new(); n + 1]; + + for (i, elem) in qki.iter_mut().enumerate().take(n + 1) { + elem.push(Pshifted.column(i).into()); + } + + for k in 0..=p - 1 { + let derivative_point = generateDerivativeControlPoint(n - p + 1, k, Pshifted, U0, p, n, dim); //generateDerivativeControlPoint(n - p + 1, k, shifted_control_points, original_knot_vector); + qki[n - p + 1].push(derivative_point); + } + + for i in n - p + 2..=n { + qki[i].resize(n - i + 1, VecD::zeros(dim)); + for k in (0..=n - i).rev() { + qki[i][k] = ((Uadj[i + p] - Uadj[i + k]) / ((p - k) as f64)) * &qki[i - 1][k + 1] + &qki[i - 1][k]; + } + } + + for (i, elem) in qki.iter_mut().enumerate().take(n + 1) { + padj.set_column(i, &elem[0]); + } + + padj +} + +fn adjust_shifted_control_points_of_both_splines( + a: &Curve, + Sshifted: &MatD, + Tshifted: &MatD, + Vadj: &VecD, + Wrev: &VecD, + Wadjrev: &VecD, +) -> (MatD, MatD) { + let p = a.degree(); + let n = a.segments(); + let dim = a.dimension(); + + let V0 = a.knots.vector(); + + // obtain adjusted, shifted control points of the second spline + let Sadj = adjust_shifted_control_points(Sshifted, V0, Vadj, p, n, dim); + + // reverse control points of the second spline + let P2shiftedrev = curve::points::reversed(Tshifted); + + // obtain adjusted, shifted control points of the reversed second spline + let Tadjrev = adjust_shifted_control_points(&P2shiftedrev, Wrev, Wadjrev, p, n, dim); + + // obtain adjusted, shifted control points of the second spline + let Tadj = curve::points::reversed(&Tadjrev); + + (Sadj, Tadj) +} + +fn generate_control_point_vector_of_merged_spline(a: &Curve, b: &Curve, Sadj: &MatD, Tadj: &MatD) -> MatD { + let p = a.degree(); + let dim = a.dimension(); + let n_points1 = a.points.count(); + let n_points2 = b.points.count(); + + let mut P_merged = MatD::zeros(dim, n_points1 + n_points2 - p); + + P_merged.columns_mut(0, n_points1).copy_from(Sadj); + + let right_cols = n_points2 + 1 - p; + P_merged + .columns_mut(P_merged.ncols() - right_cols, right_cols) + .copy_from(&Tadj.columns(Tadj.ncols() - right_cols, right_cols)); + + P_merged +} + +fn kronecker_delta(i: usize, j: usize) -> bool { + i == j +} + +/// `p` degree +/// `i` index +/// `i0` 0th order index +/// `k` derivative order +/// `U0` 0th order knot vector +/// `P0` 0th order control point matrix +fn prefactor(p: usize, i: usize, i0: usize, k: usize, P0: &MatD, U0: &VecD) -> f64 { + let n = P0.ncols() - 1; // TODO use points + + if i <= n - k { + if k == 0 { + if kronecker_delta(i, i0) { + 1. + } else { + 0. + } + } 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)) + } + } else { + 0. + } +} + +#[cfg(test)] +mod tests { + use nalgebra::{dmatrix, dvector}; + + use crate::curve::{ + generation::{generate, Generation::Manual}, + knots::Generation::Uniform, + Curve, + }; + + use super::*; + + fn test_bspline(degree: usize, points: MatD) -> Curve { + generate(Manual { degree, points: ControlPoints::new(points), knots: Uniform }).unwrap() + } + + mod knots { + use super::*; + + #[test] + fn knots() { + let c = merge(&test_bspline(1, dmatrix![-2.,-1.,0.;]), &test_bspline(1, dmatrix![0.,1.,2.;])).unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0.25, 0.5, 0.75, 1., 1.]); + } + } + + mod control_points { + use std::ops::Mul; + + use approx::assert_relative_eq; + + use crate::curve::points; + + use super::*; + + #[test] + fn no_shift_degree_1() { + let p = 1; + let mat = dmatrix![ + 0.,1.,2.; + 0.,1.,2.; + ]; + let c = merge(&test_bspline(p, points::reversed(&mat).mul(-1.)), &test_bspline(p, mat)).unwrap(); + assert_relative_eq!( + c.points.matrix(), + &dmatrix![ + -2.,-1.,0.,1.,2.; + -2.,-1.,0.,1.,2.; + ], + epsilon = f64::EPSILON.sqrt() + ); + } + + #[test] + fn no_shift_degree_2() { + let p = 2; + let mat = dmatrix![0.,1.,2.;]; + let c = merge(&test_bspline(p, points::reversed(&mat).mul(-1.)), &test_bspline(p, mat)).unwrap(); + assert_relative_eq!(c.points.matrix(), &dmatrix![-2.,-1.,1.,2.;], epsilon = f64::EPSILON.sqrt()); + } + + #[test] + fn shift_degree_1() { + let p = 1; + let mat = dmatrix![0.5,1.,2.;]; + let c = merge(&test_bspline(p, points::reversed(&mat).mul(-1.)), &test_bspline(p, mat)).unwrap(); + assert_relative_eq!(c.points.matrix(), &dmatrix![-2.,-1.,0.,1.,2.;], epsilon = f64::EPSILON.sqrt()); + } + + #[test] + fn shift_degree_2() { + let p = 2; + let mat = dmatrix![0.5,1.,2.;]; + let c = merge(&test_bspline(p, points::reversed(&mat).mul(-1.)), &test_bspline(p, mat)).unwrap(); + assert_relative_eq!(c.points.matrix(), &dmatrix![-2.,-1.,1.,2.;], epsilon = f64::EPSILON.sqrt()); + } + + #[test] + fn shift_constrain_left_degree_2() { + let p = 2; + let mat = dmatrix![0.5,1.,2.;]; + let c = merge_from(&test_bspline(p, points::reversed(&mat).mul(-1.)), &test_bspline(p, mat)).unwrap(); + + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 0.5, 1., 1., 1.]); + + assert_relative_eq!(c.points.matrix(), &dmatrix![-2.,-1.5,0.5,2.;], epsilon = f64::EPSILON.sqrt()); + } + + #[test] + fn shift_constrain_right_degree_2() { + let p = 2; + let mat = dmatrix![0.5,1.,2.;]; + let c = merge_to(&test_bspline(p, points::reversed(&mat).mul(-1.)), &test_bspline(p, mat)).unwrap(); + + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 0.5, 1., 1., 1.]); + + assert_relative_eq!(c.points.matrix(), &dmatrix![-2.,-0.5,1.5,2.;], epsilon = f64::EPSILON.sqrt()); + } + } +} diff --git a/src/manipulation/mod.rs b/src/manipulation/mod.rs new file mode 100644 index 0000000..9bac688 --- /dev/null +++ b/src/manipulation/mod.rs @@ -0,0 +1,4 @@ +pub mod insert; +pub mod merge; +pub mod reverse; +pub mod split; diff --git a/src/manipulation/reverse.rs b/src/manipulation/reverse.rs new file mode 100644 index 0000000..951ea7d --- /dev/null +++ b/src/manipulation/reverse.rs @@ -0,0 +1,24 @@ +#![cfg_attr(feature = "doc-images", +cfg_attr(all(), +doc = ::embed_doc_image::embed_image!("reverse-before", "doc-images/plots/manipulation/reverse-before.svg"), +doc = ::embed_doc_image::embed_image!("reverse-after", "doc-images/plots/manipulation/reverse-after.svg")))] +//! Reverses the curve parametrization. +//! +//! | A normal curve. | The reversed curve. | +//! |:--------------------:|:-------------------:| +//! | ![][reverse-before] | ![][reverse-after] | + +use crate::curve::Curve; + +pub fn reverse(curve: &mut Curve) -> &mut Curve { + curve.knots.reverse(); + curve.points.reverse(); + curve +} + +// TODO +pub fn reversed(curve: &Curve) -> Curve { + let mut clone = curve.clone(); + clone.reverse(); + clone +} diff --git a/src/manipulation/split.rs b/src/manipulation/split.rs new file mode 100644 index 0000000..9c98400 --- /dev/null +++ b/src/manipulation/split.rs @@ -0,0 +1,218 @@ +#![cfg_attr(feature = "doc-images", +cfg_attr(all(), +doc = ::embed_doc_image::embed_image!("split-before", "doc-images/plots/manipulation/split-before.svg"), +doc = ::embed_doc_image::embed_image!("split-after", "doc-images/plots/manipulation/split-after.svg")))] +//! Splits a curve into two independent ones. +//! +//! | A curve before splitting. | The two independent curves after splitting at `$u=1/2$`. | +//! |:-------------------------:|:--------------------------------------------------------:| +//! | ![][split-before] | ![][split-after] | +//! +//! The splitting is conducted by adding the respective knot `$p+1$`-times, which allows for splitting the knot vector. +//! The knot vector can then be re-normalized on the interval `$[0,1]$`. + +use thiserror::Error; + +use crate::{ + curve::{ + knots::{normalize, DomainKnotComparatorType, Knots}, + points::{ControlPoints, Points}, + Curve, CurveError, + }, + manipulation::insert::insert, + types::{MatD, VecD, VecHelpers}, +}; + +#[derive(Error, Debug, PartialEq)] +pub enum SplitError { + #[error("The curve cannot be disconnected. The Multiplicity of `{multiplicity}` at u = {u}` exceeds `p = {p}.")] + MultiplicityTooHigh { u: f64, p: usize, multiplicity: usize }, + #[error("Parameter `u = {u}` lies outside the interval `({lower_bound}, {upper_bound})`.")] + OutOfBounds { u: f64, lower_bound: f64, upper_bound: f64 }, + + #[error("Curve generation failed with error.")] + CurveError(#[from] CurveError), +} + +pub fn split(c: &Curve, u: f64) -> Result<(Curve, Curve), SplitError> { + split_and_normalize(c, u, (true, true)) +} + +pub fn split_and_normalize( + c: &Curve, + u: f64, + normalizeKnotVectors: (bool, bool), +) -> Result<(Curve, Curve), SplitError> { + if u <= 0.0 || u >= 1.0 { + return Err(SplitError::OutOfBounds { u, lower_bound: 0.0, upper_bound: 1.0 }); + } + + let p = c.degree(); + + let mut bs_inserted = c.clone(); + + let l = c.knots.find_idx(u, 0, DomainKnotComparatorType::LeftOrEqual); + let multiplicity = c.knots.vector().iter().skip(l).take_while(|&&x| x == u).count(); + + if multiplicity > p { + return Err(SplitError::MultiplicityTooHigh { u, p, multiplicity }); + } + + for _ in 0..p - multiplicity { + insert(&mut bs_inserted, u).unwrap(); + } + + let U = bs_inserted.knots.vector(); + let P = bs_inserted.points.matrix(); + + if multiplicity > 0 { + let left = { + // TODO reduce index calcs + let mut U_left = VecD::zeros(l + p + 1); + U_left.head_mut(l + p).copy_from(&U.head(l + p)); + U_left[l + p] = u; + + if normalizeKnotVectors.0 { + normalize(&mut U_left); + } + let bot_cols = U_left.len() - (p + 2) + 1; + let P_left: MatD = P.columns(0, bot_cols).into(); + + Curve::new(Knots::new(p, U_left), ControlPoints::new(P_left))? + }; + + let right = { + let mut U_right = VecD::zeros(U.len() + 1 - l); + U_right[0] = u; + U_right.tail_mut(U.len() - l).copy_from(&U.tail(U.len() - l)); + + if normalizeKnotVectors.1 { + normalize(&mut U_right); + } + let bot_cols = U_right.len() - (p + 2) + 1; + let P_right: MatD = P.columns(P.ncols() - bot_cols, bot_cols).into(); + + Curve::new(Knots::new(p, U_right), ControlPoints::new(P_right))? + }; + Ok((left, right)) + } else { + let left = { + let mut U_left = VecD::zeros(l + p + 1 + 1); + U_left.head_mut(l + p + 1).copy_from(&U.head(l + p + 1)); + U_left[l + p + 1] = u; + + if normalizeKnotVectors.0 { + normalize(&mut U_left); + } + + let top_cols = U_left.len() + 1 - (p + 2); + let P_left: MatD = P.columns(0, top_cols).into(); + + Curve::new(Knots::new(p, U_left), ControlPoints::new(P_left))? + }; + + let right = { + let mut U_right = VecD::zeros(U.len() + 1 - (l + 1)); // length of U - the elements that occur before the split idx + U_right[0] = u; + U_right.tail_mut(U.len() - (l + 1)).copy_from(&U.tail(U.len() - (l + 1))); + + if normalizeKnotVectors.1 { + normalize(&mut U_right); + } + let bot_cols = U_right.len() + 1 - (p + 2); + let P_right: MatD = P.columns(P.ncols() - bot_cols, bot_cols).into(); + + Curve::new(Knots::new(p, U_right), ControlPoints::new(P_right))? + }; + Ok((left, right)) + } +} + +#[cfg(test)] +mod tests { + use approx::assert_relative_eq; + use nalgebra::{dmatrix, dvector}; + use rstest::{fixture, rstest}; + + use crate::curve::{ + generation::{generate, Generation::Manual}, + knots::Generation::Uniform, + }; + + use super::*; + + #[fixture] + /// A one-dimensional, linear test curve with default degree two. + fn c(#[default(2)] degree: usize) -> Curve { + let c = + generate(Manual { degree, points: ControlPoints::new(dmatrix![1., 2., 3., 4., 5., 6.;]), knots: Uniform }) + .unwrap(); + assert_eq!(c.knots.vector(), &dvector![0., 0., 0., 0.25, 0.5, 0.75, 1., 1., 1.]); + c + } + + #[rstest] + fn cannot_split_start(c: Curve) { + let u = 0.0; + let res = split_and_normalize(&c, u, (true, true)); + assert!(res.is_err()); + assert_eq!(res.unwrap_err(), SplitError::OutOfBounds { u, lower_bound: 0.0, upper_bound: 1.0 }); + } + + #[rstest] + fn cannot_split_end(c: Curve) { + let u = 1.0; + let res = split_and_normalize(&c, u, (true, true)); + assert!(res.is_err()); + assert_eq!(res.unwrap_err(), SplitError::OutOfBounds { u, lower_bound: 0.0, upper_bound: 1.0 }); + } + + /*#[rstest] + fn cannot_split_wrong_mult() { + todo!("add multiplicity test"); + }*/ + + #[rstest] + fn split_close_to_start(c: Curve) { + let eps = f64::EPSILON; + let (left, right) = split_and_normalize(&c, 0.0 + eps, (true, true)).unwrap(); + + assert_relative_eq!(left.knots.vector(), &dvector![0., 0., 0., 1., 1., 1.], epsilon = eps.sqrt()); + assert_relative_eq!(left.points.matrix(), &dmatrix![1., 1., 1.;], epsilon = eps.sqrt()); + + assert_relative_eq!(right.knots.vector(), c.knots.vector(), epsilon = eps.sqrt()); + assert_relative_eq!(right.points.matrix(), c.points.matrix(), epsilon = eps.sqrt()); + } + + #[rstest] + fn split_close_to_end(c: Curve) { + let eps = f64::EPSILON; + let (left, right) = split_and_normalize(&c, 1.0 - eps, (true, true)).unwrap(); + + assert_relative_eq!(left.knots.vector(), c.knots.vector(), epsilon = eps.sqrt()); + assert_relative_eq!(left.points.matrix(), c.points.matrix(), epsilon = eps.sqrt()); + + assert_relative_eq!(right.knots.vector(), &dvector![0., 0., 0., 1., 1., 1.], epsilon = eps.sqrt()); + assert_relative_eq!(right.points.matrix(), &dmatrix![6., 6., 6.;], epsilon = eps.sqrt()); + } + + #[rstest] + fn normalized(c: Curve) { + let (left, right) = split_and_normalize(&c, 0.5, (true, true)).unwrap(); + assert_eq!(left.knots.vector(), &dvector![0., 0., 0., 0.5, 1., 1., 1.]); + assert_eq!(left.points.matrix(), &dmatrix![1., 2., 3., 3.5;]); + + assert_eq!(right.knots.vector(), &dvector![0., 0., 0., 0.5, 1., 1., 1.]); + assert_eq!(right.points.matrix(), &dmatrix![3.5, 4., 5., 6.;]); + } + + #[rstest] + fn unnormalized(c: Curve) { + let (left, right) = split_and_normalize(&c, 0.5, (false, false)).unwrap(); + assert_eq!(left.knots.vector(), &dvector![0., 0., 0., 0.25, 0.5, 0.5, 0.5]); + assert_eq!(left.points.matrix(), &dmatrix![1., 2., 3., 3.5;]); + + assert_eq!(right.knots.vector(), &dvector![0.5, 0.5, 0.5, 0.75, 1., 1., 1.]); + assert_eq!(right.points.matrix(), &dmatrix![3.5, 4., 5., 6.;]); + } +} diff --git a/src/types.rs b/src/types.rs new file mode 100644 index 0000000..bec4c10 --- /dev/null +++ b/src/types.rs @@ -0,0 +1,77 @@ +use nalgebra::{Dyn, Matrix, MatrixView, MatrixViewMut, OMatrix, OVector, Owned, U1}; + +pub type VecD = OVector; +pub type RowVecD = Matrix>; + +pub type VecDView<'a> = MatrixView<'a, f64, Dyn, U1, U1, Dyn>; +pub type VecDViewMut<'a> = MatrixViewMut<'a, f64, Dyn, U1, U1, Dyn>; + +pub type MatD = OMatrix; +pub type MatDView<'a> = MatrixView<'a, f64, Dyn, Dyn>; +pub type MatRowD = OMatrix; + +pub type KnotVectorDerivatives = Vec; +pub type ControlPointDerivatves = Vec; + +pub trait VecHelpers { + fn head(&self, n: usize) -> MatrixView; + fn head_mut(&mut self, n: usize) -> MatrixViewMut; + + fn segment(&self, i: usize, n: usize) -> MatrixView; + fn segment_mut(&mut self, i: usize, n: usize) -> MatrixViewMut; + + fn tail(&self, n: usize) -> MatrixView; + fn tail_mut(&mut self, n: usize) -> MatrixViewMut; +} + +impl VecHelpers for VecD { + fn head(&self, n: usize) -> MatrixView { + self.segment(0, n) + } + + fn head_mut(&mut self, n: usize) -> MatrixViewMut { + self.segment_mut(0, n) + } + + fn segment(&self, start: usize, n: usize) -> MatrixView { + self.generic_view((start, 0), (Dyn(n), U1)) + } + + fn segment_mut(&mut self, start: usize, n: usize) -> MatrixViewMut { + self.generic_view_mut((start, 0), (Dyn(n), U1)) + } + + fn tail(&self, n: usize) -> MatrixView { + self.segment(self.len() - n, n) + } + + fn tail_mut(&mut self, n: usize) -> MatrixViewMut { + self.segment_mut(self.len() - n, n) + } +} + +#[cfg(test)] +mod vec_helpers { + use nalgebra::dvector; + + use super::*; + + fn example() -> VecD { + dvector![0.0, 1.0, 2.0, 3.0] + } + + #[test] + fn head() { + assert_eq!(example().head(2).as_slice(), [0.0, 1.0]); + } + + #[test] + fn segment() { + assert_eq!(example().segment(1, 2).as_slice(), [1.0, 2.0]); + } + + #[test] + fn tail() { + assert_eq!(example().tail(2).as_slice(), [2.0, 3.0]); + } +}