Commit 25dcaea8 authored by Egor Larionov's avatar Egor Larionov

Upgraded HRBF interface to be more Rustic + docs

Renamed HRBF to Hrbf in all type occurances to comply with Rust API
guidelines.

Built a builder interface for building HRBFs. This way it is impossible
to create an invalid/empty HRBF field. The intent is now better
represented in the interface.

The above changes also force the bundled C bindings to change slightly:
building an HRBF potential or the corresponding system is now done with
on C API call instead of two.

The C API now uses enums instead of ints to select kernel types which is
more idomatic in both C and Rust.

All dependencies were upgraded to latest versions.

Crate level documentation was added with a basic representative example.
Docs were also formatted better.

A warn(missing_docs) was added to lib.rs to keep track that
documentation is complete.

Bumped crate to v0.3.0.
parent e21fae99
Pipeline #126890125 passed with stage
in 2 minutes and 42 seconds
[package]
name = "hrbf"
version = "0.2.1"
version = "0.3.0"
authors = ["Egor Larionov <[email protected]>"]
license = "MPL-2.0"
description = "An implementation of Hermite Radial Basis Functions with higher order derivatives"
......@@ -19,9 +19,9 @@ maintenance = { status = "passively-maintained" }
[dependencies]
num-traits = "0.2"
na = { package = "nalgebra", version = "0.18" }
na = { package = "nalgebra", version = "0.20" }
[dev-dependencies]
approx = "0.3"
autodiff = "0.1"
rand = "0.6"
rand = "0.7"
[package]
name = "chrbf"
version = "0.1.0"
version = "0.2.0"
authors = ["Egor Larionov <[email protected]>"]
edition = "2018"
[lib]
name = "chrbf"
......@@ -9,7 +10,7 @@ crate-type = ["cdylib", "rlib"] # rlib is needed for testing
[dependencies]
hrbf = { path = "../." }
nalgebra = "*"
na = { package = "nalgebra", version = "*" }
libc = "*"
[build-dependencies]
......
This diff is collapsed.
extern crate itertools;
extern crate libc;
extern crate nalgebra as na;
#[macro_use]
extern crate approx;
extern crate chrbf;
use chrbf::*;
use itertools::Itertools;
use libc::{c_double, size_t};
......@@ -54,8 +50,8 @@ fn test_hrbf_fit() {
unsafe {
let n: size_t = pts.len() / 3;
let hrbf = HRBF_create(n, pts.as_ptr(), 0);
assert!(HRBF_fit(hrbf, n, pts.as_ptr(), nmls.as_ptr()));
let hrbf = HRBF_create_fit(n, pts.as_ptr(), nmls.as_ptr(), KernelType::HRBF_POW3);
assert_ne!(hrbf, std::ptr::null_mut() as *mut hrbf::Pow3Hrbf<f64>);
for ((&p0, &p1, &p2), (&n0, &n1, &n2)) in pts.iter().tuples().zip(nmls.iter().tuples()) {
let p = Point3::new(p0, p1, p2);
......
//! Hermite radial basis function kernels (funcions φ)
//! We require that the first derivative of the kernel f go to zero as x -> 0.
//! This ensures that the hrbf fitting matrix is well defined.
//! I.e. in its taylor series representation, f(x) = ∑ᵢ aᵢxⁱ, we require that a₁ = 0.
//! Hermite radial basis function kernels (funcions `φ`)
//! We require that the first derivative of the kernel `f` go to zero as `x -> 0`.
//! This ensures that the HRBF fitting matrix is well defined.
//! I.e. in its taylor series representation, `f(x) = ∑ᵢ aᵢxⁱ`, we require that `a₁ = 0`.
//!
//! NOTE: Kernels like x^2, x^3, gauss, and (1-x)^4 (4x+1) all satisfy this criterion.
//! NOTE: Kernels like `x^2`, `x^3`, `exp(-x*x)`, and `(1-x)^4 (4x+1)` all satisfy this criterion.
//!
//! To ensure that derivativevs at zero are computed accurately, each kernel must provide a formula
//! for the following additional function:
//!
//! ```verbatim
//! df_l(x) = φʹ(x)/x
//! where φ is the kernel function. It must be that df(x), ddf(x) - df_l(x) -> 0 as x -> 0 for the HRBF
//! ```
//! where `φ` is the kernel function. It must be that `df(x)`, `ddf(x) - df_l(x) -> 0` as `x -> 0`
//! for the HRBF
//! derivaitves to exist.
//!
//!
//! Furthermore, if the HRBF is to be used in the optimization framework,
//! where a 3rd or even 4th derivatives are required, we also need the 3rd derivative to vanish.
//! I.e. in its taylor series representation, f(x) = ∑ᵢ aᵢxⁱ, we require that a₁ = a₃ = 0.
//! I.e. in its taylor series representation, `f(x) = ∑ᵢ aᵢxⁱ`, we require that `a₁ = a₃ = 0`.
//!
//! NOTE: The gauss and x^2 kernels satisfy this criteria, but x^3, x^2*log(x) and (1-x)^4 (4x+1) do not!
//! NOTE: The `exp(-x*x)` and `x^2` kernels satisfy this criteria, but `x^3, x^2*log(x)` and
//! `(1-x)^4 (4x+1)` do not!
//!
//! To ensure that derivatives at zero are computed accurately, each kernel must provide a formula for
//! the function df_l (as described above) and the following additional functions:
//! the function `df_l` (as described above) and the following additional functions:
//!
//! ```verbatim
//! g(x) = φʺ(x)/x - φʹ(x)/x^2
//! g_l(x) = φʺ(x)/x^2 - φʹ(x)/x^3
//! h(x,a) = φ‴(x)/x - a(φʺ(x)/x^2 - φʹ(x)/x^3)
//! ```
//!
//! to see how these are used, see mesh_implicit_surface.cpp
//! It must be that g(x), h(x,3) -> 0 as x -> 0 for the HRBF derivatives to exist.
//! It must be that `g(x), h(x,3) -> 0` as `x -> 0` for the HRBF derivatives to exist.
use num_traits::Float;
use std::marker::PhantomData;
......@@ -34,36 +43,36 @@ pub trait Kernel<T>
where
T: Float,
{
/// Main kernel function φ(x)
/// Main kernel function `φ(x)`
fn f(&self, x: T) -> T;
/// The first derivative φʹ(x)
/// The first derivative `φʹ(x)`
fn df(&self, x: T) -> T;
/// The second derivative φʺ(x)
/// The second derivative `φʺ(x)`
fn ddf(&self, x: T) -> T;
/// The third derivative of φ(x)
/// The third derivative of `φ(x)`
fn dddf(&self, x: T) -> T;
/// The fourth derivative of φ(x)
/// The fourth derivative of `φ(x)`
fn ddddf(&self, x: T) -> T;
/// Additional function to ensure proper derivatives at x = 0.
/// equivalent to df(x)/x for x != 0.
/// This function should be well defined at all values of x.
/// Additional function to ensure proper derivatives at `x = 0`.
/// equivalent to `df(x)/x` for `x != 0`.
/// This function should be well defined at all values of `x`.
fn df_l(&self, x: T) -> T;
/// Need the following functions for third and fourth hrbf derivaitves
/// Need the following functions for third and fourth HRBF derivaitves
///
/// Additional function to ensure proper derivatives at x = 0.
/// equivalent to ddf(x)/x - df(x)/(x*x) for x != 0.
/// This function should go to zero as x goes to zero.
/// Additional function to ensure proper derivatives at `x = 0`.
/// equivalent to `ddf(x)/x - df(x)/(x*x)` for `x != 0`.
/// This function should go to zero as `x` goes to zero.
fn g(&self, x: T) -> T;
/// Additional function to ensure proper derivatives at x = 0.
/// equivalent to ddf(x)/(x*x) - df(x)/(x*x*x) for x != 0.
/// This function should be well defined at all values of x.
/// Additional function to ensure proper derivatives at `x = 0`.
/// equivalent to `ddf(x)/(x*x) - df(x)/(x*x*x)` for `x != 0`.
/// This function should be well defined at all values of `x`.
fn g_l(&self, x: T) -> T;
/// Additional function to ensure proper derivatives at x = 0.
/// equivalent to dddf(x)/x - a*(ddf(x)/(x*x) - df(x)/(x*x*x)) for x != 0.
/// Additional function to ensure proper derivatives at `x = 0`.
/// equivalent to `dddf(x)/x - a*(ddf(x)/(x*x) - df(x)/(x*x*x))` for `x != 0`.
/// This function should go to zero as x goes to zero.
fn h(&self, x: T, a: T) -> T;
}
......@@ -73,14 +82,17 @@ pub trait LocalKernel<T>
where
T: Float,
{
/// Construct a new local kernel with a radius `r`.
fn new(r: T) -> Self;
}
/// Global kernel trait defines a constructor for kernels without a radial fallofff.
pub trait GlobalKernel {
/// Construct a new global kernel.
fn new() -> Self;
}
/// The `x^2` kernel.
#[derive(Copy, Clone)]
pub struct Pow2<T>(PhantomData<T>);
......@@ -127,6 +139,7 @@ impl<T: Float> Kernel<T> for Pow2<T> {
}
}
/// The `x^3` kernel.
#[derive(Copy, Clone)]
pub struct Pow3<T>(::std::marker::PhantomData<T>);
......@@ -179,6 +192,7 @@ impl<T: Float> Kernel<T> for Pow3<T> {
}
}
/// The `x^4` kernel.
#[derive(Copy, Clone)]
pub struct Pow4<T>(::std::marker::PhantomData<T>);
......@@ -225,7 +239,7 @@ impl<T: Float> Kernel<T> for Pow4<T> {
}
}
/// x^5 kernel.
/// The `x^5` kernel.
#[derive(Copy, Clone)]
pub struct Pow5<T>(::std::marker::PhantomData<T>);
......@@ -272,7 +286,7 @@ impl<T: Float> Kernel<T> for Pow5<T> {
}
}
/// Gaussian kernel.
/// Gaussian `exp(-x*x)` kernel.
#[derive(Copy, Clone)]
pub struct Gauss<T> {
r: T,
......@@ -337,8 +351,10 @@ impl<T: Float> Kernel<T> for Gauss<T> {
}
}
/// Quintic kernel. Generates a positive definite hrbf fitting matrix.
/// Third and fourth order hrbf derivatives don't exist at x = 0.
/// Quintic kernel `(1-x)^4 (4x+1)`.
///
/// Generates a positive definite HRBF fitting matrix.
/// Third and fourth order HRBF derivatives don't exist at `x = 0`.
#[derive(Copy, Clone)]
pub struct Csrbf31<T> {
r: T,
......@@ -460,6 +476,9 @@ impl<T: Float> Kernel<T> for Csrbf31<T> {
}
}
/// The `(1-x)^6 (35x^2 + 18x + 3)` kernel.
///
/// This is a higher order version of the `Csrbf31` kernel with similar properties.
#[derive(Copy, Clone)]
pub struct Csrbf42<T> {
r: T,
......
This diff is collapsed.
......@@ -14,22 +14,22 @@ fn cube() -> (Vec<Point3<f64>>, Vec<Vector3<f64>>) {
// Fit an hrbf surface to a unit box
let pts = vec![
// Corners of the box
Point3::new(0.0, 0.0, 0.0),
Point3::new(0.0, 0.0, 1.0),
Point3::new(0.0, 1.0, 0.0),
Point3::new(0.0, 1.0, 1.0),
Point3::new(1.0, 0.0, 0.0),
Point3::new(1.0, 0.0, 1.0),
Point3::new(1.0, 1.0, 0.0),
Point3::new(1.0, 1.0, 1.0),
[0.0, 0.0, 0.0],
[0.0, 0.0, 1.0],
[0.0, 1.0, 0.0],
[0.0, 1.0, 1.0],
[1.0, 0.0, 0.0],
[1.0, 0.0, 1.0],
[1.0, 1.0, 0.0],
[1.0, 1.0, 1.0],
// Extra vertices on box faces
Point3::new(0.5, 0.5, 0.0),
Point3::new(0.5, 0.5, 1.0),
Point3::new(0.5, 0.0, 0.5),
Point3::new(0.5, 1.0, 0.5),
Point3::new(0.0, 0.5, 0.5),
Point3::new(1.0, 0.5, 0.5),
];
[0.5, 0.5, 0.0],
[0.5, 0.5, 1.0],
[0.5, 0.0, 0.5],
[0.5, 1.0, 0.5],
[0.0, 0.5, 0.5],
[1.0, 0.5, 0.5],
].into_iter().map(Point3::from).collect();
let a = 1.0f64 / 3.0.sqrt();
let nmls = vec![
......@@ -55,14 +55,13 @@ fn cube() -> (Vec<Point3<f64>>, Vec<Vector3<f64>>) {
}
/// Finite difference test
fn test_derivative_fd<F, K: Kernel<f64> + Default>(x: Point3<f64>, compare: F, order: usize)
fn test_derivative_fd<F, K: Kernel<f64> + Clone + Default>(x: Point3<f64>, compare: F, order: usize)
where
F: Fn(f64, f64),
{
let (pts, nmls) = cube();
let mut hrbf = HRBF::<f64, K>::new(pts.clone());
assert!(hrbf.fit(&pts, &nmls));
let hrbf = HrbfBuilder::<_, K>::new(pts.clone()).normals(nmls).build().unwrap();
// We will test hrbf derivatives using central differencing away from zeros since autodiff
// fails in these scenarios because it can't simplify expressions with division by zero.
......@@ -100,7 +99,7 @@ where
}
}
fn test_hrbf_derivative_simple<K: Kernel<f64> + Default>(order: usize) {
fn test_hrbf_derivative_simple<K: Kernel<f64> + Clone + Default>(order: usize) {
let test_pts = [
Point3::new(0.5, 1.0, 0.5),
Point3::new(0.1, 0.0, 0.0),
......@@ -113,7 +112,7 @@ fn test_hrbf_derivative_simple<K: Kernel<f64> + Default>(order: usize) {
}
}
fn test_hrbf_derivative_random<K: Kernel<f64> + Default>(order: usize) {
fn test_hrbf_derivative_random<K: Kernel<f64> + Clone + Default>(order: usize) {
use rand::distributions::Uniform;
let mut rng: StdRng = SeedableRng::from_seed([3u8; 32]);
......@@ -124,11 +123,10 @@ fn test_hrbf_derivative_random<K: Kernel<f64> + Default>(order: usize) {
}
}
fn test_hrbf_fit<K: Kernel<f64> + Default>() {
fn test_hrbf_fit<K: Kernel<f64> + Clone + Default>() {
let (pts, nmls) = cube();
let mut hrbf = HRBF::<f64, K>::new(pts.clone());
assert!(hrbf.fit(&pts, &nmls));
let hrbf = HrbfBuilder::<_, K>::new(pts.clone()).normals(nmls.clone()).build().unwrap();
for (p, n) in pts.into_iter().zip(nmls) {
let p_u = p + 0.001 * n;
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment