[2D MVP] Milestone 3: Snapping in 2d viewer

Description

For setting the 2d vector for calibration we want to provide the user with snapping for convenience.

Goals

[ ] As a user I want to be able to snap to lines in the 2d viewer
[x] Snapping should occur as the cursor moves
[x] Cursor should snap to Corners, Intersections and Curves/Lines
[ ] There is a method that can be called externally that will return the snapped point
[ ] It is possible to enable or disable snapping

Tasks

[x] Add skeleton infrastructure for testing snapping behaviour
[x] Test polyfil to process paths into line segments
[x] Make functions to translate between container rect and svg coordinate system for svgImage
[x] Make cursor for debugging
[x] Find out which SVG elements need to be supported
[x] Add support for polyline
[x] Add support for path
[x] Choose a suitable acceleration structure for Line intersection tests
[ ] Get some more SVGs from ODA for testing
[x] Extract node-points from primitives for closest point tests and add to collector
[x] Add logic to compute intersections and add to collector
~Add acceleration structure for points~
[x] Integrate point tests into snapping logic and add type to results structure
[x] Add check to make sure svg element has width and height properly set (otherwise the mapping breaks)
[x] Make cursor update on zoom
~PolyCollector should not add duplicate points~
[x] Add Bezier Curve closest point test
[x] Add Curve & Line intersection test
[ ] Add Curve & Curve intersection test
[ ] Add curve self-intersection test
[ ] Add rect support
[ ] Add circle support
[ ] Add ellipse support
[ ] Add line support
[x] Add polyline support
[ ] Add polygon support
[x] Add path support

A few thoughts on this...

It seems there are a few approaches:

Image Space

This is what is used in the broadphase in the 3D viewer and involves applying an image feature detector to a patch from the Canvas.

This is fast and we know it works from the 3D viewer.

The main issue is that it cannot select the centre of lines or intersections (it can only select edges of the stroke, and corners of stroke interactions), and the precision is limited to whatever the rasterization resolution is, which depends on the zoom level.

In the 3D viewer anything snapped via this method is then refined using geometry tests.

DOM with getPointAtLength

This method will return points in the viewport for any SVGGeometryElement, which covers almost everything except reference items (embedded images).

This is the only DOM method currently support by all browsers.

The issue with this approach is that we get very little information about the underlying primitive. The only option is to sample all paths to an arbitrary precision, and perform tests on these segments.

This could end up very expensive, especially if we ever need to rebuild the set (e.g. to make higher precision one on zooming).

DOM with getPathData polyfill

The SVG Working Group's proposed getPathData with normalisation, will provide a minimal set of 2D primitives that it would be practical to write a series of tests against.

This is about as efficient as it is possible to get. The segments can be arranged into a bounding hierarchy for quick closest-point tests and intersection tests. The only primitive that would be computationally expensive to deal with would be the cubic beziers.

The problem is that this API is not supported by Chromium, and indeed Chromium decommissioned the alternative API to this almost 10 years ago, and hasn't seen fit to introduce the new one!

There is a polyfill available: https://github.com/jarek-foksa/path-data-polyfill, which works by effectively parsing bits of the SVG itself.

Relying on a polyfill that does its own parsing isn't great. Though one compensatory factor is that most of the SVGs we will handle will be produced by bouncer, so if we can make sure it works with whatever's output by ODA we should cover most use cases.

WASM based parser

The nuclear option is to get a reliable SVG parser in another language (e.g. the one from the Chromium source), build our snapping services around it, and cross-compile it into WASM.

This allows us to avoid the spotty support from browser manufacturers, and will give as good performance as we can get other than relying on user-agent implementations.

Other things to consider are that DOMParser isn't available in Web Workers, so for 2 & 3 we need to parse the SVG in the main thread (though I believe it's then possible to hand the graph to a web worker, but am not sure).

A brief update on this...

I have implemented the Harris Detector, which we chose for the 3D viewer.

It works well as a broadphase for SVGs. For example, here is the feature map configured to show Edges:

It also works with PNGs, at roughly the native resolution:

It does not work with blurry rasterisations of PNGs however. This is expected because the bilinear interpolation is designed to smooth out edges. We need to look at changing the sampling of PNGs, so that when zoomed in we perhaps use point-based sampling.

Importantly, it also does not work with single pixel lines...

This is expected because the Harris Detector is designed to detect the edge of regions.

The mechanics behind it are that the first step is to derive the image gradients, which we do with a Sobel operator. However, when running over a line 1 pixel wide, each side of the kernel of this operator will cancel the other out.

(Zooming in on a line we see it detected as a region with two edges...)

We can either try to get an image gradient that is sensitive to lines, or run a second detector dedicated to lines.

Additionally, the current implementation does not consider the opacity. This is important to consider to make sure we don't miss features because black lines are considered to have the same colour as "nothing".

Similarly to the 3D viewer, probably the best approach here is to run the detector again on the Opacity Map, so we specifically detect the edges of the primitives, and then that map is combined with the grayscale feature map, and Lines map.

Though, it may also be possible to do something clever with the opacity and negative numbers to make the Sobel end up outputting something sensible...

@sebjf if you want to compare expected behaviour, I think a good reference point may be trying the corner snapping in revit (on the 2d floor plans of course). It should behave fairly similarly to those situations.

If it makes your life any easier. I don't think we need to do snapping on raster images, as in we can probably just disable the feature. (@AsiteZane please clarify if I'm talking rubbish!) whole point of snapping is to get the accurate point on the drawing outside of pixel accuracy, and if the drawing is a png/jpg I'm not sure if there's much meaning behind doing it.

Agreed, we shouldn't do it for raster images.

Brief update/notes for self on this!

The polyfill appears to work on at least some images. Above, the lines (only lines, not curves) extracted from the svg are rendered to a separate canvas. This canvas is scaled (not the primitives), so is in the svg-space. We can see the lines are being extracted correctly, and we can map cursor positions between the viewer and the svg's local space.

We next need to determine a suitable acceleration structure.

Snapping Goals

Based on Revit, the features we'd like to identify are:

Endpoints of Lines and Curves
Intersection points of Lines and Curves
Closest point to lines

We are not interested in Text, or arrowheads.

Acceleration Structures

The user can zoom to arbitrary degrees, so cursor sensitivity defined in pixels will change in svg space dramatically (unlike for the 3d viewer where its all done in screen space).

BVH

A bounding volume hierarchy is a known quantity as we use it in number of places. BVHs can have highly efficient build times that do not depend on insertion order, depending on the construction method (one of the only structures for which this is true).

The advantages of BVHs (assuming AABBs - in this case we'd implement an R-Tree) are that the intersection tests are going to be very, very efficient in 2D. It can also be quite shallow if the number of nodes per level is relatively high (e.g. 4-8), which will help with memory overheads.

The main advantage is that it can handle curves more easily than other methods such as the Quadtree, as we only need to determine the maximum bounds for a curve.

BVHs can be searched at arbitrary distances using a Minkoswki sum depending on zoom level easily. They can also be self-intersected to find intersection points - although we would want to do this only once and store the intersections in another hierarchy.

Though the intersection tests themselves are going to be highly efficient, the tree is not as an efficient partition as something like a quadtree, due to the large potential overlaps between primitives leading a more expensive search for each query.

Quadtree

In Quadtrees, space is regularly subdivided as needed to reach a condition that minimises the number of primitives at each node. (This could be a count of primitives per node, or a maximum depth.) When lines are considered, they are split into q-edges, which are the segments of the lines that sit within the node.

Searching the structure is fast, as we can search a region of a particular size simply by picking the correct level in the tree to start from.

What is particularly interesting about this structure, is that if two lines intersect anywhere, they will always do so at a leaf node that has them (their q-edges) in common. So, after building the tree we can just take the leaf nodes as a set of optimised pairs to look for intersections with.

The main complication with quadtrees are that we would have to figure out how to make them work with Cubic Beziers (CBs). The same principles that apply to lines (i.e. q-lines, and the termination criteria) should be the same, but we'd need to work out what a q-edge of a curve looks like in memory, for example.

The other issue with quadtrees is the possibly considerable memory consumption in JavaScript, though this is likely to be true of any structure that can't be implemented as an array (BVHs are the only structure that could be implemented as an array).

We also need to consider the building time of the tree - though such trees could also be precomputed and cached, either by the client in indexeddb, or perhaps even on .io.

Additional notes on types of quadtree: https://people.cs.vt.edu/~shaffer/Papers/SametCVPR85.pdf

(For completeness, quadtree types are distinguished by their termination criteria..)

kd-tree

It is worth mentioning kd-trees, though they won't be a solution for everything because they don't support closest distance tests for lines/curves.

kd-trees split different dimensions at different levels. For 2D, we'd not use actual kd-trees, but either use a BSP (binary space partitioning), or create effectively something like an irregular or point quadtree, where a split occurs on both dimensions at each point.

Searching 2d-trees is very fast, but the insertion order determines the depth, which is not ideal.

Spatial Hash

A spatial hash would only be useful for looking at closest points, however in this regard it is by far the fastest to search.

The problem with such a hash is that it has a fixed cell size, meaning when zooming out more cells would need to be checked, counteracting some of its advantages.

A spatial hash would need all the intersection points pre-defined, which means we'd need to use another structure to find the intersections.

Continuing notes...

A proof of concept snapping implementation has been created using QuadTrees, that snaps a cursor in the image container to a line from the SVG.

There are many caveats to this:

It doesn't snap to points explicitly
It doesn't detect intersections
Its not optimised
The SVGs output by ODA don't work because the PolyLines are not picked up,
There are some slight offsets in the test code for the cursor
etc

However, the poc does fundamentally work, as can be seen...

The next step is to clean up the test visualisations, so we can be sure the results are accurate when converting between the coordinate systems (i.e. we place a point, then zoom in and see it maintains accuracy greater than sampling resolution when it was placed).

After this, the most important thing is to profile the code to check its performance, as if we need to do something like build a WASM library in C, or rely on async snapping via WebWorkers, we need to know that now.

If we are confident that this implementation can scale, then we can add in end-point snapping, and begin to detect intersections.

From what we can see of the performance so far, with the map with 80,000 line primitives the quadtree is built and sampled very quickly; though sampling it should be noted is only occurring on clicks at the moment, not each mouse move. The tree depth limit is important. With a depth of 24 for example the tree building process becomes unusable. With a tree depth limit of 14 (what its set to currently), it builds almost instantly and results in just over half a million nodes with the largest leaf node having 25 lines.

When testing with the carpark stress test, as 12 level quadtree takes about 5 seconds to produce. The intersection tests range from 1 ms to 30 ms, depending on the range.

This is because the quadtree with a max limit of 12 can end up with over 196 elements in the leaf nodes. A depth of 12 creates a little over 5 million nodes.

A binary r-tree by comparison would have about 1.6 million nodes across 19 levels.

Given the likelihood of encountering pathological svgs that mean the quadtree will top out the maximum depth in many cells, it is worth creating an r-tree implementation for comparison as well, as given the high leaf counts it may not take that much longer to traverse.

With an STR R-Tree implementation given 10 children per node, the tree takes about 1.5 seconds to build. The intersection tests range from 0 to 5 ms, depending on the range.

While in theory the quadtree should over better spatialisation, the STR - Sort-Tile-Recursive algorithm - offers a close approximation in 2D, without the extraneous nodes.

Behaving better in both build and query time then the R-Tree appears to be the best structure, certainly for the line-distance snap.

We will likely need further tests to determine if it is best for the points. Most likely it will still have best performance, but only when we have a separate tree constructed to make a best fit to the points. It is worth comparing it with a kd-tree however - the SRT building heuristic is not far off a kd-tree construction.

Another update...

A kd-tree has been implemented for closest points, which can be queried for the Revit Sample SVG and most others in between 0 and 0.3 ms. For the carpark, it is typically the same, though some locations can take up to 5 ms.

This performance is as expected. There is an optimisation available, which is to store the exact point that the split is made on, at each node, instead of having only leaf & branch nodes. This may allow earlier termination of the traversal. However the current profiles do not suggest this is necessary.

The main issue is with the construction time of the tree, which is about 17 seconds for the car park. (For the Revit Sample it is about 141 ms.)

This may be too long in absolute terms. In the first instance we can look at possible optimisations in layout for the structure. Other options include delivering a pre-built object, or offloading to a web worker (though it would still be 17 seconds before snapping would "work").

The next steps are to extract intersection points from the RTree, and build some logic to swap the cursors etc based on the intersection point type.

A helper class was constructed for building a predefined list of intersection points. However this immediately revealed some performance problems: the Revit Sample had over 2.9 million potential intersections, though this list took only a few ms to build.

The Carpark however failed to iterate all potential intersections at all. Simply counting the intersections took too long, and adding them to a list caused an overflow.

A better approach could be to drop the kd-tree entirely and perform both closest-edge, closest-node and closest-intersections using the R-Tree.

This would be done by tweaking the current traversal to return all leaf nodes that overlap the query circle. The closest-nodes would be extracted from the primitives in the same way as they'd be when building the point-acceleration structure. A second stage could then look for intersections as well.

Counting simply the number of overlapping nodes, the Revit Sample has up to a couple hundred nodes at max extents. The carpark however has over 38,000 leaf nodes at maximum extents.

Depending on how many nodes the browser can effectively check and maintain real-time performance (its probably much higher than we expect - note the closest edge tests are still running in this case, and take just over 5 ms for 36,000 nodes). One approach also could simply be to not snap if the number of nodes is above a certain threshold, as in this case there will be so many snapping points the user wouldn't be able to control it effectively regardless.

Updating the RTree to handle all intersections, the node and edge queries are as fast as before, meaning the kdTree is superfluous.

For intersection tests, an NxN intersection check of overlapping leaf nodes runs below a few ms in most cases, but in the case of the car park at maximum extents, can result in noticable lag...

IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 1919.6000000238419, numNodes: 19178, numIntersectingPairs: 9391760}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 97.30000001192093, numNodes: 5428, numIntersectingPairs: 424314}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 110.09999996423721, numNodes: 5599, numIntersectingPairs: 430918}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 108.39999997615814, numNodes: 5663, numIntersectingPairs: 439360}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 102.19999998807907, numNodes: 5832, numIntersectingPairs: 481451}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 360.2000000476837, numNodes: 8418, numIntersectingPairs: 2404347}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 101.5, numNodes: 5444, numIntersectingPairs: 432086}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 93.5, numNodes: 5340, numIntersectingPairs: 419875}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 95.09999996423721, numNodes: 5447, numIntersectingPairs: 420509}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 101.39999997615814, numNodes: 5724, numIntersectingPairs: 432318}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 63.80000001192093, numNodes: 4830, numIntersectingPairs: 323671}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 213.0999999642372, numNodes: 7002, numIntersectingPairs: 1439782}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 875.9000000357628, numNodes: 12975, numIntersectingPairs: 6126426}
IntersectionTestResults {closestEdge: Vector2, closestNode: Vector2, queryTime: 3642, numNodes: 25217, numIntersectingPairs: 19336215}

Times above are in ms.

numIntersectingPairs demonstrates that a a pre-pass is necessary (as otherwise, the number of pairs undergoing geometry tests would be much higher), but that the model still has a prohibitive number of actual primitives whose AABBs overlap. (These figures do not include actual intersection tests between primitives.)

The figure above illustrates the zoom level and approximate search radius at that level. That region contains 22608 individual primitives/17 million potential intersections.

With regards to the elements we need to support, it seems ODA predominantly uses the Path and Polyline parameters.

Still, the full list of SVGElement subclasses is as follows,

SVGElement
SVGSVGElement 
SVGCircleElement 
SVGClipPathElement 
SVGDefsElement 
SVGDescElement 
SVGEllipseElement 
SVGFEBlendElement 
SVGFEColorMatrixElement 
SVGFEComponentTransferElement 
SVGFECompositeElement 
SVGFEConvolveMatrixElement 
SVGFEDiffuseLightingElement 
SVGFEDisplacementMapElement 
SVGFEDistantLightElement 
SVGFEDropShadowElement 
SVGFEFloodElement 
SVGFEFuncAElement 
SVGFEFuncBElement 
SVGFEFuncGElement 
SVGFEFuncRElement 
SVGFEGaussianBlurElement 
SVGFEImageElement 
SVGFEMergeElement 
SVGFEMergeNodeElement 
SVGFEMorphologyElement 
SVGFEOffsetElement 
SVGFEPointLightElement 
SVGFESpecularLightingElement 
SVGFESpotLightElement 
SVGFETileElement 
SVGFETurbulenceElement 
SVGFilterElement 
SVGForeignObjectElement 
SVGGElement 
SVGImageElement 
SVGLineElement 
SVGLinearGradientElement 
SVGMarkerElement 
SVGMaskElement 
SVGMetadataElement 
SVGPathElement 
SVGPatternElement 
SVGPolygonElement 
SVGPolylineElement 
SVGRadialGradientElement 
SVGRectElement 
SVGStopElement 
SVGSwitchElement 
SVGSymbolElement 
SVGTextElement 
SVGTextPathElement 
SVGTSpanElement 
SVGUseElement 
SVGViewElement

Out of these, the following are edge-like geometric types (i.e. they can stroke lines or curves, and have nodes that can be snapped to). This includes primitives like circles, but not area fills like gradients, for example.

This list corresponds to the Basic Shapes defined by the SVG specification (and the SVGPathElement)

SVGRectElement
SVGCircleElement 
SVGEllipseElement
SVGLineElement
SVGPolylineElement
SVGPolygonElement
SVGPathElement

All the above shapes can be represented by an SVGPathElement and the functions to do so are defined in the SVG spec.

It may be worth however performing intersection tests against the primitives themselves, as these might be more efficient than against b-curves.

Honorary mentions include:

SVGDefsElement
SVGMarkerElement
SVGSymbolElement

These define elements or templates that are referenced elsewhere. For example, the SVGMarkerElement is like a def specifically for arrowheads and similar.

Since ODA doesn't seem to use the basic shapes much, if at all, (apart from polyline), the most efficient approach is to degenerate them to lines/bcurves to begin with.

Rect is likely as efficiently handled as 4x lines, in fact may end up considerably more efficient when we consider that a full rect would cover a larger area in the RTree.

Polyline/Polygon and Lines are just lines with multiple points.

The one shape that might be worth handling directly is Circle, since intersection tests with this will be much much faster than the 4x b-curve equivalent. Ellipses are probably best handled as bcurves.

Every time we add support for a specific shape the complexity of the solution increases significantly, as it needs an intersection test with every other distinct primitive supported. So lines and b-curves (the minimum) will be three tests (line-line, line-bcurve, bcurve-bcurve), as well as the closest-point tests for each individually.

Even adding circle support would increase this to six tests (line-line, line-bcurve, bcurve-bcurve, line-circle, bcurve-circle, circle-circle), and so on and so forth.

The next step for this is to add support for Bezier Curves. SVG specifies that Bezier Curves are cubic. The three tests required are

Closest Point
Line-Cubic Intersection
Cubic-Cubic Intersection

Out of these, the most algorithmically complex is the Closest Point.

The most performant method expresses the distance as a quintic polynomial and finds the roots. As there is no analytic method for finding the roots of a quintic (more on that later), an iterative method is used. The intervals to searched are culled based on [Sturms theorem]()https://en.wikipedia.org/wiki/Sturm%27s_theorem), with bisection to find the actual root within potential intervals.

There are a number of numeric methods based on finding intervals on the curve itself (e.g. here) but they are all slower than operating with polynomial roots.

Perhaps unintuitively, Line-Cubic Intersection and Cubic-Cubic Intersection are simpler. This is because the equation of the intersection point can be derived by substitution, and we only need to solve for the curve or line "parameter". For the line case, this involves solving a cubic function.

For the curve-curve case, it is a quantic function.

Equations exist for solving polynomials of order 4 or below, but not above. However, naively writing them out will result in stability issues due to quantisation errors. The solvers must be carefully implemented. (For example, these described by David Eberly).

Something worth considering before choosing an implementation, is that the Closest Point does not require the same level of precision as the intersection tests, owing to the fact the cursors precision will be dependent on the current zoom level.

There is a question of whether if we were to combine hueristics designed for cubic Beziers with some simpler numerical approach (Chen's method is designed for arbitrary order Beziers), the result might actually be more performant, so it is worth finding a few more algorithms to read first.

This blog post by nictvg is also really good, as it derives the unweighted cubic bezier equations, which is helpful, but also covers a few unexpected implementation details. For example, Sturms algorithm expects very high precision, which we don't control in JavaScript, so we might need another method of segmenting the intervals.

Following on from the above, there is now an implementation of closest point on Bezier, which is the most complex test we will have to do for snapping.

Some notes are below for documentation purposes.

The approach taken is the root-finding method, which seems to be the agreed upon optimal method. The idea behind this is to express the distance between a point and a curve as a polynomial. The equation that approximates this distance is given originally by Zhou as $\ (p-q(u)) \cdot q'(u) $.

Geometrically, this is the dot product of the vector between the query point p and a point the curve at u, and the tangent of the curve. Since the dot product is 0 when the vectors are perpendicular, this evaluates to 0 when the shortest line is perpendicular to the curve: a necessary condition for $q(u)$ to be the closest point.

(This is all well introduced by Chen et al.)

Since this expanded equation consists of dot products, the resulting polynomial is univariate for u.

(nictvg has kindly done the first part of working out the actual coefficients from the control & query points. The equation can also be split into two such that we can compute half the coefficients and store them for future tests ($-q(u) \cdot q'(u)$), then when we update the query point we have only three dot products to recompute $p \cdot q'(u)$.)

Below is a plot of a Bezier with varying query points (red circles).

The 'top' line is the distance between the curve and the query point, evaluated numerically for various points along u. The 'bottom' line is $(p-q(u)) \cdot q'(u)$ (also evaluated numerically, at the same points).

We can see from the above how this equation approximates the derivate of the distance. And, as one would expect with the derivative, there is a root wherever there is a local minima.

It is important to note however that it is not the actual derivative,

The above compares $\sqrt(|p-q(u)|)$ with that given above. As can be seen the roots are in the same place but the shape differs. Importantly, the sign is inverted. Chen et al propose pruning intervals based on the signs of the interval bounds, as they can be used to ascertain the shape:

However their test is inverted for us since they assert that $f'(u) = q(u)$ where $f(u) = || p-q(u) ||^2$, which we were unable to demonstrate as shown above.

Root Finder

The polynomial is quintic, and so we need an iterative solver of some sort. Most solutions use Sturms theorem, however Rouillier and Zimmerman propose a compelling method based on Vincent's theorem.

As described in the above, pretty much all methods operate as follows:

Bisect the polynomial into intervals
Count the roots in each interval
If that section has only one root, find the actual root through a naïve search
Repeat until all intervals have been found

The "outer loop" of this is the problem of root isolation, which gives a set of intervals containing only one root each. Once we have these intervals, a very simple algorithm - such as bisection comparing the signs of each split - can be used to find the true root within the interval. The simple algorithm can be used safely because there is guaranteed to only be one root, so it's not possible for the algorithm to 'miss' one.

What is particularly interesting to us about Rouillier and Zimmerman's approach is that it works on the principle of transforming polynomials, both to get a polynomial to which Descartes Rule-of-Signs can be applied to count roots, and to get a polynomial representing a specific interval.

(In this case the intervals are defined by a pair of integers - k & c. These are sort of like rationals, and applying a transformation to the polynomial parameterised by k & c ($P_{kc}(P(x))$), gives a new polynomial that represents just that interval (see Krandick). At the end, k & c are resolved to actual interval bounds on the original polynomial for finding the true roots.

These transformations can be fairly complex. However, they are true 'transformations', in the sense that applying the function to a polynomial results in another polynomial of the same order.

Coupled with the fact we only ever have to deal with one order (5th), we can take advantage of this and use a symbolic maths package (such as Matlabs Symbolic Toolbox) to get each function as a set of basic arithmetic operations applied to the polynomials coefficients, per power.

For example, define a 5th order polynomial where all the terms are symbolic variables:

P(x) =
 a5*x^5 + a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0

We can then write out the equations needed by Rouillier and Zimmerman, for example a Taylor Shift on P:

>> T(x) = P(x+1)
T(x) =
a0 + a1*(x + 1) + a2*(x + 1)^2 + a3*(x + 1)^3 + a4*(x + 1)^4 + a5*(x + 1)^5

And then have the package rearrange to get the equations for each coefficient:

>> collect(T)
ans(x) =
a5*x^5 + (a4 + 5*a5)*x^4 + (a3 + 4*a4 + 10*a5)*x^3 + (a2 + 3*a3 + 6*a4 + 10*a5)*x^2 + (a1 + 2*a2 + 3*a3 + 4*a4 + 5*a5)*x + a0 + a1 + a2 + a3 + a4 + a5

We can then just copy the resulting equations into JavaScript and have those functions available in Js, without the need to use any math packages or other imports.

This only works when the parameters of the functions are known however - for example, above we define the Taylor Shift $T_1(P(x))$ with a value of 1, but we wouldn't be able to define it for $T_n$.

The approach of Rouillier and Zimmerman however is to search the intervals in such a way as the polynomial representing each interval can be transformed into another adjacent interval in the search, by a small subset of operations.

Specifically, $H_i(P(x)) = P(xi)$ and $HTranslate_i(P(x)) = H_i(T_1(P(x))$

As these operations are applied, they change the polynomial to represent an interval over different k & c.

For example, transforming $H_2(P)$ transforms $P_{kc}$ into $P_{kc'}$ such that $k = k + 2$ and $c = 2^2 c$

For HTranslate, with a Taylor Shift of 1, the transforms are $ k = k + i $ and $ c = (c + 1)2^k `$

With these, we can make certain steps through the interval, for example to go from an interval defined by k=1, c=1 to k=2,c=2, we'd use $H_1$. To go from k=3,c=1 to k=2,c=1, we'd use $HTranslate_{-1}$: because $c' = (1 + 1) 2^{-1} = 1, and \space k' = k - 1 = 2$.

All steps through k & c have to be able to fit into these two functions in such a way, however Rouillier and Zimmerman demonstrate that this is the case when the search is conducted via a depth first traversal of a binary tree.

Fortunately, the scalars of x on H and HTranslate can be parameterised without exploding the complexity of the operations, and HTranslate always has a Taylor Shift of one. The other pre-transform before root counting ($T_1R(Q)$ from definition 3 in the paper), can also be found using the symbolic package in the same way as above.

With these in place we can isolate all the roots, and then it is a case of performing a bisection algorithm (such as Algorithm 1 from Chen et al.) to find the true root(s). These represent the points at u on q closest to p, so thereafter we can just take the closest point to p (or the end nodes).

Results

We can see this now functioning on a number of different curves, including those with cusps & self intersections (note Num Roots now applies Chen's pruning).

The number of iterations is between 3 and 9 for the found intervals. The total query time, including any initial coefficient calculations that cannot be cached, is between 0 and 0.1 ms for all the above.

An update, line-curve intersections have been added based on substituting the parametric equation of a curve into the implicit equation of a line.

For curve/curve tests, it is more complex than anticipated. The assertion before that curve/curves are quartic only was only for quadratic Beziers, but we deal with cubics.

There are a number of what are known as clipping algorithms for finding intersections, however some papers claim that for low order curves the root finding method is still optimal in time.

The root finding method is similar to the line/curve test, but in this case one of the curves is converted into implicit form, a non-trivial operation, but one for which Sederberg and Goldman in the above say they show is possible.

We would then need to find the roots of a 9th order polynomial. Fortunately, we have a working root finding class. We would need to create versions of the baked methods for 9th order polynomials in addition to 5th order ones (and also not use Chen's filtering) but otherwise the algorithm is identical. We can use the same approach as before to make the baked methods from symbolic expressions.

This seems to be the best approach as is both suggested to be fastest, and uses a lot of existing functionality, whereas other searching algorithms are complex and would need to be written from scratch.

3drepo / 3drepo.io