Skip to content

Commit

Permalink
Replaced Vec based Polynomials with Matrix based ones
Browse files Browse the repository at this point in the history
  • Loading branch information
@gammelalf
gammelalf committed Jun 10, 2022
1 parent f459fa9 commit f256396
Show file tree
Hide file tree
Showing 5 changed files with 299 additions and 71 deletions.
84 changes: 36 additions & 48 deletions src/bezier.rs
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
use std::ops::{Deref, DerefMut, Add};
use nalgebra::{RealField, Scalar, Field, Vector2, RowVector2, RowVector3, RowVector4, RowDVector};
use nalgebra::{RealField, Scalar, Field, Vector2, RowVector2, RowVector3, RowVector4, RowDVector, Matrix2xX};
use num::{Num, One};
use smallvec::{SmallVec, smallvec};
use crate::graham_scan::convex_hull;
use crate::bounding_box::BoundingBox;
use crate::polynomial::Polynomial;
use crate::npolynomial::{Polynomial, Polynomial1xX, Polynomial2xX};

#[derive(Clone, Debug, PartialEq)]
pub struct BezierCurve<K: Scalar>(pub CurveInternal<K>);
Expand Down Expand Up @@ -66,9 +66,11 @@ impl <K: RealField + Scalar> BezierCurve<K> {
let mut points: SmallVec<[Vector2<K>; 6]> = SmallVec::new();
points.push(self[0].clone());
points.push(self[self.len()-1].clone());
let [dx, dy] = self.derivative();
dx.roots().into_iter()
.chain(dy.roots().into_iter())
let d = self.derivative();
let roots_x = Polynomial(d.0.row(0)).roots();
let roots_y = Polynomial(d.0.row(1)).roots();
roots_x.into_iter()
.chain(roots_y.into_iter())
.filter(|t| &K::zero() <= t && t <= &K::one())
.for_each(|t| {
points.push(self.castlejau_eval(t));
Expand Down Expand Up @@ -303,28 +305,26 @@ impl <K: Field + Scalar> BezierCurve<K> {

/* Stuff using polynomials */
impl <K: Field + Scalar> BezierCurve<K> {
pub fn polynomial(&self) -> [Polynomial<K>; 2] {
pub fn polynomial(&self) -> Polynomial2xX<K> {
let zero = K::zero();
let one = K::one();
match &self[..] {
[a] => {
[0, 1].map(|i| Polynomial(vec![a[i].clone()]))
Polynomial(Matrix2xX::from_iterator(1, a.into_iter().map(Clone::clone)))
}
[a, b] => {
let p_a = a * RowVector2::new(one.clone(), zero.clone() - one.clone());
let p_b = b * RowVector2::new(zero, one);
let p = p_a + p_b;
[Polynomial(vec![p[(0, 0)].clone(), p[(0, 1)].clone()]),
Polynomial(vec![p[(1, 0)].clone(), p[(1, 1)].clone()])]
Polynomial(Matrix2xX::from_iterator(2, p.into_iter().map(Clone::clone)))
}
[a, b, c] => {
let two = one.clone() + one.clone();
let p_a = a * RowVector3::new(one.clone(), zero.clone() - two.clone(), one.clone());
let p_b = b * RowVector3::new(zero.clone(), two.clone(), zero.clone() - two);
let p_c = c * RowVector3::new(zero.clone(), zero, one);
let p = p_a + p_b + p_c;
[Polynomial(vec![p[(0, 0)].clone(), p[(0, 1)].clone(), p[(0, 2)].clone()]),
Polynomial(vec![p[(1, 0)].clone(), p[(1, 1)].clone(), p[(1, 2)].clone()])]
Polynomial(Matrix2xX::from_iterator(3, p.into_iter().map(Clone::clone)))
}
[a, b, c, d] => {
let three = one.clone() + one.clone() + one.clone();
Expand All @@ -334,49 +334,44 @@ impl <K: Field + Scalar> BezierCurve<K> {
let p_c = c * RowVector4::new(zero.clone(), zero.clone(), three.clone(), zero.clone() - three);
let p_d = d * RowVector4::new(zero.clone(), zero.clone(), zero, one);
let p = p_a + p_b + p_c + p_d;
[Polynomial(vec![p[(0, 0)].clone(), p[(0, 1)].clone(), p[(0, 2)].clone(), p[(0, 3)].clone()]),
Polynomial(vec![p[(1, 0)].clone(), p[(1, 1)].clone(), p[(1, 2)].clone(), p[(1, 3)].clone()])]
Polynomial(Matrix2xX::from_iterator(4, p.into_iter().map(Clone::clone)))
}
_ => {
let mut ps = bernstein_polynomials::<K>(self.degree())
.into_iter()
.zip(self.iter())
.map(|(p, a)| {
let p = RowDVector::from_iterator(self.len(), p.0.into_iter()); // TODO chain with 0 else panic
a * p
a * p.0
});
let p = if let Some(mut p) = ps.next() {
if let Some(mut p) = ps.next() {
for q in ps {
p += q;
}
p
return Polynomial(p);
} else {
unreachable!();
};
[0, 1].map(|i| Polynomial(p.row(i).iter().map(Clone::clone).collect()))
}
}
}

pub fn derivative(&self) -> [Polynomial<K>; 2] {
pub fn derivative(&self) -> Polynomial2xX<K> {
let zero = K::zero();
let one = K::one();
match &self[..] {
[_] => {
[Polynomial(vec![]), Polynomial(vec![])]
Polynomial(Matrix2xX::zeros(0))
}
[a, b] => {
let p = b - a;
[Polynomial(vec![p[(0, 0)].clone()]),
Polynomial(vec![p[(1, 0)].clone()])]
Polynomial(Matrix2xX::from_iterator(1, p.into_iter().map(Clone::clone)))
}
[a, b, c] => {
let two = one.clone() + one.clone();
let p_a = (b - a) * RowVector2::new(one.clone(), zero.clone() - one.clone());
let p_b = (c - b) * RowVector2::new(zero, one);
let p = (p_a + p_b) * two;
[Polynomial(vec![p[(0, 0)].clone(), p[(0, 1)].clone()]),
Polynomial(vec![p[(1, 0)].clone(), p[(1, 1)].clone()])]
Polynomial(Matrix2xX::from_iterator(2, p.into_iter().map(Clone::clone)))
},
[a, b, c, d] => {
let two = one.clone() + one.clone();
Expand All @@ -385,37 +380,30 @@ impl <K: Field + Scalar> BezierCurve<K> {
let p_b = (c - b) * RowVector3::new(zero.clone(), two.clone(), zero.clone() - two);
let p_c = (d - c) * RowVector3::new(zero.clone(), zero, one);
let p = (p_a + p_b + p_c) * three;
[Polynomial(vec![p[(0, 0)].clone(), p[(0, 1)].clone(), p[(0, 2)].clone()]),
Polynomial(vec![p[(1, 0)].clone(), p[(1, 1)].clone(), p[(1, 2)].clone()])]
Polynomial(Matrix2xX::from_iterator(3, p.into_iter().map(Clone::clone)))
}
_ => {
let mut degree = zero;
let mut ps = bernstein_polynomials::<K>(self.degree()-1)
.into_iter()
.enumerate()
.map(|(i, p)| {
degree = degree.clone() + one.clone(); // hack converting degree from usize to K

let point = &self[i+1] - &self[i];
let p = RowDVector::from_iterator(self.len() - 1, p.0.into_iter()); // TODO chain with 0 else panic
point * p
let a = &self[i+1] - &self[i];
a * p.0
});
let p = if let Some(mut p) = ps.next() {
if let Some(mut p) = ps.next() {
for q in ps {
p += q;
}
p
return Polynomial(p);
} else {
unreachable!();
};
[0, 1].map(|i| Polynomial(p.row(i).iter().map(Clone::clone).collect()))
}
}
}
}

pub fn tangent(&self, t: K) -> Vector2<K> {
let [dx, dy] = self.derivative();
Vector2::new(dx.evaluate(t.clone()), dy.evaluate(t.clone()))
self.derivative().evaluate(t.clone())
}

pub fn normal(&self, t: K) -> Vector2<K> {
Expand All @@ -427,8 +415,8 @@ impl <K: Field + Scalar> BezierCurve<K> {
}
}

pub fn bernstein_polynomials<N>(degree: usize) -> Vec<Polynomial<N>>
where N: Num + Clone
pub fn bernstein_polynomials<N>(degree: usize) -> Vec<Polynomial1xX<N>>
where N: Field + Scalar
{
let one = N::one();
let zero = N::zero();
Expand All @@ -439,22 +427,22 @@ pub fn bernstein_polynomials<N>(degree: usize) -> Vec<Polynomial<N>>
Vec::with_capacity(degree+1),
);

powers.0.push(Polynomial(vec![one.clone()])); // TODO don't alloc p(x) = 1
powers.1.push(Polynomial(vec![one.clone()]));
powers.0.push(Polynomial(vec![zero.clone(), one.clone()]));
powers.1.push(Polynomial(vec![one.clone(), zero - one]));
powers.0.push(Polynomial(RowDVector::from_row_slice(&[one.clone()]))); // TODO don't alloc p(x) = 1
powers.1.push(Polynomial(RowDVector::from_row_slice(&[one.clone()])));
powers.0.push(Polynomial(RowDVector::from_row_slice(&[zero.clone(), one.clone()])));
powers.1.push(Polynomial(RowDVector::from_row_slice(&[one.clone(), zero - one])));

for i in 1..degree {
powers.0.push(&powers.0[i] * &powers.0[1]);
powers.1.push(&powers.1[i] * &powers.1[1]);
powers.0.push(powers.0[i].mul(&powers.0[1].0));
powers.1.push(powers.1[i].mul(&powers.1[1].0));
}

// Combine powers into Bernstein polynomials
let mut base = Vec::with_capacity(degree+1);
let pascal = pascal_triangle::<N>(degree);
for (i, coeff) in pascal.into_iter().enumerate() {
let mut b = &powers.0[i] * &powers.1[degree-i];
b *= coeff; // Use MutAssign to multiply inplace and avoid reallocation
let mut b = powers.0[i].mul(&powers.1[degree-i].0);
b.0 *= coeff; // Use MutAssign to multiply inplace and avoid reallocation
base.push(b);
}

Expand Down
42 changes: 21 additions & 21 deletions src/lib.rs
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,15 @@
pub mod bounding_box;
pub mod polynomial;
pub mod npolynomial;
pub mod graham_scan;
pub mod bezier;

#[cfg(test)]
mod tests {
use smallvec::smallvec;
use nalgebra::Vector2;
use nalgebra::{RowDVector, Vector2, RowVector2, RowVector3};
use crate::bezier::{bernstein_polynomials, BezierCurve, pascal_triangle};
//use crate::graham_scan;
use crate::polynomial::Polynomial;
use crate::npolynomial::Polynomial;

#[test]
fn bezier_split() {
Expand Down Expand Up @@ -46,8 +46,7 @@ mod tests {
Vector2::new(0.0, 66.0),
Vector2::new(50.0, 100.0),
]);
assert_eq!(curve.x_derivative(), curve.x_polynomial().derive());
assert_eq!(curve.y_derivative(), curve.y_polynomial().derive());
assert_eq!(curve.derivative(), curve.polynomial().derive());
}

/*#[test]
Expand Down Expand Up @@ -85,7 +84,8 @@ mod tests {

#[test]
fn polynomials() {
assert_eq!(Polynomial(vec![1.0, 2.0, 3.0]).derive(), Polynomial(vec![2.0, 6.0]));
assert_eq!(Polynomial(RowVector3::new(1.0, 2.0, 3.0)).derive(),
Polynomial(RowVector2::new(2.0, 6.0)));

let cmp_float = |x: &f64, y: &f64| x.partial_cmp(y).expect("A wild NaN appeared");

Expand All @@ -97,7 +97,7 @@ mod tests {
(2.0, 3.0),
(4.0, 5.0),
] {
let p = &Polynomial(vec![-n, 1.0]) * &Polynomial(vec![-m, 1.0]);
let p = Polynomial(RowVector2::new(-n, 1.0)).mul(&RowVector2::new(-m, 1.0));
let mut roots = p.roots(); roots.sort_by(cmp_float);
assert_eq!(roots, vec![n, m])
}
Expand All @@ -115,23 +115,23 @@ mod tests {

#[test]
fn bernstein() {
assert_eq!(bernstein_polynomials::<i32>(0), vec![
Polynomial(vec![1]),
assert_eq!(bernstein_polynomials::<f32>(0), vec![
Polynomial(RowDVector::from_row_slice(&[1.0])),
]);
assert_eq!(bernstein_polynomials::<i32>(1), vec![
Polynomial(vec![1, -1]),
Polynomial(vec![0, 1]),
assert_eq!(bernstein_polynomials::<f32>(1), vec![
Polynomial(RowDVector::from_row_slice(&[1.0, -1.0])),
Polynomial(RowDVector::from_row_slice(&[0.0, 1.0])),
]);
assert_eq!(bernstein_polynomials::<i32>(2), vec![
Polynomial(vec![1, -2, 1]),
Polynomial(vec![0, 2, -2]),
Polynomial(vec![0, 0, 1]),
assert_eq!(bernstein_polynomials::<f32>(2), vec![
Polynomial(RowDVector::from_row_slice(&[1.0, -2.0, 1.0])),
Polynomial(RowDVector::from_row_slice(&[0.0, 2.0, -2.0])),
Polynomial(RowDVector::from_row_slice(&[0.0, 0.0, 1.0])),
]);
assert_eq!(bernstein_polynomials::<i32>(3), vec![
Polynomial(vec![1, -3, 3, -1]),
Polynomial(vec![0, 3, -6, 3]),
Polynomial(vec![0, 0, 3, -3]),
Polynomial(vec![0, 0, 0, 1]),
assert_eq!(bernstein_polynomials::<f32>(3), vec![
Polynomial(RowDVector::from_row_slice(&[1.0, -3.0, 3.0, -1.0])),
Polynomial(RowDVector::from_row_slice(&[0.0, 3.0, -6.0, 3.0])),
Polynomial(RowDVector::from_row_slice(&[0.0, 0.0, 3.0, -3.0])),
Polynomial(RowDVector::from_row_slice(&[0.0, 0.0, 0.0, 1.0])),
]);
}
}
2 changes: 1 addition & 1 deletion src/main.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@ use crate::svg::SVG;

mod bounding_box;
mod graham_scan;
mod polynomial;
mod npolynomial;
mod bezier;
mod svg;

Expand Down
Loading

0 comments on commit f256396

Please sign in to comment.