~raph/bezoid

b09aa211ed8a9b8cae2d11d06a334abc4ace84ab — Raph Levien 1 year, 6 months ago 4645edd
Use sophisticated numerical integration

Use Gauss-Legendre quadrature to characterize the curve, in favor of
crude Euler integration.

Rendering is still done by brute force.
3 files changed, 79 insertions(+), 17 deletions(-)

M src/bez.rs
M src/bezoid.rs
M src/grapher.rs
M src/bez.rs => src/bez.rs +2 -2
@@ 93,8 93,8 @@ impl Widget<AppData> for Bez {
        let c = CubicBez::new(Point::ORIGIN, p1, p2, Point::new(1.0, 0.0));
        ctx.stroke(a2 * c, &Color::rgb8(64, 64, 128), 1.0);
        let params: CurveParams = data.into();
        let curve = params.compute();
        grapher::plot_xy(ctx, &curve.pts, Point::new(OFFSET_X, OFFSET_Y), SCALE);
        let pts = params.render_to_pts();
        grapher::plot_xy(ctx, &pts, Point::new(OFFSET_X, OFFSET_Y), SCALE);
    }
}


M src/bezoid.rs => src/bezoid.rs +62 -3
@@ 1,6 1,7 @@
//! The math for the bezoid curve family.

use druid::kurbo::{Affine, CubicBez, Point, Vec2};
use druid::kurbo::common as coeffs;

use crate::appdata::AppData;



@@ 13,7 14,6 @@ pub struct CurveParams {
}

pub struct CurveResult {
    pub pts: Vec<Point>,
    pub th0: f64,
    pub th1: f64,
    pub chord: f64,


@@ 25,6 25,7 @@ pub const BEZ_CHORD_EXP: f64 = 1.4;
const N: usize = 1000;
const N_RECIP: f64 = 1.0 / (N as f64);

#[allow(unused)]
fn compute_basis(bias: f64) -> Vec<f64> {
    let mut sum = 0.0;
    let mut result = (0..=N)


@@ 81,6 82,13 @@ impl CurveParams {
    /// We won't keep this because we're going to do another stage of numerical
    /// integration, but it's useful as we develop.
    pub fn compute_thetas(&self) -> Vec<f64> {
        /*
        let i24 = self.integrate(0.0, 1.0, 24);
        for order in &[3, 5, 7, 9, 11] {
            let err = (self.integrate(0.0, 1.0, *order) - i24).hypot();
            println!("error at {}: {}", order, err);
        }
        */
        (0..=N)
            .map(|i| {
                let t = (i as f64) * N_RECIP;


@@ 94,8 102,59 @@ impl CurveParams {
    }

    pub fn compute(&self) -> CurveResult {
        let integral = self.integrate(0.0, 1.0, 24);
        let th_chord = integral.atan2();
        let chord = integral.hypot();
        let th0 = th_chord - self.compute_theta(0.0);
        let th1 = self.compute_theta(1.0) - th_chord;
        CurveResult {
            th0,
            th1,
            chord,
        }
    }

    pub fn integrate(&self, t0: f64, t1: f64, order: usize) -> Vec2 {
        let c = match order {
            3 => coeffs::GAUSS_LEGENDRE_COEFFS_3,
            5 => coeffs::GAUSS_LEGENDRE_COEFFS_5,
            7 => coeffs::GAUSS_LEGENDRE_COEFFS_7,
            9 => coeffs::GAUSS_LEGENDRE_COEFFS_9,
            11 => coeffs::GAUSS_LEGENDRE_COEFFS_11,
            24 => coeffs::GAUSS_LEGENDRE_COEFFS_24,
            _ => panic!("don't have coefficients for {}", order),
        };
        let mut result = Vec2::ZERO;
        let tm = 0.5 * (t1 + t0);
        let dt = 0.5 * (t1 - t0);
        for (wi, xi) in c {
            let t = tm + dt * xi;
            let th = self.compute_theta(t);
            result += *wi * Vec2::from_angle(th);
        }
        dt * result
    }

    pub fn render_to_pts(&self) -> Vec<Point> {
        let thetas = self.compute_thetas();
        integrate_curve(&thetas)
        let n = thetas.len();
        let mut p = Point::ORIGIN;
        let scale = 1.0 / (n - 1) as f64;
        let mut pts = (0..n)
            .map(|i| {
                let this_p = p;
                if i < thetas.len() - 1 {
                    let th = 0.5 * (thetas[i] + thetas[i + 1]);
                    p += scale * Vec2::from_angle(th);
                }
                this_p
            })
            .collect::<Vec<_>>();
        let a = Affine::new([p.x, p.y, -p.y, p.x, 0.0, 0.0]).inverse();
        for p in &mut pts {
            *p = a * *p;
        }
        pts
    }
}



@@ 117,12 176,12 @@ pub fn integrate_curve(thetas: &[f64]) -> CurveResult {
    for p in &mut pts {
        *p = a * *p;
    }
    println!("integral: {:?}", p.to_vec2());
    let th_chord = p.to_vec2().atan2();
    let chord = p.to_vec2().hypot();
    let th0 = th_chord - thetas[0];
    let th1 = thetas[n - 1] - th_chord;
    CurveResult {
        pts,
        th0,
        th1,
        chord,

M src/grapher.rs => src/grapher.rs +15 -12
@@ 71,23 71,26 @@ impl Widget<AppData> for Grapher {
        let thetas = params.compute_thetas();
        plot(ctx, &thetas, Point::new(50.0, 400.0), 600.0, 100.0);
        let curve = bezoid::integrate_curve(&thetas);
        plot_xy(ctx, &curve.pts, Point::new(100.0, 300.0), 400.0);
        let pts = params.render_to_pts();
        plot_xy(ctx, &pts, Point::new(100.0, 300.0), 400.0);
        let cb = curve.infer_bezier(data);
        let bez = cb.into_path(1e-3);
        let a = Affine::translate((100.0, 300.0)) * Affine::FLIP_Y * Affine::scale(400.0);
        let bez = a * bez;
        ctx.stroke(bez, &Color::rgb8(128, 128, 255), 1.0);

        let a2 = Affine::translate(Vec2::new(800.0, 250.0)) * Affine::FLIP_Y * Affine::scale(100.0);
        let solver = crate::solve::Solver;
        let p1 = Point::new(data.px1, data.py1);
        let p2 = Point::new(data.px2, data.py2);
        let l1 = Line::new(Point::ORIGIN, p1);
        ctx.stroke(a2 * l1, &Color::WHITE, 1.0);
        let l2 = Line::new(Point::new(1.0, 0.0), p2);
        ctx.stroke(a2 * l2, &Color::WHITE, 1.0);
        let params = solver.solve(p1, p2);
        let curve = params.compute();
        plot_xy(ctx, &curve.pts, a2 * Point::ORIGIN, 100.0);
        if false {
            let a2 = Affine::translate(Vec2::new(800.0, 250.0)) * Affine::FLIP_Y * Affine::scale(100.0);
            let solver = crate::solve::Solver;
            let p1 = Point::new(data.px1, data.py1);
            let p2 = Point::new(data.px2, data.py2);
            let l1 = Line::new(Point::ORIGIN, p1);
            ctx.stroke(a2 * l1, &Color::WHITE, 1.0);
            let l2 = Line::new(Point::new(1.0, 0.0), p2);
            ctx.stroke(a2 * l2, &Color::WHITE, 1.0);
            let params = solver.solve(p1, p2);
            let pts = params.render_to_pts();
            plot_xy(ctx, &pts, a2 * Point::ORIGIN, 100.0);
        }
    }
}