Implement Offset Algorithm

[thehappycheese]

• I agree to follow the project's code of conduct.
• I added an entry to CHANGES.md if knowledge of this change could be valuable to users.

---

Defines an OffsetCurve trait and implements the associated function .offset_curve(distance) for some geometry types. The offset_curve function offsets the edges of a geometry by a specified distance, where the sign of distance determines if the edges of the geometry are offset to the 'left' or 'right'.

'Offset curve' is like a 'one-sided buffer' and, unlike buffer, it is only meaningful for some geometry types. The table below summarises the different outputs.

A 'clipped offset curve' has had parts of the output removed to ensure that all edges in the output are at least distance from the input, even if the input is self-overlapping. So far, this clipping step has not been implemented.

The following table is a work in progress:

Input	Offset Curve	Offset Curve - Clipped	Buffer
Point	Undefined	Undefined	Circle (Polygon)
Line	Line	Line	Capsule (Polygon)
LineString	LineString (possibly self overlapping)	MultiLineString	Sausage (Polygon)
Polygon	Option<Polygon> (always closed, but possibly malformed)	Option<Polygon> (possibly unclosed, possibly malformed)	Polygon
...	...	...	...

GEOS and JTS

• geos::operation::buffer::OffsetCurve
• org.locationtech.jts.operation.buffer.OffsetCurveBuilder

Choice of Algorithm

• Issue Geometry buffering #641 suggests a more advanced algorithm (Voronoi-based straight skeleton approach which is meant for buffering, but might also work for offset_curve)
• Issue Support (Convex) Minkowski Sum construction #883 suggests a Minkowski Sum which is a more powerful type of buffering operation

I am not sure if the algorithm I am trying to implement here is any better or worse than those, but it is one I have got working in the past in python so I figure I have the best chance of success with it. It loosely follows the algorithm described by
Xu-Zheng Liu, Jun-Hai Yong, Guo-Qin Zheng, Jia-Guang Sun. An offset algorithm for polyline curves. Computers in Industry, Elsevier, 2007, 15p. inria-00518005

Trait Implementations

• Line<T>
• LineString<T>
• MultiLineString<T>
• Triangle<T>
• Rectangle<T>
• Polygon<T>
• MultiPolygon<T>
• Point<T>
  • MultiPoint<T>
  • Geometry<T>
  • GeometryCollection<T>

For coordinate types

• where T:CoordFloat
• where T:CoordNum (??? seems tricky)

Note:
Offset for Point Can't be defined without confusion.
The user may expect the result to be one of the following:

1. Return the input unchanged,
2. a translation (but in what direction?)
3. a circle (like a buffer, but as previously explained,
   the offset operation is different to the buffer operation).
4. an empty LineString (weird, yes, but this is actually what GEOS does 🤔)

Issues to Fix:

• Change name from Offset to OffsetCurve to match terminology in GEOS
  and jts
• Match offset sign / direction with GEOS and JTS if possible
• Make proper use of SimpleKernel / RobustKernel
• Prevent execution when the specified offset distance is zero
• Mitre-limit is implemented so that a LineString which doubles-back on
  itself will not produce an elbow at infinity
  • Make Mitre Limit configurable?
• Option to clip output such that non-adjacent line segments in the output
  are not self-intersecting.

Hopefully I have started off my contribution in the right direction? :)

Many thanks for your time

[urschrei]

Hi @thehappycheese, welcome and thanks for the contribution! I've left a load of review comments – please don't feel that you need to address them all at once. In particular, it looks like you're having some trouble using the Kernel trait, and this might be a good place to start. I and other contributors are happy to help here, and we can do it inline in the PR, and / or on Discord: https://discord.gg/Fp2aape

[urschrei]

Orientation is an enum, so you have to use a match for control flow.

something like

match Kernel::orient2d(*a, *b, *c) {
    Orientation::Collinear => { // keep going }
    _ => { // not collinear, bail out or test different combination?? }
}

As you note, the use of arbitrary "magic" threshold values without considerable justification and accompanying comments explaining the reasoning behind it for future maintainers is unlikely to be accepted (I've just finished removing an algorithm containing a lot of these that I had the misfortune of writing myself…), so if we can use existing mechanisms that would be much better.

[urschrei]

I'm a bit confused about this. Would the dot product sign, i.e. the Orientation suffice? It looks like you're just checking whether the result is -1, 1, or 0…

[thehappycheese]

Apologies, looks like my inept git-fu has disconnected all your comments from the code 😕 I think maybe re-basing the branch on my fork was the wrong thing to do.

cross_product simplifies line_segment_intersection_with_parameters(): While it is used for a collinear check at the start of that function (which could be replaced with Kernel::dot_product_sign or Kernel::orient2d) it is used twice more further along to calculate the parameters t_ab and t_cd so I do not believe it can be substituted for an orientation check.

let t_ab = cross_product_2d(ac, cd) / ab_cross_cd;
let t_cd = -cross_product_2d(ab, ac) / ab_cross_cd;

line_segment_intersection_with_parameters() is currently only used via a wrapper function line_segment_intersection_with_relationships(); in that function it checks each t_ab and t_cd to find which of three ranges they lie in:

• less than zero
• between zero and one
• greater than one

It's possible there is some tricky way to express this check in terms of orient2d(), I'll have a think about it.

It does feel weird to have cross_product() out in it's own module. In the past I have worked with libraries that have a Vector2 type where this kind of function would be implemented. It is simple enough to be inline perhaps, but I wrote up all the tests and documentation because in the past I have made many mistakes with putting negatives in the wrong place and if you make this kind of mistake an even number of times you might not realize since it cancels out 😋

[urschrei]

A minor nit: there's no need to qualify this as "very basic", since the docs will make the limitations clear.

[thehappycheese]

👍 I will try to reword it to be more specific, or be more clear about any limitations in the rest of the docs

[urschrei]

There's no need to specify this as geo operates only on coordinates in the real plane (with the exception of some explicitly documented convenience algorithms for spherical coordinates). There is an explicit assumption that only projected / cartesian coordinates are used. That this algorithm may be suitable for coordinates near the equator is a coincidence as far as the library is concerned.

[urschrei]

Is there a specific difficulty you've encountered while trying to implement a check against this?

[urschrei]

Where could this occur?

[urschrei]

We will need a check for this

[thehappycheese]

Have added this check to LineString but it probably isn't worth it for Line. I have added comments to the relevant code in my next commit.

[urschrei]

This is a fairly significant limitation, but if the algorithm can't handle it, then we should consider a check for self-crossing LineStrings (I believe we already have one) and bail out early if this is encountered.

[urschrei]

My expectation is that a positive (or outward) offset on a LineString produces a Polygon (as per JTS, GEOS, and derived libraries). Could you explain the reasoning here?

[dabreegster]

I have no context on how JTS/GEOS/etc does things, but offsetting a LineString and getting a LineString has many use cases. If you start with something representing the center of a road and want the left or right edge, or the center of a lane, for example. To buffer around a LineString and get a polygon, you could offset inwards (the hole) and outwards (the exterior), right?

[urschrei]

@dabreegster I'm struggling to visualise that. Do you have an example? The standard implementations of offset (in plain language, buffer) operations always produce a polygon from a Line or LineString as far as I know. The Shapely docs are representative here: https://shapely.readthedocs.io/en/stable/manual.html#object.buffer

[michaelkirk]

Disclaimer that I haven't looked at this implementation, but just a note on terminology...

JTS has an "offset curve" which is a line string to line string transformation.

a line lying on one side of another line at a given offset distance

http://lin-ear-th-inking.blogspot.com/2022/01/jts-offset-curves.html?m=1

[thehappycheese]

Hi @urschrei thank you so much for all your feedback I really appreciate it, I will do my best to address each one :)

Regarding the sign of the specified offset, I agree that it would be strange if a negative offset resulted in an inset for polygons. I am certainly open to changing this, I will add it to the TODO list in my PR.

On this @dabreegster has already explained my reasoning really well; To me a buffer operation different to an offset; I am interested in offsetting road center lines, not producing enclosed polygons.

I recognize that a buffer operation is more generally useful, but I also just need this signed offset function. This offset trait can be extended by a a few clipping and stitching steps to a full buffer operation, this is the end goal of the paper I am following.

[dabreegster]

I'm struggling to visualise that. Do you have an example?

https://a-b-street.github.io/docs/tech/map/geometry/index.html#projecting-a-polyline, and for the inspiration behind it, https://wwwtyro.net/2019/11/18/instanced-lines.html. The JTS offset curve looks like the same thing. I don't know what this linestring -> linestring operation should be called, but I have plenty of use cases for it, as well as buffering (to get a polygon)

[urschrei]

Thanks all, it looks like it's just confusion about terminology, mostly on my part. It might be worth discussing whether the current trait name conveys the functionality clearly enough, though maybe it's just a question of ensuring that the module-level docstring does this.

[frewsxcv]

@thehappycheese By the way, your screenshot link doesn't appear to work: https://github.com/thehappycheese/powerbi-visual-geojson-map-1/raw/main/git_docs/hero_image.jpg (404)

[thehappycheese]

Keen to align terminology with suggestions / other libraries :)

offset_curve is nice because it makes it more obvious why this type of signed offset cannot be implemented for a Point object. Point objects can only have a "buffer" or "outset" applied, not an 'offset_curve'.

Ah apologies for the broken link I picked it from one of my public repos but didn't realize the image link is to a private repo; this one should work?

This is from a custom PowerBI visual I created. I will eventually make that code public when I get a chance to clean up /republish repo. I want the PowerBI visual to dynamically connect to a rust server for a specialized kind of chainage-based geo-coding.

[thehappycheese]

I have deployed an interactive tester thing here https://geo-tester.thehappycheese.com/

Source code here: https://github.com/thehappycheese/geo-offset-test-server

Hoping it will help me find edge cases as I start to implement the miter limit an clipping step of the algorithm. To be honest though after playing with it for a bit, I am a little worried that some edge cases are going to be pretty challenging.

[michaelkirk]

As an aside, Wow, that's a really nice demo!

We've talked in the past about having interactive tools like this as a kind of documentation for some of our algorithms, but I think think this is the first actual working example. 👏 👏👏

[michaelkirk]

Hi @thehappycheese - this is an exciting new feature that I think a lot of people would like to use. Are you still interested in working on it? Is there anything I can do to help move it along?

[thehappycheese]

Hi @thehappycheese - this is an exciting new feature that I think a lot of people would like to use. Are you still interested in working on it? Is there anything I can do to help move it along?

Hi @michaelkirk 😃

I have to admit I stalled out on this project; There are a few contributing factors

1. While the basic offset algorithm is done, the clipping of the output part is not started.
   • Clipping the output geometry is pretty tricky, I'm feeling a bit apprehensive about tackling it
   • And also, my main motivating use-case for this feature does not actually require this additional clipping step
   • Perhaps this could just be an offset_curve_raw and in the future there could be an offset_curve_clipped?
2. I have introduced several new data structures and traits which overlap/conflict with existing code elsewhere in the project... On the one hand it makes sense to push forward and get the offset_curve function working then resolve these later. On the other hand, perhaps I should be trying to build up the required infrastructure with smaller separate pull requests first.
   • eg struct LineMeasured and the not-yet-implemented LineStringMeasured which are needed to avoid too many length calculations
   • eg fn line_segment_intersection_with_parameters
   • eg trait VectorExtensions
     • This one is particularly problematic... I think I should be trying to get these functions into Kernel or something
     • But I rather like the ergonomics of impl VectorExtensions for Coord but it is a really different approach
3. There are now two implementations of trait OffsetCurveOld (old) and trait OffsetCurve (new)
   • The new one uses trait LineStringOffsetSegmentPairs which was intended to make things tidy inside the trait OffsetCurve which I think it does... but it came at the cost of much headaches for me working out how to make iterators work.
   • Anyway I succeeded and now I think I will delete trait OffsetCurveOld
   • BUT the iterator trait LineStringOffsetSegmentPairs needs one more feature;
     • support wrapping over to the first segment at the end to work on closed shapes
     • I stopped before implementing that and took a break over Christmas and have not made it back here yet in the new year :/

I am still keen to get back into this but it might take me a few more weeks to get going again. In the meantime, if anyone is keen to jump in I would be happy to collaborate in any way possible :)

[thehappycheese]

My current draft PR is stalled because it got a bit too messy with too many new features at once.
Partly due to actually needing new features, and partly because I am unsure how or hesitant on using existing features.

In this comment I am trying to figure out if I can break off smaller bite-sized contributions and have a go at those first.
I'll probably add to this comment later

Required for Simple Offset

• In my current code I used a new trait impl VectorExtensions for Coord for the vector operations described in the dot points below. But VectorExtensions doesn't gel with the rest of the crate. To move forward I need to find replacements for the following functions

  • Desireable: 2D Wedge Product (aka exterior product, aka perpendicular product, aka 2d cross product) ( See )

    • Simplifies line intersection calculations see ( See example )
    • Currently implemented directly on geotypes::Point but its is not right it takes three arguments which is not what we need ( See )
    • It is related to Kernel::orient2d (which returns an enum) but we need the numeric value
    • Currently implemented as geotypes::Line::determinant() which is a bit wierd to use (Line::new(a, b).determinant() vs a.wedge(b))
    • nalgebra crate has perp for Vector2 like a.perp(&b)

  • Desireable: Dot Product

    • Currently implemented directly on geotypes::Point ( See )
    • Feels like it should live together with some of these other ops... should they all be in the impl Point or should it be in a trait?
    • nalgebra crate has dot for Vector2 like a.dot(&b)

  • Desireable: Vector Magnitude and MagnitudeSquared

    • Available as Line::line_euclidean_length, but this requires unnecessary construction of Line object
    • On further inspection separating these functions may not be needed... looks like with a bit of refactoring the EuclideanDistance trait could be used
    • I cant find a MagnitudeSquared... It is used for (perhaps unnecessary) micro-optimization to avoid calculating the square root when it is not needed.

  • Desireable: left and right functions for explicit 90 degree rotation around origin (for Coord or Point)

    • My current implementation uses Coord::<VectorExtensions>::left() and Coord::<VectorExtensions>::right() ( See )
      • These are "fast" and "precise" implementations because they just swap the x and y ordinates and negate one or the other
    • Alternatively there could be associated functions for Rotate trait ::left() and ::right()?
      • Actually I don't see a '::rotate_about_origin()'? is there a reason we only have ::rotate_around_point??
      • Rotate::rotate_about_origin(90.0) seems like it might be a reasonable implementation
      • Rotate::left() doesnt make sense for anything other than a point
    • Alternatively it could be a special constructor for struct AffineTransform ::left() and ::right().
      • Only really makes sense for Point/Coord
    • I still think a trait that explicitly allows 'left' and 'right' function directly on Coord and/or Point would help with clarity and testing

• Required: Specialized Line segment intersection;

  • In my current code I am using a special implementation that can
    • report if the intersection Real | Virtual( negative_side | positive_side )
    • the parameters t1 and t0 of the intersection points
  • Alternatively could augment existing algorithm/line_intersection.rs or add alternative implementation to that module
  • There is a private cartesian_intersect predicate in simplify_vw.rs ( See ) which may be redundant

• Required: trait LinesIter (already got this one)

  • Need a pairwise iterator to iterate over pairs of lines.
    • The standard library offers tuple_windows which clones elements
      • Line is Copy so this should be fine, I remember having trouble using it though
  • Maybe need a wrap_one implementation which would iterate over all segments plus the first one again
    • Not sure yet, but I have a feeling with wrap_one it would be trivial to implement offset_curve for closed shapes using almost the same code as open shapes

Required for Clipped Offset

• Required: split impl SplitBy<LineString> for LineString {fn split_by(self, other)->MultiLineString}
  • These discussions might help: LineString intersection API #985 and LineString slicing #986
  • In old python implementation I did this in two steps; first an algorithm that finds the 'parameters' where the intersection occurs ( See ) then an split-at-parameters function ( see ) note sure if that's the way to go but it seemed pretty efficient
• Required: 'CookieCut' operation; separate LineString by intersection with a circle
  • Required: a MultiCookieCut operation which does the same thing, but efficiently handles more than one circle at a time
    In a previous python implementation I used a Interval/MultiInterval subtraction algorithm to implement a function called remove_circles_from_linestring It worked fine but I think it was worse than O(n^2).
• Required: Efficient Method to stitch together a collection of LineStrings where the endpoints are touching
  • Looks like something similar LineAssembly in bool_ops/assembly.rs
• Required: something called the 'General closest point pair' by the paper referenced in my original post See Definition 6 and Figure 12 I got it implemented in python here
• Desireable: MeasuredLineString (and/or MeasuredLine?)
  • Reduce repeated magnitude calculations by storing segment length and/or overall line string length
  • Can mitigate the need for this by implementing supporting algorithms using "Segment-Fraction Line String Parameterization" where the integer part and fraction part of the parameter are used to refer to the segment number and segment fraction respectively.
  • Discussion here may be relevant Prepared Geometries #803

Point vs Coord, Z / M

• This old closed issue explains why we have both Coord and Point Coordinate vs. Point #15
• Since then discussion has been on implementing Geometries as Traits ( eg Introduce trait for accessing geometries #67 and introduce Coord/Point trait #901 and Trait-based Geometries #838)
  • The motivation being that a Point trait might be implemented on arbitrary structs that include Z and/or M ordinates in addition to X & Y.
• I found some recent discussion about geo_traits PointTrait with the suggestion being that it might harmonize Point / Coord?

With all that in mind, is the potential future trait Point the place to implement all the strictly 2D linear-algebra functions like dot, cross, perp, left/right etc?