Experiments with offset curves

[raphlinus]

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.

[raphlinus]

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.

[ctrlcctrlv]

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]

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.

[ctrlcctrlv]

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.

[derekdreery]

I would be interested in this.

[derekdreery]

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

[raphlinus]

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.

[derekdreery]

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)
}

[derekdreery]

Really excited to land this and start playing with it.

[derekdreery]

typo (accurage -> accurate)

[derekdreery]

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

[derekdreery]

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

[derekdreery]

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.

[derekdreery]

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.

[derekdreery]

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.

[derekdreery]

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

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

[derekdreery]

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

[derekdreery]

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

[derekdreery]

I would be interested in this.

[derekdreery]

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),
}

[derekdreery]

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?

[derekdreery]

same comment as fit_to_bezpath re. accuracy parameter

[derekdreery]

same comment as fit_to_bezpath re. accuracy parameter

[derekdreery]

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.

[raphlinus]

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.

[derekdreery]

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]

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!).

[raphlinus]

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.