Paper.js: Optimize propagateWinding(), and improve confidence.
@iconexperience I think we should devote an issue to the attempt to reduce the amount of winding "sampling" per curve. At the moment, it checks three evenly distributed locations along a chain, and then takes the average winding. This is weird and probably also a bit slow due to the call to getWinding(), but has worked well so far for us.
If I reduce the amount of samples to just one in the middle of the chain, I get only one failing test. But if I increase the samples to 4 or more, I get the exact same test failing. Even with 10 samples, it's the only test that fails. This has me thinking that using a value of 3 is actually hiding another issue in this one test, and we may be better off finding out what that is?
lehni
All 93 comments
I've created the propagate-winding topic branch for easier testing, see 265cc71d2fdb78e2540ac9ab4bbc00462c1119bd
lehni
on 11 Jun 2016
What I find really weird is that the average is taken from the three points. So if two points have a winding of 2 and one has a winding of 0, the overall winding is 1? I think a majority voting makes more sense.
Also, even if we use three points, this can be optimized. If the first two windings are the same, the loop can stop. This should be the majority of cases, so this already should improve performance of propagate winding by about 30%.
iconexperience
on 11 Jun 2016
There used to be code that sorted the windings and took the middle one. It turned out to be the average always in ever case ever encountered, so I think I switched to it as the code is easier. But it's definitely misleading. The scenario you describe probably never occurs. I have a hunch the encountered winding values are never away from each other by more than one.
lehni
on 13 Jun 2016
But I don't think this matters anymore. I think we can just sample once in the middle of the chain, and use that. And this one issue that I'm observing has almost certainly another cause.
lehni
on 13 Jun 2016
@iconexperience and I think this has everything to do with you improving the getWinding() code as much as you did!
lehni
on 13 Jun 2016
So just to confirm: The simplification of only doing one sampling was possible, and no new issues arose!
lehni
on 13 Jun 2016
BTW, there are more things that can/should be optimized with propagateWinding() and getWinding(). Here is the first one:
In getWinding(), if we cast a horizontal ray through a curve and the curve is outside the x-tolerance of the point (in code:if ((v[0] < xBefore && v[2] < xBefore && v[4] < xBefore && v[6] < xBefore) || (v[0] > xAfter && v[2] > xAfter && v[4] > xAfter && v[6] > xAfter))` two optimizations are possible:
- We do not have to calculate the
xvalue of the crossing, because we are only interested if it is left or right ofpxWe could just usex = v[0]. - We do not need a y-monotone curve in this case. Any curve where
pyis betweenv[1]andv[7]can be treated as a monotone curve, with the windingw = v[1] > v[7] ? 1 : v[1] < v[7] : -1 : 0. And we can ignore the curve if(py < v[0]) && py < v[7] || (py > v[0] && py > v[7]). Here is an illustration:
At py=py1 the curve can be ignored, at py=py2 the curve can be treated as a y-monotone curve.
As you can see we need only very few monotone curves to determine the winding, most curves can be used as is.
The reason why I mention this is that I want to get rid of the ugly handling of horizontal curves and replace this with a method to cast a horizontal ray for more vertical curves and a vertical ray for more horizontal curves. For vertical rays we need x-monotone curves, so splitting a curve into monotone curves on demand will probably be more efficient than precalculating x-monotone and y-monotone curves for all curves.
iconexperience
on 13 Jun 2016
Oh, did I mention that I already have a working implementaion of this? :smile:
/**
* Adds the winding contribution of a curve to the already found windings. The curve does not have
* to be a monotone curve.
*
* @param v the values of the curve
* @param px x coordinate of the point to be examined
* @param py y coordinate of the point to be examined
* @param windings an array of length 2. Index 0 is the windings to the left, index 1 to the right.
* @param coord The coordinate direction of the cast ray (0 = x, 1 = y)
*/
function addWindingContribution(v, vPrev, px, py, windings, isOnCurve, coord) {
var epsilon = 2e-7;
var pa = coord ? py : px; // point's abscissa
var po = coord ? px : py; // point's ordinate
var vo0 = v[1 - coord],
vo3 = v[7 - coord];
if ((vo0 > po && vo3 > po) ||
(vo0 < po && vo3 < po)) {
// if curve is outside the ordinates' range, no intersection with the ray is possible
return v;
} else if (vo0 === vo3) {
// if curve does not change in ordinate direction, windings will be added by adjacent curves
return vPrev;
}
var aBefore = pa - epsilon;
var aAfter = pa + epsilon;
var va0 = v[coord],
va1 = v[2 + coord],
va2 = v[4 + coord],
va3 = v[6 + coord];
var roots = [];
var a = (va0 < aBefore && va1 < aBefore && va2 < aBefore && va3 < aBefore) ||
(va0 > aAfter && va1 > aAfter && va2 > aAfter && va3 > aAfter)
? (va0 + va3) / 2
: po === vo0 ? va0
: po === vo3 ? va3
: Curve.solveCubic(v, coord ? 0 : 1, po, roots, 0, 1) === 1
? Curve.getPoint(v, roots[0])[coord ? 'y' : 'x']
: (va0 + va3) / 2;
var winding = vo0 < vo3 ? 1 : -1;
var prevWinding = coord ? vPrev[0] < vPrev[6] ? 1 : -1 : vPrev[1] < vPrev[7] ? 1 : -1 ;
var prevAEnd = vPrev[coord ? 7 : 6];
var prevAStart = vPrev[coord ? 1 : 0];
if (a != null) {
if (po !== vo0) {
// standard case, the ray crosses the curve, but not at the start point
if (a < aBefore) {
windings[0] += winding;
} else if (a > aAfter) {
windings[1] += winding;
} else {
isOnCurve[0] = true;
windings[0] += winding;
windings[1] += winding;
}
} else if (a >= aBefore && a <= aAfter) {
if (prevAEnd >= aBefore && prevAEnd <= aAfter) {
if (winding !== prevWinding) {
if (prevAStart < v[6]) { // ToDo: This should be done with comparing tangens
windings[1] += 2 * winding;
} else {
windings[0] += 2 * winding;
}
}
} else {
if (a > prevAEnd) {
windings[1] += winding;
} else {
windings[0] += winding;
}
}
} else if (prevAEnd >= aBefore && prevAEnd <= aAfter) {
if (winding !== prevWinding) {
if (a < aBefore) {
windings[0] += 2 * winding;
} else if (a > aAfter) {
windings[1] += 2 * winding;
}
}
} else if ((pa - prevAEnd) * (pa - a) < 0) {
if (a < aBefore) {
windings[0] += winding;
} else if (a > aAfter) {
windings[1] += winding;
}
} else if (winding !== prevWinding) {
if (a < aBefore) {
windings[0] += winding;
} else if (a > aAfter) {
windings[1] += winding;
}
}
}
return v;
}
A little optimization is certainly possible and the way the windings are increased need some explanation. I will try to make a chart explaining each of the cases.
iconexperience
on 13 Jun 2016
And here is the function to split a single curve into x- or y-monotone curves:
splitToMonoCurves: function(v, coord) {
var vMono = [];
// getLength is a rather expensive check, therefore we test two cheap preconditions first
if (v[0] === v[6] && v[1] === v[7] && Curve.getLength(v) === 0)
return vMono;
var o0 = v[1 - coord],
o1 = v[3 - coord],
o2 = v[5 - coord],
o3 = v[7 - coord];
if (o0 >= o1 === o1 >= o2 && o1 >= o2 === o2 >= o3 || Curve.isStraight(v)) {
vMono.push(v);
} else {
var a = 3 * (o1 - o2) - o0 + o3,
b = 2 * (o0 + o2) - 4 * o1,
c = o1 - o0,
tMin = 4e-7,
tMax = 1 - tMin,
roots = [],
n = Numerical.solveQuadratic(a, b, c, roots, tMin, tMax);
if (n === 0) {
vMono.push(v);
} else {
roots.sort();
var t = roots[0],
parts = Curve.subdivide(v, t);
vMono.push(parts[0]);
if (n > 1) {
t = (roots[1] - t) / (1 - t);
parts = Curve.subdivide(parts[1], t);
vMono.push(parts[0]);
}
vMono.push(parts[1]);
}
}
return vMono;
}
Haha, you're crazy! And it works?
lehni
on 13 Jun 2016
What's this windings array, and how can it be combined with my changes from today?
lehni
on 13 Jun 2016
Yes, like a charm, at least for my cases. It's a bit messy because the values passed to addWindingContribution () can either be of a monotone curve or of a standard curve, depending on where the curve is located relative to the point.
The windings array contains the left and right winding. I have a getWinding() implemetation that can replace the current getWinding() function. The only change necessary is that the curves instead of the y-monotone curves will be passed as an argument. I need to clean that function up a little before posting here.
iconexperience
on 13 Jun 2016
Here is another thing I noticed about path propagation that currently causes some of my cases to fail. If two cruves are very close to each other (< gemoetical epsilon), the winding is 0 for one side and 2 for the other. For unite operations we replace this winding with 1, which works ok. But it can happen that the two curves diverge and at the end point the distance is greater than geometrical epsilon. In these cases the end point is no intersection and the winding gets propagated to adjacent curves, which causes incorrect results.
Here is an illustration of the problem:
At p1 the distance is < epsilon, so the winding is 1. But at p2 the distance is greater epsilon, so the winding is propagated to the next curve, which actually should have a winding of 2. So maybe in cases where the winding is 2 to one side and zero to the other, we should not propagate the winding. Not propagating the winding will not work, because the winding cannot change along a chain. So maybe we should use some quality metrics when determining the winding. If the quality is poor, we will try another point.
iconexperience
on 13 Jun 2016
Yeah it make sense that this kind of situation would make troubles... I'm not sure how to avoid those situations though. Do we need to sample a bit away form the contour instead of on it? Since the paths are oriented, we could step away bit a bit more than epsilon?
Also, please note we're not replacing the winding anymore for unite operations. I store the original winding plus the information wether the point is on an outer contour or not. Did you see those changes? I wrote about it today.
lehni
on 13 Jun 2016
I think that I have found the definitive solution for determining the winding number of a single point.
For determining the winding we cast a ray through a point and determine the crossings of the ray with the path. If the direction of the path is downward at the crossing, the crossing has the winding contribution of -1, If the direction is upward, the winding contribution is +1. For each side of the point the absolute value of the sum of winding contributions is calculated. The greater of both sums is the winding number of the point. Here is an example:
In this example the ray crosses the path three times on the right side of the point and once on the left side. The sum of winding contributions to the left is +1, the sum to the right is (-1) + (+1) + (-1) = -1. The overall winding number of the point is then calculated as max(abs(-1), abs(+1)) = 1.
Things get a little bit more interesting if the point is _on_ the path. A point p is defined to be on the path, if the distance between the crossing point of the path and the ray and the point p is less than epsilon, which in our case is defined as 2e-7. The path itself is considered to part of the the _inside_ of the path. For a point on the path the winding contribution is added to to the left _and_ right sum of winding contributions. In the following example the left winding is +1 and the right winding is (+1) + (-1) = 0. The overall winding number is max(abs(+1), abs(0)) = 1.
If the path is self touching at the point, the windings can cancel each other. Here are two examples of self touching paths. In the first example the winding number of the point is zero, because the winding contributions of the two crossings cancel each other. In the second example the overall winding number is 1 altough the windings at the point cancel each other, because there is an additional crossing left and right of the point.
Since the path itself is considered part of the inside, the winding number in the first example is not correct (0 mean _outside_ of the path for both possible fill rules). Therefore, if a point is on the path and the overall winding number is zero, the winding number has to be adjusted. This is done by using the windings that a point fully inside the path would have. For a path with clockwise orientation the windings are then be set to windingLeft = -1, windingRight = +1. For counter clockwise orientation of the path the winding are set to windingLeft = +1, windingRight = -1 (@lehni this is an important modification to the current algorithm).
iconexperience
on 18 Jun 2016
If there is more than one path, the windings for each side are calculated individually for each path and then added. Here is a standard case:
The winding sums are windingLeft = +1 and windingRight = -1 for the outer path and windingLeft = -1 and windingRight = +1 for the inner path. Therefore the overall winding sums are windingLeft = 0 and windingRight = 0. The overall winding number results in winding = max(abs(0), abs(0)) = 0, which means the point is outside the combination of both paths.
For comparison here is the same case, but with an opposite orientation for the inner path.
Here the windings are windingLeft=-2 and windingRight=+2. The overall winding number is winding = max(abs(-2), abs(+2)) = 2. This means the point is inside the combination of both paths.
If paths overlap at p, the left and wind windings are calculated for each path individually and then added. If the winding number for the combined path is 0, it is not adjusted. In the following example the orange part is an overlap of the paths. The left and right winding of a point on that part is 0, and therefore the overall winding is 0, too. This is the correct winding number because that part of the paths is not inside the combined path.
So the algorithm for finding the winding number of a point looks like this:
windingLeft = 0
windingRight = 0
for each path
calculate windingLeftPath
calculate windingRightPath
if windingLeftPath == 0 and windingRightPath == 0 and point is on path {
if path is clockwise {
windingLeftPath = -1
windingRightPath = +1
} else {
windingLeftPath = +1
windingRightPath = -1
}
}
windingLeft += windingLeftPath
windingRight += windingRightPath
}
winding = max(abs(windingLeft), abs(windingRight))
Technically the crossings of the ray with the path is done by first converting all curves of a path to a set of curves that are monotone in y-direction.
A y-monotone curve has a winding of +1 if its start point has a higher y value than the end point and a winding of -1 if the start point has a lower y value than the end point. A special case are curves where start and end points have the same y value (these curves must be lines), which have a winding of 0.
For calculating the winding number of a point each y-monotone curve is checked if it crosses a horizontal ray through the point's y-value. If there is a crossing, the x-value of the crossing is determined. If the x-value of the crossing is on the left side of the point, the winding of the y-montone curve is added to windingLeft. If it is on the right side of the point, it is added to windingRight. If the x-value of the crossing is the same as the x-value of the point (with a tolerance of epsilon), the winding of the curve is added to both windingLeft _and_ windingRight.
There are special cases when the horizontal ray crosses exactly at the start point of a y-montone curve (which usually is the end point of the previous y-monotone curve). In this case the ray crosses two y-monotone curves, one curve at the end point and one at the start point. This double counting of crossings would result in incorrect winding numbers, because two crossings would be counted althrough the path is only crossed once.. To prevent double counting of winding we usually ignore crossings of the ray with a curve if the crossing is at the start point of the curve. This is not always sufficient, so sometimes an adjustment is necessary if the crossing is at the end point of the curve.
The following table shows all possible cases of a ray crossing the path at the end/start point of two curves. In order to achieve the correct winding contribution of the crossing, in some cases an adjustment is necessary.
| id | type | correct winding | winding at curve end | necessary adjustment |
| :-: | :-: | :-: | :-: | :-: |
| 1 | 
| +1/+1 | +1/+1 | 0/0 |
| 2 | 
| 0/+1 | 0/+1 | 0/0 |
| 3 | 
| +1/+1 | 0/+1 | +1/0 |
| 4 | 
| +1/+1 | 0/+1 | +1/0 |
| 5 | 
| +1/+1 | +1/+1 | 0/0 |
| 6 | 
| 0/+1 | 0/+1 | 0/0 |
| 7 | 
| 0/0 | +1/+1 | -1/-1 |
| 8 | 
| 0/0 | 0/+1 | 0/-1 |
| 9 | 
| 0/0 | 0/+1 | 0/-1 |
| 10 | 
| 0/0 | 0/+1 | 0/-1 |
| 11 | 
| 0/0 | +1/+1 | -1/-1 |
| 12 | 
| 0/0 | 0/+1 | 0/-1 |
It may seem odd that for types 7 to 12 the left and right winding are zero. Actually,the correct values depend on whether the part of the path is convex or concave. As you can see in the following illustration, the winding numbers for this type should be 0/0 for concave situations and +1/-1 for convex situations.
Since we do not know if the path is convex or concave when evaluating the y-monotone curves, a winding of 0/0 will be used for both cases. If this results in an overall winding of 0/0for the whole path, the winding will be adjusted as described previously.
From the table above the following simple rules can be used to determine the the winding contribution of a crossing for a point p:
if crossing not at start point of curve {
if crossing left of p {
windingLeft++
} else if crossing right of p {
windingRight++
} else {
windingLeft++
windingRight++
}
} else { // crossing is at start point of curve
if winding has changed from previous curve {
undo winding contribution from previous curve
} else if there was a horizontal curve after last non-horizontal curve {
if crossing at p or p between this crossing and previous crossing {
add previous winding contribution, but on opposite side
}
}
}
Amazing work, @iconexperience, and beautifully illustrated! I need some time to digest this all, but fully trust you with doing the right thing here! In the meantime I shall look at that other boolean issue. I think I have an idea how to solve it.
lehni
on 19 Jun 2016
@hkrish you should check this out, I think you would approve!
lehni
on 19 Jun 2016
@iconexperience I was thinking that the addWinding() method would make the most sense as an inner function of getWinding() since it is only used there and needs to pass a lot of data back and forth through the array "hacks" (windings, onCurveWinding).
This would allow us to simply use outer scope variable, and name them something like windingLeft, windingRight, isOnPath, and modify them from this inner function.
But I'll wait with modifying more of the code until you've implemented the suggestions above, because it sounds like a lot may change again.
lehni
on 19 Jun 2016
@lehni This is not completely finished. As you saw in my code I am casting a horizontal or vertical ray, depending on the tangent of the curve at p. I still need make a little illustration what _left_, _right_, _+1_ and _-1_ mean in case of a vertical ray.
Regarding the getWinding() / addWindingContribution() functions: I separated these so addWindingContribution() is the one to handle a monotone curve (or a curve that can be treated as monotone). I thought that 'getWinding()' may become a little too big, but we can see at the end if these two should better be combined into one.
iconexperience
on 19 Jun 2016
Nice write up and well illustrated! I need some time to read this through.
hkrish
on 19 Jun 2016
@iconexperience finally finding a bit of time for this again. Where do we stand currently in all this?
lehni
on 6 Jul 2016
I have an implementation up and running, just need to find the time to clean it up and create a PR.
iconexperience
on 6 Jul 2016
@iconexperience the implementation you mentioned here is the one I've already merged in now, yes?
lehni
on 16 Jul 2016
@iconexperience you mentioned your new implementation of getWinding() here: https://github.com/paperjs/paper.js/issues/1075#issuecomment-233200314 Could you share the code with me so I can test this and close this issue also?
Update: I meant propagateWinding()
lehni
on 21 Jul 2016
You mean the implementation of getInteriorPoint()? Here it is:
getInteriorPoint: function() {
var bounds = this.getBounds(),
point = bounds.getCenter(true);
if (!this.contains(point)) {
var curves = this.getCurves(),
y = point.y,
intercepts = [],
monoCurves = [],
winding;
for (var i = 0, l = curves.length; i < l; i++) {
var v = curves[i].getValues();
if ((v[1] <= y || v[3] <= y || v[5] <= y || v[7] <= y) &&
(v[1] >= y || v[3] >= y || v[5] >= y || v[7] >= y)) {
var monoParts = Curve.getMonoCurves(v, 0);
for (var j = 0; j < monoParts.length; j++) {
var monoV = monoParts[j];
if (y >= monoV[1] && y <= monoV[7]
|| y >= monoV[7] && y <= monoV[1]) {
winding = monoV[1] > monoV[7] ? 1 : monoV[1] < monoV[7] ? -1 : 0;
if (winding) {
monoCurves.push({values: monoV, winding: winding});
}
}
}
}
}
if (!monoCurves.length) {
return point;
}
for (var i = 0, l = monoCurves.length; i < l; i++) {
var v = monoCurves[i].values;
var roots = [];
var x = y === v[1] ? v[0]
: y === v[7] ? v[6]
: Curve.solveCubic(v, 1, y, roots, 0, 1) === 1
? Curve.getPoint(v, roots[0]).x
: (v[0] + v[6]) / 2;
intercepts.push(x);
}
intercepts.sort(function(a, b) {
return a - b;
});
point.x = (intercepts[0] + intercepts[1]) / 2;
}
return point;
}
My implementation of getWinding() should be the same as the one I posted before (the one that you reformatted).
iconexperience
on 21 Jul 2016
Sorry, I meant propagateWinding()... onPathCount & co, as mentioned here: https://github.com/paperjs/paper.js/issues/1075#issuecomment-233201975
But this issue is about propagateWinding(), so I got them confused again.
lehni
on 21 Jul 2016
Here is propagateWinding():
function propagateWinding(segment, path1, path2, monoCurves, operator) {
var chain = [],
curve,
curveLength,
start = segment,
totalLength = 0,
winding;
function roughLength(crv) {
var p1 = crv.point1,
h1 = crv.handle1,
h2 = crv.handle2,
p2 = crv.point2,
dx = p2.x + h2.x - p1.x - h1.x,
dy = p2.y + h2.y - p1.y - h1.y;
return h1.length + h2.length + Math.sqrt(dx * dx + dy * dy);
}
do {
curve = segment.getCurve();
if (curve) {
curveLength = roughLength(curve);
chain.push({segment: segment, curve: curve, length: curveLength});
totalLength += curveLength;
}
segment = segment.getNext();
} while (segment && !segment._intersection && segment !== start);
var chainLength = chain.length;
var length = totalLength * 0.47;
for (var i = 0; i < chainLength; i++) {
curveLength = chain[i].length;
if (length <= curveLength + Numerical.MACHINE_EPSILON) {
curve = chain[i].curve;
var path = curve._path,
parent = path._parent;
for (var j = 0; j < 3; j++) {
var t = 0.5;//j == 0 ? 0.95 : j == 2 ? 0.05 : 0.47;
var pt = curve.getPointAtTime(t),
hor = Math.abs(curve.getTangentAtTime(t).normalize().y) < 0.5;
if (parent instanceof CompoundPath)
path = parent;
var w = !(operator.subtract && path2 && (
path === path1 && path2._getWinding(pt, hor) ||
path === path2 && !path1._getWinding(pt, hor)))
? getWinding(pt, monoCurves, hor, operator)
: {winding: 0, confidence: true};
if (!winding || winding.confidence < w.confidence) {
winding = w;
}
if (winding.confidence === 1) break;
}
break;
}
length -= curveLength;
}
for (var j = chainLength - 1; j >= 0; j--) {
var seg = chain[j].segment;
seg._winding = winding.winding;
seg._contour = winding.contour;
}
}
You must change getWinding()so it returns this:
return {
winding: max(windingL, windingR),
contour: !windingL ^ !windingR,
confidence: onPathCount <= 1 ? 1 : 0
};
For performance improvement I implemented roughLength(), but as it turned out Chrome cannot optimize this so it wasn't faster at all.
iconexperience
on 21 Jul 2016
Thanks! is onPathCount the same as onPathWinding? Or is it a counter that's counted up each time you add to / remove from onPathWinding? I don't have that value here.
lehni
on 21 Jul 2016
Right, this is different in my getWinding(). onPathCount is incremented each time if the point is on a path. You will have to add this variable and insert if (isOnPath) onPathCount++; before this line:
https://github.com/paperjs/paper.js/blob/develop/src/path/PathItem.Boolean.js#L547
iconexperience
on 21 Jul 2016
@iconexperience thanks! I'll experiment with this a bit. Perfect!
lehni
on 21 Jul 2016
So you're picking the curve at length = totalLength * 0.47, and then you sample three points within that one curve, at t = 0.047, 0.5, 0.95. Any reason why you don't do this across the full chain instead, at length = totalLength * 0.47, * 0.5, * 0.95?
lehni
on 21 Jul 2016
Sorry, I don't remember if there is a specific reason for that. Maybe there was a tiny curve at 0.5 in one case which caused a problem, so I selected 0.47? But probably that 0.47 is not necessary anymore, because getWinding() now produces much more reliable results than before, when I introduced that.
iconexperience
on 21 Jul 2016
I'm not talking about the difference between 0.5 vs 0.47. I'm talking about the distribution: Why not spread it out across the full chain, instead of just the one curve in the middle of the chain, which is what your current code is doing.
This spreads it out evenly, at 0.25, 0.5, 0.75 across the full chain.
I've implemented it as a sub-function for now so I can toy around with different offset more easily than if it was in one loop, but eventually this could be optimized again:
function propagateWinding(segment, path1, path2, curves, operator) {
// Here we try to determine the most likely winding number contribution
// for the curve-chain starting with this segment. Once we have enough
// confidence in the winding contribution, we can propagate it until the
// next intersection or end of a curve chain.
var chain = [],
start = segment,
totalLength = 0;
function getWindingAt(length) {
for (var j = 0, l = chain.length; j < l; j++) {
var entry = chain[j],
curveLength = entry.length;
if (length <= curveLength) {
var curve = entry.curve,
path = curve._path,
parent = path._parent,
t = curve.getTimeAt(length),
pt = curve.getPointAtTime(t),
// Determine the direction in which to check the winding
// from the point (horizontal or vertical), based on the
// curve's direction at that point.
dir = abs(curve.getTangentAtTime(t).normalize().y) < 0.5
? 1 : 0;
if (parent instanceof CompoundPath)
path = parent;
// While subtracting, we need to omit this curve if it is
// contributing to the second operand and is outside the
// first operand.
return !(operator.subtract && path2 && (
path === path1 && path2._getWinding(pt, dir) ||
path === path2 && !path1._getWinding(pt, dir)))
? getWinding(pt, curves, dir)
: { winding: 0, confidence: true };
}
length -= curveLength;
}
}
do {
var curve = segment.getCurve(),
length = curve.getLength();
chain.push({ segment: segment, curve: curve, length: length });
totalLength += length;
segment = segment.getNext();
} while (segment && !segment._intersection && segment !== start);
// Sample the point at 0.25, 0.5 and 0.75 of the length of the chain
// to get its winding: #1075#issuecomment-233201975
var winding = getWindingAt(totalLength * 0.25);
if (!winding.confidence) {
winding = getWindingAt(totalLength * 0.5);
if (!winding.confidence) {
winding = getWindingAt(totalLength * 0.75);
}
}
// Now assign the winding to the entire curve chain.
for (var j = chain.length - 1; j >= 0; j--) {
var seg = chain[j].segment;
seg._winding = winding.winding;
seg._contour = winding.contour;
}
}
If you look at the example here, you will see that there is a higher chance of getting the correct winding if we use samples near the start and end of the curve. Many cases that caused trouble were exactly of that kind.
But I am not saying that I found the best solution, it just solved a lot of issues that I had.
iconexperience
on 21 Jul 2016
I see... But generally speaking, without optimizing for particular cases, isn't it better to distribute evenly across the full chain, and sample in more places perhaps, and keep looping until we find a location with good confidence? Most of the time, the first location will be fine, probably?
lehni
on 21 Jul 2016
You ma be right. I tried a lot of things and this gave the best results for me. If we have time to test and think more, I am sure we will find a better solution.
BTW, here is one of the cases that cause me trouble.
There is a crossing at A, so the green path must run at the outside. But if you zoom in at B, you can see that the green path is on the inside. But there is no crossing at C, which does not really fit. These are the cases where a "confidence" factor really makes sense. Of course it would also be nice to have a better fat line clipping that can find intersections on almost identical curves, but this may be difficult (but maybe possible?).
iconexperience
on 21 Jul 2016
I have a plan for improving tracePath() to handle exactly such cases much better. In short, I plan on keeping a tree structure in memory that branches each time tracePath() crosses at an intersection, where the current segment would also be a valid choice.
If it ends up in no man's land (visited segment encountered, etc), it then walks back to the last crossing, marks the visited segments as not visited and removes them from the resulting path gain, and then goes down the other path.
Going down the other path then does not mean creating a new branch, it means staying on the path of the previous branch, to which the new segments will be added.
If this too ends up in no man's land, it can just repeat the steps described above until no branches are left to choose from, and only then it'll give up.
I'm pretty confident that this should help us with handling many of the tricky edge cases that we're still encountering. It's also what I think will help fix the tricky scenario in #1091
lehni
on 21 Jul 2016
Yes, a tree structure would be really great. Meanwhile maybe I can look at the fat line clipping. I think I remember that there was an addition on how to handle tangential cases, which I did not undestand at that time. But I have a much better undestanding of the algorithm now.
iconexperience
on 21 Jul 2016
Very cool, that would be great!
lehni
on 21 Jul 2016
@iconexperience I've just implemented the branches mechanism described above: 445d8ae22f78f3cc1f98ee15f1c81dd485ddde1f.
The work's not finished yet, as I am still observing issues with some edge cases with it, but it does seem to have solved some other cases at my end, and not produced new ones.
I am a bit confused about the sketch you posted in https://github.com/paperjs/paper.js/issues/1073#issuecomment-234305530, since that already seems to work in v0.10.2, and still does with my new code. What am I missing?
lehni
on 27 Jul 2016
@lehni This is a misunderstanding, I did not want to say that the example fails, but I think it was one of the tests that failed when I implemented the new getWinding() code. Even if it does not fail, you have to admit that this case is a bit awkward with the path changing sides without any intersections found :)
iconexperience
on 27 Jul 2016
Is this a problem of isCrossing() then? C should be a crossing, but isn't seen as one? Not fully understanding your analysis of the case.
lehni
on 27 Jul 2016
Maybe, or a problem of Curve.getIntersections(). But either way, I did not see a way to improve isCrossing() or getIntersections(), so I needed to find another solution, and the confidence value helped. I did not really approach this problem from an academic point of view. What I needed was a quick solution to make some edge cases work, and adding the confidence value seemed to help.
iconexperience
on 28 Jul 2016
But the right intersections seem to be found, no? It's probably just that it doesn't think that it's crossing. Maybe because the lines are almost parallel....
lehni
on 28 Jul 2016
It'd be good to know how the new code I pushed today works for your edge cases.
lehni
on 28 Jul 2016
Maybe the correct intersection was found, but not at the correct location, so isCrossing() was not able to give the correct result? I do not know. At the moment, these almost identical curves are the problem that worries me the most. If your new tracing algorithm can resolve most of these edge cases, this would be just great.
I will test your new code tomorrrowtoday.
iconexperience
on 28 Jul 2016
@iconexperience the reliability of the found intersection is easy to verify, by checking the distance of the found points on both paths. I have a hunch it's the right location, but an imprecision in isCrossing()...
lehni
on 28 Jul 2016
@iconexperience I can confirm my hunch:
I changed the length multiplier in my version of isCrossing() from 1/64 to 1/32, and the intersection turns into a crossing.
I believe what https://bitbucket.org/andyfinnell/vectorboolean does is probably what we end up doing as well there: Move away from the intersection until all tangents are unambiguous...
I'll have to loko into this more, as part of wrapping up #1074
lehni
on 28 Jul 2016
@iconexperience your version of isCrossing() suffers from the same problem as mine...
lehni
on 28 Jul 2016
I finally found some time to extract a failing example in issue #1174.
The problem is that if you reduce the length multiplier further, you will not recognize cases as in the example in #1074. If you increase the length multiplier, you will falsely report non-crossing intersections as crossings, as you can see in #1174.
We have a dilemma here.
iconexperience
on 23 Sep 2016
@lehni Today I learned how the winding information has to be used correctly for our unite operations. As you know all segments that have a winding of 1 or are part of the contour are also part of the result in a unite operation.
As it turns out, the definition of the contour is actually quite simple. It is that part of a path where one side of the path is filled and the other side isn't. The contour depends on the fill rule. For even-odd a segment is on the contour if winding % 2 == 0 is true on one side and false on the other. For non-zero the segment is on the contour if winding == 0 is true on one side and false on the other.
Actually, this could probably be the only condition to decide if a segment should be part of the result in a unite operation, without checking the segments's winding number.
Currently we set onContour: !windingL ^ !windingR, which means we are doing it correctly for non-zero winding rule, but not for even-odd. But I don't think we have to change this now.
iconexperience
on 4 Jan 2017
@iconexperience very interesting! Why do you think we don't need to change it for even-odd?
We could change isValid() quite easily, and include your comments there for better explanation
lehni
on 4 Jan 2017
Just not now, because non-zero is our default fill rule. We should correctly implement it for even-odd eventually.
BTW, I am planning to post my new getWinding() and propagateWinding() code that I promised long time ago tomorrow.
iconexperience
on 4 Jan 2017
Sounds good! Regarding isValid(), I just tried the test suite with this change, and I am getting 54 failing tests:
function isValid(seg) {
var winding = seg._winding || {};
return !!(seg && !seg._visited && (!operator
|| operator.unite
? winding.onContour
: operator[winding.winding]));
}
Yes, there is something wrong, either with my assumption or with the implementation.
iconexperience
on 4 Jan 2017
Since it's not causing problems currently, we may as well stick with the current implementation and not further wonder about this :)
lehni
on 5 Jan 2017
While we're at it:
- I still have
onPathCountin the current code base, but it's not in use anywhere. Are you using this still to determine confidence? Where have we left that discussion? - I remember the discussion about the edge case of two nested circles touching... That's still undecided too, right?
lehni
on 5 Jan 2017
Regarding the onPathCount, this is not used anymore by me, only onPath. I am not sure where the dicussion was last time, but let me post the new code before we continue.
I have solved the case with the two nested circles for me simply by checking if a point is on the path. If so, contains() always returns true. You will see in the new code.
iconexperience
on 5 Jan 2017
@lehni I have pushed my changes to the new branch improvedWinding2 in my fork. The changes are not really that big, but believe me, this handles a lot of endge cases better than the current develop branch.
iconexperience
on 5 Jan 2017
Great! I will look at it shortly.
lehni
on 5 Jan 2017
@iconexperience is this pending here?
//ToDo:
quality = 0;
Also, what's the 100 * based on here?
if (a > pa - 100 * epsilon && a < pa + 100 * epsilon) {
1e-12 * 100 = 1e-10. Is that the epsilon that works best there?
lehni
on 5 Jan 2017
Some more: There's an unfinished comment here. Did it get deleted?
// If the winding one
And lastly, going back to these offsets:
var offsets = [0.48, 0.1, 0.9];
Did you compare with even ones, as discussed further up? e.g. at 1/4, 1/2 and 3/4?
lehni
on 5 Jan 2017
Yes, the whole calculation of the quality factor is unfinished, but it already works well. I wrote this specific line probably a few months ago, so I have to figure out again what it means.
iconexperience
on 5 Jan 2017
The broken comment was supposed to explain what we are doing with the windings. I started it shortly before the commit but forgot to finish.
iconexperience
on 5 Jan 2017
Here is a little more background on the quality factor: One weak point of our algorithm are cases where multiple curves are very close to each other. In these cases often a wrong winding is detected. Since I still do not know how to handle these cases correctly, I am trying to find locations on the wining-chains where the curves have a certain distance to each other. That's why each time our horizontal/vertical ray crosses a curve very close to the specified point, the quality factor is divided by two. A quality factor of 0.5 indicates that the point is on one curve, which means the winding will probably be determined correctly. A quality factor of less than 0.5 indicates that the ray crosses at least two curves very close to the specified point, which makes it difficult to determine the correct winding.
iconexperience
on 5 Jan 2017
Here is an example that results in a wrong winding number with the current algorithm. Both paths are clockwise, but note the flap at the bottom of the green path.
If we cast a vertical ray through a certain point this can happen:
The three crossings near the point are not counted, because their distance to the point is less than WINDING_EPSILON. Therefore we have a winding of 2 to the top and a winding of -1 to the bottom. This results in an overall winding of 2 and the condition for being on the contour is not fulfilled (at least not for non-zero fill rule). This makes the segment not valid in tracePaths with a unite operation, although it should be.
iconexperience
on 5 Jan 2017
@iconexperience thanks for clarifying! Should I merge this in already, and we carry on from there?
lehni
on 5 Jan 2017
Oh, and I am still curious about the origins of the epsilon * 100 (= 1e-10) described above!
lehni
on 5 Jan 2017
The epsilon * 100 was the threshold were most glitches disappeared. I tested this on a lot of pairs of almost identical paths, which are quite hard to get right. As I mentined above, it's not finished but works surprisingly well.
iconexperience
on 5 Jan 2017
I would be interested in a test case with these edge cases, also to build some unit tests out of them. Do you have them packaged?
lehni
on 5 Jan 2017
Unfortunately I never really extracted the test cases. All I did was to look at the results of a few thousand path operations and tweak the code until all glitches disappeared.
Here is one test case that is very sensitive to scaling:
var p1 = new Path({segments:[[890.91, 7.52, 0, 0, 85.40999999999997, 94.02], [1024, 351.78, 0, -127.03999999999996, 0, 127.14999999999998], [890.69, 696.28, 85.54999999999995, -94.03999999999996, 0, 0], [843.44, 653.28, 0, 0, 75.20000000000005, -82.66999999999996], [960, 351.78, 0, 111.75, 0, -111.63], [843.65, 50.51999999999998, 75.07000000000005, 82.63999999999999, 0, 0]], closed:true})
var p2 = new Path({segments:[[960, 352, -0.05999999999994543, 111.67999999999995, 0, 0], [1024, 352, 0, 0, -0.05999999999994543, 127.07], [890.69, 696.28, 85.5, -93.98000000000002, 0, 0], [843.44, 653.28, 0, 0, 75.14999999999998, -82.61000000000001]], closed:true})
var scale = 0.2647557802943046;
p1.scale(scale, [0, 0]);
p2.scale(scale, [0, 0]);
var res = p1.unite(p2);
p1.strokeColor = "green";
p2.strokeColor = "red";
res.fillColor = "rgba(255, 255, 0, 0.3)";
I do think that we need further tests for these edge cases... I also would like to look a bit more into these seemingly random values, e.g. the best epsilon for the quality test, and the three length factors at which we sample. I've brought this up before, but these seem really strange to me: 0.48, 0.1, 0.9 I am not really comfortable with code that avoids 0.5 because segments occurs more commonly there than at 0.48... So I would be curious to test different approaches there and be able to compare how it performs on these tests.
We also can't really assume that there is no bias in the path operations we run, e.g. resulting from the way the paths are crafted. It could for example be that their scale, or size of cusps, loops, etc, may favor certain epsilons over others, while in a more broad reality they wouldn't hold such an advantage.
lehni
on 6 Jan 2017
Having said all that, I will merge this in now since it already performs better, and then further work can happen on the develop branch. It would be good to finish these TODOs and comments though, to make sure the code is ready for the master branch.
lehni
on 6 Jan 2017
Please note also that your commit changes the formatting in many places, which makes it harder to see the actual change. It would be good if the code formatting could remain the same.
(Really appreciating the work of course, just pointing out something that would make merging easier)
For example, there is a reason that pre-compile directives aren't aligned with the rest of the code, because they are applied on build time, e.g.:
/*#*/ if (__options.nativeContains || !__options.booleanOperations) {
@iconexperience I've merged it in with reformatting now, see 7583e6ed5f1be63b1df6a03e6e3ae60f044cd6b0, but for some reason, it didn't keep the reference to your commit. I am still learning about Git as well. ARGH.
lehni
on 6 Jan 2017
Ok, I've "hacked" your commit in now. Historic consistency maintained. :)
lehni
on 6 Jan 2017
@iconexperience with these changes merged in, your test in https://github.com/paperjs/paper.js/issues/1073#issuecomment-270866689 still fails. Is this expected?
lehni
on 6 Jan 2017
And more food for thought... Out of curiosity, I tried these values for offsets:
var offsets = [0.5];
This is the same as we used before merging this in. With this, I now get 5 failing cases. This makes me wonder if some of the code you removed from getWinding() should have staid, as it may produce more reliable results in edge cases?
Out of curiosity I then also tried this, which is the same we used to do a while back:
var offsets = [0.5, 0.25, 0.75];
With this I still get 4 failing cases, and that's those where many values for t end up directly on 1. I think we could simply add some clamping there to avoid those values for t. But I am also curious as to why this now fails, and was working better before.
lehni
on 6 Jan 2017
...and some more: If I do clamp the values for t like this, then I get these 4 failing cases down to 1:
t = Numerical.clamp(curve.getTimeAt(length),
Numerical.CURVETIME_EPSILON,
1 - Numerical.CURVETIME_EPSILON),
I do feel that we need to get to the bottom of this more. I am only comfortable with a new solution when we can give it [0.5, 0.25, 0.75] as offsets and it performs well... Otherwise we're just masquerading a potential issue behind obscure offsets.
lehni
on 6 Jan 2017
I made an interesting discovery:
I compared the code before 3c2588fdec1035d3960b46d5c7271d6ed5e2b7ff with after, and started bringing things back, then starting removing them again, to see what helps with my edge cases. It turns out this simplified version of the prior code inserted at PathItem.Boolean.js#L599, together with the above clamping, is all that's required to fix those:
if (onPath && !windingL && !windingR) {
// If the point is on the path and the windings canceled
// each other, we treat the point as if it was inside the
// path. A point inside a path has a winding of [+1,-1]
// for clockwise and [-1,+1] for counter-clockwise paths.
// If the ray is cast in y direction (dir == 1), the
// windings always have opposite sign.
var add = path.isClockwise(closed) ^ dir ? 1 : -1;
windingL += add;
windingR += add;
}
@iconexperience how do you feel about bringing this back?
lehni
on 6 Jan 2017
@iconexperience I forgot to mention a second change above:
I removed this line: https://github.com/paperjs/paper.js/blob/3c2588fdec1035d3960b46d5c7271d6ed5e2b7ff/src/path/PathItem.Boolean.js#L461
And brought this back instead, after https://github.com/paperjs/paper.js/blob/develop/src/path/PathItem.Boolean.js#L490
} else {
windingL += winding;
windingR += winding;
onPath = true;
}
And this too fixed an edge case.
lehni
on 6 Jan 2017
@iconexperience I've commit the suggested changes. Could you test this against your code-base and let me know if it still performs well? We can undo any of it again, of course. This still fails for me though: https://github.com/paperjs/paper.js/issues/1073#issuecomment-270866689. Interestingly enough, that edge case works in the current release: Sketch
lehni
on 6 Jan 2017
@iconexperience I created the separate issue #1239 for the new edge case you posted above, and it's already fixed. I closed this issue now, but I am going to reopen as I am still waiting on feedback from you regarding the changes I made to the windings code, outline here onwards:
https://github.com/paperjs/paper.js/issues/1073#issuecomment-270888992
More importantly here onwards:
https://github.com/paperjs/paper.js/issues/1073#issuecomment-270897288
lehni
on 25 Jan 2017
@iconexperience still waiting for feedback, and I would love to close this...
lehni
on 4 Feb 2017
@iconexperience also please note that I consider the code handling quality incomplete, especially due to the lack of finished comments... I would love to improve things there and fix the pendnig // TODO comments.
lehni
on 4 Feb 2017
@iconexperience as outlined in https://github.com/paperjs/paper.js/issues/1261#issuecomment-282520231, I had to remove the code reintroduced in https://github.com/paperjs/paper.js/issues/1073#issuecomment-270897288 again to solve another issue. I am still curious to hear back if this code was wrongly reintroduced or not...
lehni
on 26 Feb 2017
@lehni I am a bit confused right now. All the open issues seem to work with the current version, so I am not really sure what to look at. Of course, we have to discuss the "quality" factor, but what else needs to be looked at?
iconexperience
on 26 Feb 2017
@lehni Here is at least an answer to one open questions. The factor 0.48 was introduced when there were aquite a few more bugs in the boolean code. At that time the windings at a corner (directly on a segment point) could not be calculated reliably. And in many regular shapes (squares for example) there is a corner exactly at 0.5. A corner is much more unlikely to appear at 0.48, so that's the reason why I chose that value. I am quite confident that the problems have been solved meanwhile, so I would be comfortable with t=0.5.
iconexperience
on 27 Feb 2017
@iconexperience I think you gave all the answers that I needed. I was out hiking today and had an idea for an approach that may help solve #1263, as well as generally speed things up. I need to put the idea into writing, but to realize it I will most probably need your help with the winding part. More about it shortly.
lehni
on 27 Feb 2017
I am going to keep this issue open for the remaining work on handling quality in edge case.
lehni
on 8 Mar 2017
Related issues
sansumbrella
路
3Comments
Mark12870
路
3Comments
tom10271
路
5Comments
JohnMcLear
路
5Comments
eugene-rozov
路
6Comments
Most helpful comment
I have an implementation up and running, just need to find the time to clean it up and create a PR.