Commit db0fa714 authored by Egor Larionov's avatar Egor Larionov

Add C API. Improve Docs. Smarter Traits.

parent ae660bc4
[package]
name = "chrbf"
version = "0.1.0"
authors = ["Egor Larionov <[email protected]>"]
[lib]
name = "chrbf"
crate-type = ["cdylib", "rlib"] # rlib is needed for testing
[dependencies]
hrbf = { path = "../." }
nalgebra = "*"
libc = "*"
[build-dependencies]
cbindgen = { git = "https://github.com/eqrion/cbindgen.git" }
[dev-dependencies]
approx = "0.1"
itertools = "0.7"
extern crate cbindgen;
use std::env;
fn main() {
let crate_dir = env::var("CARGO_MANIFEST_DIR").unwrap();
let mut config: cbindgen::Config = Default::default();
config.include_guard = Some(String::from("CHRBF_BINDINGS_H"));
config.namespace = Some(String::from("hrbf"));
config.line_length = 80;
config.tab_width = 2;
config.language = cbindgen::Language::Cxx;
cbindgen::generate_with_config(&crate_dir, config)
.expect("Unable to generate bindings")
.write_to_file("hrbf.h");
}
extern crate libc;
extern crate hrbf;
extern crate nalgebra as na;
use libc::{size_t, c_int, c_double};
use na::{
Point3,
DimName,
VectorN,
Point,
U3,
DefaultAllocator
};
use na::allocator::Allocator;
use hrbf::{
Pow3HRBF,
Pow5HRBF,
GaussHRBF,
Csrbf31HRBF,
Csrbf42HRBF,
HRBFTrait,
};
/// Create an HRBF field object. At this stage, the field is zero. You need to fit it to some data
/// using `HRBF_fit` to produce a useful field.
///
/// * `num_sites` is the number of sites needed to represent this HRBF.
/// * `sites_ptr` is the pointer to a contiguous array of 3 dimensional vectors of `double`s. the
/// size of this array must be `3*num_sites`.
/// * `kern_type` is an integer from 0 to 4 representing a particular kernel to use for fitting.
/// Use:
/// `0` for the cubic `x^3` kernel, (Default)
/// `1` for the quintic `x^5` kernel,
/// `2` for Gaussian `exp(-x^2)` kernel,
/// `3` for CSRBF31 kernel: `(1-x)^4 (4x + 1)` and
/// `4` for CSRBF42 kernel: `(1-x)^6 (35x^2 + 18x + 3)`.
#[no_mangle]
pub unsafe extern "C" fn HRBF_create(num_sites: size_t,
sites_ptr: *const c_double,
kern_type: c_int) -> *mut HRBFTrait<f64>
{
let sites = ptr_to_vec_of_points::<f64, U3>(num_sites, sites_ptr);
let hrbf: Box<HRBFTrait<f64>> =
match kern_type {
1 => Box::new(Pow5HRBF::new(sites)),
2 => Box::new(GaussHRBF::new(sites)),
3 => Box::new(Csrbf31HRBF::new(sites)),
4 => Box::new(Csrbf42HRBF::new(sites)),
_ => Box::new(Pow3HRBF::new(sites)),
};
Box::into_raw(hrbf)
}
/// Free memory allocated by HRBF_create.
#[no_mangle]
pub unsafe extern "C" fn HRBF_free(hrbf: *mut HRBFTrait<f64>) {
let _ = Box::from_raw(hrbf);
}
/// Given an HRBF object, fit it to the given set of points and normals.
/// NOTE: Currently this function only works for the same number of points as there are sites.
#[no_mangle]
pub unsafe extern "C" fn HRBF_fit(hrbf: *mut HRBFTrait<f64>,
num_pts: size_t,
pts_ptr: *const c_double,
nmls_ptr: *const c_double) -> bool {
let pts = ptr_to_vec_of_points::<f64, U3>(num_pts, pts_ptr);
let nmls = ptr_to_vec_of_vectors::<f64, U3>(num_pts, nmls_ptr);
(*hrbf).fit(&pts, &nmls)
}
/// Evaluate the given HRBF field at a particular point.
#[no_mangle]
pub unsafe extern "C" fn HRBF_eval(hrbf: *const HRBFTrait<f64>,
x: c_double,
y: c_double,
z: c_double) -> c_double
{
(*hrbf).eval(Point3::new(x,y,z))
}
/// Evaluate the gradient of the given HRBF field at a particular point.
/// `out` is expected to be an array of 3 doubles representing a 3D vector.
#[no_mangle]
pub unsafe extern "C" fn HRBF_grad(hrbf: *const HRBFTrait<f64>,
x: c_double,
y: c_double,
z: c_double,
out: *mut c_double)
{
let g = (*hrbf).grad(Point3::new(x,y,z));
*out.offset(0) = g[0];
*out.offset(1) = g[1];
*out.offset(2) = g[2];
}
/// Evaluate the hessian of the given HRBF field at a particular point.
/// `out` is expected to be an array of 9 doubles representing a 3x3 matrix.
#[no_mangle]
pub unsafe extern "C" fn HRBF_hess(hrbf: *const HRBFTrait<f64>,
x: c_double,
y: c_double,
z: c_double,
out: *mut c_double)
{
let h = (*hrbf).hess(Point3::new(x,y,z));
for i in 0..3 {
for j in 0..3 {
*out.offset(3*i + j) = h[(i as usize,j as usize)];
}
}
}
/// Advanced. Build the system for fitting an HRBF function. Don't do the actual fitting,
/// simply return the necessary fit matrix A and right-hand-side b, such that $A^{-1}*b$
/// gives the hrbf coefficients.
///
/// * `num_sites` is the number of sites needed to represent this HRBF.
/// * `sites_ptr`` is the pointer to a contiguous array of 3 dimensional vectors of `double`s. The
/// size of this array must be `3*num_sites`.
/// * `num_pts` is the number of data points to fit to.
/// * `pts_ptr`` is the pointer to a contiguous array of 3 dimensional vectors of `double`s. The
/// size of this array must be `3*num_pts`. This array represents the point positions to fit to.
/// * `nmls_ptr`` is the pointer to a contiguous array of 3 dimensional vectors of `double`s. The
/// size of this array must be `3*num_pts`. This array represents the normals to fit to.
/// * `kern_type` is an integer from 0 to 4 representing a particular kernel to use for fitting.
/// Use:
/// `0` for the cubic `x^3` kernel, (Default)
/// `1` for the quintic `x^5` kernel,
/// `2` for Gaussian `exp(-x^2)` kernel,
/// `3` for CSRBF31 kernel: `(1-x)^4 (4x + 1)` and
/// `4` for CSRBF42 kernel: `(1-x)^6 (35x^2 + 18x + 3)`.
/// * `A` output matrix expected to be a contiguous array of `4*num_sites*4*num_pts` doubles. This
/// array represents a `num_sites x num_pts` sized matrix and is stored in column major order.
/// * `b` output vector expected to be a contiguous array of `4*num_pts` doubles.
#[no_mangle]
#[allow(non_snake_case)]
pub unsafe extern "C" fn HRBF_fit_system(num_sites: size_t,
sites_ptr: *const c_double,
num_pts: size_t,
pts_ptr: *const c_double,
nmls_ptr: *const c_double,
kern_type: c_int,
A: *mut c_double,
b: *mut c_double)
{
let sites = ptr_to_vec_of_points::<f64, U3>(num_sites, sites_ptr);
let pts = ptr_to_vec_of_points::<f64, U3>(num_pts, pts_ptr);
let nmls = ptr_to_vec_of_vectors::<f64, U3>(num_pts, nmls_ptr);
let hrbf: Box<HRBFTrait<f64>> =
match kern_type {
1 => Box::new(Pow5HRBF::new(sites)),
2 => Box::new(GaussHRBF::new(sites)),
3 => Box::new(Csrbf31HRBF::new(sites)),
4 => Box::new(Csrbf42HRBF::new(sites)),
_ => Box::new(Pow3HRBF::new(sites)),
};
let mut offsets = Vec::new();
offsets.resize(num_pts, 0.0);
let (A_na, b_na) = (*hrbf).fit_system(&pts, &offsets, &nmls);
for col in 0..4*num_sites {
for row in 0..4*num_pts {
*A.offset((num_sites*col + row) as isize) = A_na[(row,col)];
}
}
for i in 0..4*num_pts {
*b.offset(i as isize) = b_na[i];
}
}
/// Helper routines for converting C-style data to Rust-style data.
/// `n` is the number of vectors to output, which means that `data_ptr` must point to an array of
/// `n*D::dim()` doubles.
unsafe fn ptr_to_vec_of_vectors<T: na::Real, D: DimName>(n: size_t, data_ptr: *const T) -> Vec<VectorN<T, D>>
where DefaultAllocator: Allocator<T, D>
{
let mut data = Vec::with_capacity(n);
for i in 0..n as isize {
data.push(VectorN::from_fn(|r,_| *data_ptr.offset(3*i + r as isize)));
}
data
}
unsafe fn ptr_to_vec_of_points<T: na::Real, D: DimName>(n: size_t, data_ptr: *const T) -> Vec<Point<T, D>>
where DefaultAllocator: Allocator<T, D>
{
ptr_to_vec_of_vectors(n, data_ptr).into_iter().map(|v| Point::from_coordinates(v)).collect()
}
extern crate nalgebra as na;
extern crate libc;
extern crate itertools;
#[macro_use] extern crate approx;
extern crate chrbf;
use na::{Point3, Vector3};
use itertools::Itertools;
use libc::{size_t, c_double};
use chrbf::*;
#[test]
fn test_hrbf_fit() {
// Fit an hrbf surface to a unit box
let pts: &[c_double] = &[
// Corners of the box
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
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,
];
let a: c_double = 1.0/c_double::sqrt(3.0);
let nmls: &[c_double] = &[
// Corner normals
-a, -a, -a,
-a, -a, a,
-a, a, -a,
-a, a, a,
a, -a, -a,
a, -a, a,
a, a, -a,
a, a, a,
// Side normals
0.0, 0.0, -1.0,
0.0, 0.0, 1.0,
0.0, -1.0, 0.0,
0.0, 1.0, 0.0,
-1.0, 0.0, 0.0,
1.0, 0.0, 0.0,
];
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()));
for ((&p0, &p1, &p2), (&n0, &n1, &n2)) in pts.iter().tuples().zip(nmls.iter().tuples()) {
let p = Point3::new(p0, p1, p2);
let n = Vector3::new(n0, n1, n2);
let p_u = p + n*0.001;
let p_l = p - n*0.001;
assert_relative_eq!(HRBF_eval(hrbf, p[0], p[1], p[2]), 0.0, max_relative=1e-6, epsilon=1e-11);
assert!(HRBF_eval(hrbf, p_l[0], p_l[1], p_l[2]) < 0.0);
assert!(HRBF_eval(hrbf, p_u[0], p_u[1], p_u[2]) > 0.0);
let mut g = [0.0f64; 3];
HRBF_grad(hrbf, p[0], p[1], p[2], g.as_mut_ptr());
assert_relative_eq!(g[0], n[0], max_relative=1e-6, epsilon=1e-11);
assert_relative_eq!(g[1], n[1], max_relative=1e-6, epsilon=1e-11);
assert_relative_eq!(g[2], n[2], max_relative=1e-6, epsilon=1e-11);
}
}
}
......@@ -30,7 +30,7 @@ use num_traits::{Float};
use std::marker::PhantomData;
/// Kernel trait declaring all of the necessary derivatives.
pub trait Kernel<T> : Copy + Clone + Default
pub trait Kernel<T>
where T: Float
{
/// Main kernel function φ(x)
......@@ -72,10 +72,16 @@ pub trait LocalKernel<T> where T: Float {
fn new(r: T) -> Self;
}
/// Global kernel trait defines a constructor for kernels without a radial fallofff.
pub trait GlobalKernel {
fn new() -> Self;
}
#[derive(Copy, Clone)]
pub struct Pow2<T>(PhantomData<T>);
impl<T: Float> GlobalKernel for Pow2<T> { fn new() -> Self { Pow2(PhantomData) } }
/// Default constructor for `Pow2` kernel
impl<T: Float> Default for Pow2<T> { fn default() -> Self { Pow2(PhantomData) } }
......@@ -95,6 +101,8 @@ impl<T: Float> Kernel<T> for Pow2<T> {
#[derive(Copy, Clone)]
pub struct Pow3<T>(::std::marker::PhantomData<T>);
impl<T: Float> GlobalKernel for Pow3<T> { fn new() -> Self { Pow3(PhantomData) } }
/// Default constructor for `Pow3` kernel
impl<T: Float> Default for Pow3<T> { fn default() -> Self { Pow3(PhantomData) } }
......@@ -120,6 +128,8 @@ impl<T: Float> Kernel<T> for Pow3<T> {
#[derive(Copy, Clone)]
pub struct Pow4<T>(::std::marker::PhantomData<T>);
impl<T: Float> GlobalKernel for Pow4<T> { fn new() -> Self { Pow4(PhantomData) } }
/// Default constructor for `Pow4` kernel
impl<T: Float> Default for Pow4<T> { fn default() -> Self { Pow4(PhantomData) } }
......@@ -140,6 +150,8 @@ impl<T: Float> Kernel<T> for Pow4<T> {
#[derive(Copy, Clone)]
pub struct Pow5<T>(::std::marker::PhantomData<T>);
impl<T: Float> GlobalKernel for Pow5<T> { fn new() -> Self { Pow5(PhantomData) } }
/// Default constructor for `Pow5` kernel
impl<T: Float> Default for Pow5<T> { fn default() -> Self { Pow5(PhantomData) } }
......
This diff is collapsed.
......@@ -7,7 +7,6 @@ extern crate approx;
use num_traits::Float;
use hrbf::*;
use hrbf::kernel::*;
use nalgebra::{Matrix3, Vector3, Point3};
fn rel_compare(a: f64, b: f64) {
......@@ -61,7 +60,7 @@ fn cube() -> (Vec<Point3<f64>>, Vec<Vector3<f64>>) {
}
/// Finite difference test
fn test_derivative_fd<F,K: Kernel<f64>>(x: Point3<f64>, compare: F, order: usize)
fn test_derivative_fd<F,K: Kernel<f64> + Default>(x: Point3<f64>, compare: F, order: usize)
where F: Fn(f64, f64)
{
let (pts, nmls) = cube();
......@@ -102,7 +101,7 @@ fn test_derivative_fd<F,K: Kernel<f64>>(x: Point3<f64>, compare: F, order: usize
}
}
fn test_hrbf_derivative_simple<K: Kernel<f64>>(order: usize) {
fn test_hrbf_derivative_simple<K: Kernel<f64> + 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>>(order: usize) {
for &x in test_pts.iter() { test_derivative_fd::<_,K>(x, rel_compare, order); }
}
fn test_hrbf_derivative_random<K: Kernel<f64>>(order: usize) {
fn test_hrbf_derivative_random<K: Kernel<f64> + Default>(order: usize) {
use self::rand::{SeedableRng, StdRng};
use self::rand::distributions::{IndependentSample, Range};
......@@ -128,7 +127,7 @@ fn test_hrbf_derivative_random<K: Kernel<f64>>(order: usize) {
}
}
fn test_hrbf_fit<K: Kernel<f64>>() {
fn test_hrbf_fit<K: Kernel<f64> + Default>() {
let (pts, nmls) = cube();
let mut hrbf = HRBF::<f64,K>::new(pts.clone());
......
......@@ -42,13 +42,13 @@ fn test_kernel<F,K: Kernel<Num>>(ker: K, x0: f64, compare: F)
}
}
fn test_kernel_simple<K: Kernel<Num>>(kern: K) {
fn test_kernel_simple<K: Kernel<Num> + Copy>(kern: K) {
for &x in [0.0, 1.0, 0.5, ::std::f64::consts::PI].iter() {
test_kernel(kern, x, ulp_compare);
}
}
fn test_kernel_random<K: Kernel<Num>>(kern: K) {
fn test_kernel_random<K: Kernel<Num> + Copy>(kern: K) {
use self::rand::{SeedableRng, StdRng};
use self::rand::distributions::{IndependentSample, Range};
......
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