Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Experiments with offset curves #230

Merged
merged 23 commits into from
Apr 1, 2023
Merged

Experiments with offset curves #230

merged 23 commits into from
Apr 1, 2023

Conversation

raphlinus
Copy link
Contributor

This branch has some experimental work on offset curves. It is not polished enough to merge, but contains substantial work that is likely to be useful for offset curve / stroke functionality in kurbo, which we want at some point. See Parallel curves of cubic Béziers blog post for the research work that went into this. There's some overlap between the Rust code here and the JavaScript code in that blog post, but also some divergence.

Here is a brief tour of what's in this branch:

  • There is a quartic solver that is quite robust and in fairly good shape (in common.rs). That said, it has no way of reporting complex roots, which we need. The most general approach would be to take a dep on num-complex, but that may not be necessary, see below about factoring into quadratics.

  • There's an implementation of curve fitting (cubic_fit) in offset.rs. This version lags the JavaScript version in some important ways: it only looks at real roots, and it doesn't have a good way of dealing with d0 or d1 becoming negative.

  • There's no sophisticated error metric between the cubic approximation and the true offset. What's done in place is just arc length, but this is inadequate for real use. The JavaScript code has a good solution based on intersecting rays normal to the true curve with the approximation, but that code is also not 100% perfect as it uses equally spaced samples on the t parametrization.

Some other limitations before this can be considered robust. It normalizes the chord to (0, 0) - (1, 0), which is not robust when the chord of the original curve is at or near zero. Probably the best way to deal with this is to not normalize at all, but work with area and moment in the original non-scaled coordinate space. There are other issues with determining tangent directions at or near a cusp - it should probably do something like apply L'Hôpital's rule (using second derivative) when the first derivative is near zero.

About complex roots: what the JavaScript code does is factor the quartic into two quadratics. I've observed that one of the quadratics always has real roots, so the solver never goes into the case of needing complex coefficients for the quadratics themselves. Then, when there are complex roots, finding the real part of a complex conjugate pair is trivial, for x^2 + c1 x + c0, it is just -c1/2.

Following Orellana and de Michele 2020.

WIP
Add a suite of tests from the paper. Also some cleanup (getting rid of println).
Use abs and copysign rather than an abs to make the trig case give the highest-value root. This is probably a minor improvement, and makes the code slightly less true to the paper.
@raphlinus
Copy link
Contributor Author

I should also add, one of the things missing here is cusp-finding. The JavaScript code includes an implementation, but I'm not happy with it to the point where I think it's up to kurbo's standards. There are some open questions. Is it really necessary to find all cusps? It might be better to just sample the curve at n points and look for sign-crossings, perhaps using ITP to find the exact cusp. This has the risk of missing a double cusp. That should be ok most of the time (you just fit a cubic as best you can), but the risk is that on subdividing you discover an additional sign crossing. I see two ways to deal with that: find the newly uncovered cusps and subdivide there; or just reverse the direction of the tangent vector when fitting.

Of course, another viable alternative is to robustly find all cusps before doing any fitting.

Strengthen the cubic fit logic, including better error measuring, and start cusp finding. However, it doesn't fully solve the parallel curve problem yet, and I have some refactoring in mind.
This commit gets a new `ParamCurveFit` trait in place, and uses it to compute offsets, but it is not breaking at cusps (which seems to be necessary, fitting through a cusp produces artifacts). Also needs some cleanup.
Use ITP method to find cusps. The logic probably doesn't handle double and higher multiplicity cusps.

Also arguably we should search for double cusps, but this is likely a fairly subtle improvement.
Also significant cleanup of `CurveOffset`, as a lot of that has moved into the fit module.
Mostly small changes to CubicOffset.
This computes an optimized fit (minimum number of segments, minimum error given that number), but without handling cusps.

Also includes an example that plots a quartic polynomial.
Make the ParamCurveFit trait return the raw integrals, for refinement into moments in the curve fitting logic. This is easier and more efficient for new curves, and also facilitates more efficient and robust downstream processing.
The moment_integrals function is probably correct but likely can be optimized.
Implementation of SimplifyBezPath including ParamCurveFit.

This also fixes some issues with the trait and fit implementation.

Also a bit of cleanup, comment fixes, and micro-optimization.
@ctrlcctrlv
Copy link

ctrlcctrlv commented Dec 30, 2022

I hope you don't mind, but I forked this branch to master of kurbo.rlib in MFEK#1 and am going to try to use it in production against a large dataset today. Well aware it's not done and API might change, but since it's work that I have critically needed and also offered several people to do, just done for me without me having to do a thing and I also don't have to pay for it,¹ I'm eager to see if it makes FRB American Cursive small enough for a non-janky distribution channel (like archive.org) because the fonts all together are 168 207 megabytes.

  • ¹ Perhaps I should start being nicer to @davelab6, as he's just paid me again in a roundabout way it seems. 😂

@raphlinus
Copy link
Contributor Author

Let me know what you find. I've done my best to make it robust, but haven't done any kind of large scale testing of it. And yes, I would like to land this in kurbo proper, but would like to resolve some of the things marked as "discussion question" or "TODO" in the current code.

Comment on lines +80 to +83
/// Set up a new Bézier path for simplification.
///
/// Currently this is not dealing with discontinuities at all, but it
/// could be extended to do so.

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I tried to work around this issue by implementing a ContinuousBeziersIterator which can split up a piecewise Bézier path into all of its continuous subpaths, and that work is in MFEK/math.rlib#30. The reason I did it in MFEK/math.rlib and not as a PR to this is because I don't want to step on roll over your toes while you're actively working on this, and because I'm just more familiar with MFEK project's internal API's than kurbo internals.

I am reasonably confident it works well for most paths and so am now at a place where I can do actual tests with it. When I'm done, I'll heavily refactor it as a lot of the code was just PoC, and try to make a nice API for kurbo, if there's interest in that.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would be interested in this.

@richard-uk1
Copy link
Collaborator

Would you be willing to land the code with a note on the docs saying that it's considered 'experimental' quality?

@raphlinus
Copy link
Contributor Author

Would you be willing to land the code with a note on the docs saying that it's considered 'experimental' quality?

Yes, I would. I'll write up the major to-do items. I hope to be able to do this today.

@richard-uk1
Copy link
Collaborator

I have some code that I'm using to test this PR. I'll post it below (it depends on a library I made to make my life easier working with vello - https://github.com/derekdreery/snog).

I really want to use the code here for something in production at the company I'm working at now, so I'm planning to really stress test it and fix edge cases as I find them. The first one will be finding and splitting at cusps.

use snog::{
    kurbo::{self, Affine, BezPath, Circle, CubicBez, Line, Point, Rect, Size},
    peniko::{Brush, Color, Fill, Stroke},
    App, ControlFlow, ElementState, Event, MouseButton, RenderCtx, Screen,
};

fn main() {
    let data = Data::new();
    App::new_with_data(data)
        .with_render(Data::render)
        .with_event_handler(Data::event_handler)
        .run()
}

#[derive(Debug, Copy, Clone)]
struct MovingPnt {
    pnt: Pnt,
    /// We need this because the user might not click exactly the center of the point.
    move_start: Point,
}

struct Data {
    curve: CubicBez,
    solid_stroke: Stroke,
    offset_stroke: Stroke,
    lines_stroke: Stroke,
    points_brush: Brush,
    lines_brush: Brush,
    curve_brush: Brush,
    global: Affine,
    global_inv: Affine,
    /// The last recorded location of the pointer.
    ptr_loc: Point,
    /// Which control/end point are we moving, if any
    moving_curve_pt: Option<MovingPnt>,
}

impl Data {
    fn new() -> Self {
        Self {
            curve: CubicBez {
                p0: Point::new(0.3, 0.4),
                p1: Point::new(0.4, 0.4),
                p2: Point::new(0.3, 0.7),
                p3: Point::new(0.4, 0.7),
            },
            solid_stroke: Stroke::new(2.),
            offset_stroke: Stroke::new(1.),
            lines_stroke: Stroke::new(2.).with_dashes(0., [5., 5.]),
            points_brush: Color::GREEN.into(),
            lines_brush: Color::GRAY.into(),
            curve_brush: Color::WHITE.into(),
            global: Affine::IDENTITY,
            global_inv: Affine::IDENTITY,
            ptr_loc: Point::ZERO,
            moving_curve_pt: None,
        }
    }

    fn set_global(&mut self, global: Affine) {
        self.global_inv = global.inverse();
        self.global = global;
    }

    fn event_handler(&mut self, event: Event, control_flow: &mut ControlFlow) {
        match event {
            Event::CloseRequested => *control_flow = ControlFlow::ExitWithCode(0),
            Event::CursorMoved { pos } => {
                self.ptr_loc = pos;
                self.cursor_moved();
            }
            Event::MouseInput {
                state: ElementState::Pressed,
                button: MouseButton::Left,
            } => self.left_btn_down(),
            Event::MouseInput {
                state: ElementState::Released,
                button: MouseButton::Left,
            } => self.left_btn_up(),
            Event::Resized { screen } => {
                println!("{screen:?}");
                self.on_resize(screen);
            }
            _ => (),
        }
    }

    /// The cursor moved
    ///
    /// `self.ptr_loc` is the new location of the cursor.
    fn cursor_moved(&mut self) {
        // nothing yet
    }

    /// The user pressed the left mouse button
    fn left_btn_down(&mut self) {
        if let Some(pnt) = self.nearest_point(self.ptr_loc) {
            self.moving_curve_pt = Some(MovingPnt {
                pnt,
                move_start: self.ptr_loc,
            });
        }
    }

    /// The user released the left mouse button
    fn left_btn_up(&mut self) {
        if let Some(mv) = self.moving_curve_pt {
            let offset = self.global_inv * self.ptr_loc - self.global_inv * mv.move_start;
            match mv.pnt {
                Pnt::P0 => self.curve.p0 += offset,
                Pnt::P1 => self.curve.p1 += offset,
                Pnt::P2 => self.curve.p2 += offset,
                Pnt::P3 => self.curve.p3 += offset,
            }
        }
        self.moving_curve_pt = None;
    }

    /// Find the point nearest to the given point, or `None` if not near any points
    fn nearest_point(&self, to: Point) -> Option<Pnt> {
        const CLOSE_DIST: f64 = 20.;
        self.points()
            .into_iter()
            // reverse so that the point drawn last is selected first (it will be on top)
            .rev()
            .map(|(pnt, point)| (pnt, (self.global * point - to).hypot()))
            .filter(|(_, dist)| *dist < CLOSE_DIST)
            .fold((None, f64::INFINITY), |(cur_pt, cur_dist), (pt, dist)| {
                if cur_dist > dist {
                    (Some(pt), dist)
                } else {
                    (cur_pt, cur_dist)
                }
            })
            .0
    }

    fn points(&self) -> [(Pnt, Point); 4] {
        [
            (Pnt::P0, self.curve.p0),
            (Pnt::P1, self.curve.p1),
            (Pnt::P2, self.curve.p2),
            (Pnt::P3, self.curve.p3),
        ]
    }

    /// Get the current curve, taking into account any points that are being moved.
    fn curve(&self) -> CubicBez {
        let mut curve = self.curve;
        if let Some(mv) = self.moving_curve_pt {
            let offset = self.global_inv * self.ptr_loc - self.global_inv * mv.move_start;
            match mv.pnt {
                Pnt::P0 => curve.p0 += offset,
                Pnt::P1 => curve.p1 += offset,
                Pnt::P2 => curve.p2 += offset,
                Pnt::P3 => curve.p3 += offset,
            }
        }
        curve
    }

    fn on_resize(&mut self, screen: Screen) {
        let Size { width, height } = screen.size();
        let offset = if width > height {
            ((width - height) * 0.5, 0.)
        } else {
            (0., (height - width) * 0.5)
        };
        let size = width.min(height);
        let viewport = Rect::from_origin_size(offset, Size::new(size, size));
        self.set_global(Affine::map_unit_square(viewport));
    }

    fn render(&mut self, mut ctx: RenderCtx) {
        self.draw_scene(&mut ctx);
    }

    fn draw_scene(&self, ctx: &mut RenderCtx) {
        let curve = self.curve();
        let offset_left = offset_path(curve, 0.03);
        let offset_right = offset_path(curve, -0.03);
        self.draw_line(
            ctx,
            Line {
                p0: curve.p0,
                p1: curve.p1,
            },
        );
        self.draw_line(
            ctx,
            Line {
                p0: curve.p3,
                p1: curve.p2,
            },
        );
        ctx.stroke(
            &self.solid_stroke,
            Affine::IDENTITY,
            &self.curve_brush,
            None,
            &(self.global * curve),
        );
        ctx.stroke(
            &self.offset_stroke,
            Affine::IDENTITY,
            &Color::RED,
            None,
            &(self.global * offset_left),
        );
        ctx.stroke(
            &self.offset_stroke,
            Affine::IDENTITY,
            &Color::GREEN,
            None,
            &(self.global * offset_right),
        );
        self.draw_point(ctx, curve.p0, false);
        self.draw_point(ctx, curve.p1, true);
        self.draw_point(ctx, curve.p2, true);
        self.draw_point(ctx, curve.p3, false);
    }

    fn draw_point(&self, ctx: &mut RenderCtx, point: Point, square: bool) {
        let size = 10.;
        if square {
            let shape = Rect::from_center_size(self.global * point, Size::new(size, size));
            ctx.fill(
                Fill::NonZero,
                Affine::IDENTITY,
                &self.points_brush,
                None,
                &shape,
            );
        } else {
            let shape = Circle::new(self.global * point, size * 0.5);
            ctx.fill(
                Fill::NonZero,
                Affine::IDENTITY,
                &self.points_brush,
                None,
                &shape,
            );
        };
    }

    fn draw_line(&self, ctx: &mut RenderCtx, line: Line) {
        ctx.stroke(
            &self.lines_stroke,
            Affine::IDENTITY,
            &self.lines_brush,
            None,
            &(self.global * line),
        );
    }
}

#[derive(Debug, Copy, Clone)]
enum Pnt {
    P0,
    P1,
    P2,
    P3,
}

fn offset_path(c: CubicBez, d: f64) -> BezPath {
    let co = kurbo::offset::CubicOffset::new(c, d);
    kurbo::fit_to_bezpath_opt(&co, 1e-3)
}

Cleans up warnings and clippy lints. Also makes the copyright strings consistent and adds docstrings to each source file.
Use core/alloc types, add shims to libm as needed.
Add a bunch of explanation, including how to use, how the algorithm works, and limitations that we hope to address.
@raphlinus raphlinus mentioned this pull request Mar 29, 2023
@raphlinus raphlinus marked this pull request as ready for review March 29, 2023 21:02
Copy link
Collaborator

@richard-uk1 richard-uk1 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Really excited to land this and start playing with it.

// SPDX-License-Identifier: Apache-2.0 OR MIT

//! A test case for the cubic Bézier path simplifier, generating a very
//! efficient and accurage SVG file to plot a quartic polynomial. This
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

typo (accurage -> accurate)

use kurbo::{offset::CubicOffset, CubicBez, Shape};

fn main() {
println!("<svg width='800' height='600' xmlns='http://www.w3.org/2000/svg'>");
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the future, we could have interactive examples using Xilem.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Having robust intersection logic in kurbo in the future would be awesome.

use kurbo::{BezPath, Point};

fn plot_fn(f: &dyn Fn(f64) -> f64, d: &dyn Fn(f64) -> f64, xa: f64, xb: f64, n: usize) -> BezPath {
let width = 800.;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A lot of these examples use a lot of hard-coded numbers to output a specific SVG image. I think this is fine, but it would be good to have examples that either work on a standard size (say 100x100) or are parameterized by size, at some point in the future.

result.push(root);
}
return result;
return solve_quadratic(c0, c1, c2).iter().copied().collect();
Copy link
Collaborator

@richard-uk1 richard-uk1 Mar 30, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I feel like ArrayVec should make it easier to expand size (perhaps with a From impl, but this would require arithmetic trait bounds like impl<N1, N2> From<ArrayVec<N2>> for ArrayVec<N1> where N2 < N1).

Obviously not relevant to this PR, just general musings.

/// Compute epsilon relative to coefficient.
///
/// A helper function from the Orellana and De Michele paper.
fn eps_rel(raw: f64, a: f64) -> f64 {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Interesting that this function has a discontinuity at 0. If a is very close to zero, the number will be massive, but at 0 it goes back to raw. I assume this does that job that the paper authors require as-is.

// TODO: handle case d_2 is very small?
} else {
// I think this case means no real roots, return empty.
todo!()
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you want to leave this as a panic or change it?

I've never hit this while planing with the cubic offset.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Update, I have hit it. I've included the specific curves/offset/accuracy that hit this todo.

src/fit.rs Outdated
/// than one such discontinuity, any can be reported, as the function will
/// be called repeatedly after subdivision of the range.
///
/// Do not report cusps at the endpoints of the range. (TODO: talk about
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There's a TODO here, I'd be inclined to leave it unless it's easy to do.

Comment on lines +80 to +83
/// Set up a new Bézier path for simplification.
///
/// Currently this is not dealing with discontinuities at all, but it
/// could be extended to do so.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would be interested in this.

@richard-uk1
Copy link
Collaborator

richard-uk1 commented Mar 30, 2023

I managed to hit the todo at common.rs:454 with the following curve

CubicBez {
    p0: (0.36071499208860763, 0.4284909018987342),
    p1: (0.40765427215189876, 0.2682357594936709),
    p2: (0.3130191851265823, 0.6503065664556962),
    p3: (0.27402590981012653, 0.7838162579113924),
}

With the offset being one of 0.02 or -0.02 (not sure which). Tolerance for fitting was 1e-3.

This shouldn't be a block to merging, but once you have I'll create a ticket with this example in.


You can get some incorrect offsets by setting the control points outside the end points, with all points co-linear (the above curve is an example of this). These are pathological curves, and for my use cases, I'm happy to give garbage results when people use garbage curves.


Ooh, another interesting one. The following curve (with the same offset/tol as before) caused a stack overflow.

CubicBez {
    p0: (0.20421023447401776, 0.6660408745247148),
    p1: (0.4, 0.4),
    p2: (0.20270516476552597, 0.42017585551330794),
    p3: (0.39073193916349824, 0.7806846086818757),
}

///
/// When a higher degree of optimization is desired (at considerably more runtime cost),
/// consider [`fit_to_bezpath_opt`] instead.
pub fn fit_to_bezpath(source: &impl ParamCurveFit, accuracy: f64) -> BezPath {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It might be useful to say what accuracy is measuring. Is it the maximum pointwise distance between the two curves, or is it the area between the curves, or something else?

///
/// Returns the cubic segment and the square of the error.
/// Discussion question: should this be a method on the trait instead?
pub fn fit_to_cubic(
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same comment as fit_to_bezpath re. accuracy parameter

/// optimal subdivision points. Its result is expected to be very close to the optimum
/// possible Bézier path for the source curve, in that it has a minimal number of curve
/// segments, and a minimal error over all paths with that number of segments.
pub fn fit_to_bezpath_opt(source: &impl ParamCurveFit, accuracy: f64) -> BezPath {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

same comment as fit_to_bezpath re. accuracy parameter

/// optimal subdivision points. Its result is expected to be very close to the optimum
/// possible Bézier path for the source curve, in that it has a minimal number of curve
/// segments, and a minimal error over all paths with that number of segments.
pub fn fit_to_bezpath_opt(source: &impl ParamCurveFit, accuracy: f64) -> BezPath {
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this still fails when the offset path has a cusp (original path curvature = offset). Not sure if you want to fix it before or after landing but I thought I'd mention.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, this is expected, and could probably be explained better. The "source curve" here is actually the offset, so when that has a discontinuity it won't work well. It's easy for people to imagine the "source curve" being the original cubic. I'll try to write a bit more clarification before landing this, then hope to do more followup.

@richard-uk1
Copy link
Collaborator

Another comment (I cannot remember exactly where in the code it applies to):

You mention about whether to handle parts of the curve that are not continuously differentiable or to leave them to client code. How would leaving this to client code work? Would you provide a list of all the discontinuous points, and the derivatives at those points?

@raphlinus
Copy link
Contributor Author

Thanks, it's genuinely useful to have cases where the root-finding fails. It's not surprising it's cases where the points are co-linear. Even so, the algorithm should be made robust against those, as one of the goals here is to expand stroked paths to outlines, and, even though those cases are somewhat pathological, it is possible to characterize what the output should look like (and it definitely shouldn't panic!).

Mostly just doc improvements.
@raphlinus
Copy link
Contributor Author

Ok, I responded to the doc comments. I agree the examples can be made much more complete; this was mostly to help me develop the algorithm.

Regarding simplification of Bézier paths with corners, this is a deeper philosophical question: should kurbo provide core geometric primitives that can be used by a sophisticated client, or more "user-friendly" methods such as "simplify this path"? In the former case, I imagine that the client will break the source path into continuous path ranges itself, decide for each range what to do (some might be straight lines, not needing simplification, or a single curved segment), and call the provided function. In the latter case, the fit_to_bezpath function is too low level.

Probably the best way to go forward is to implement a path->path simplify function in terms of the existing fit_to_bezpath and existing SimplifyCubic implementations, in particular, not reporting cusps through break_cusp. The sophisticated client would still be able to call the primitive.

@raphlinus raphlinus merged commit d18e302 into master Apr 1, 2023
@raphlinus raphlinus deleted the offset_experiment branch April 1, 2023 18:59
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants