CavalierContours
CavalierContours copied to clipboard
jbuckmccready
Integrating new curve types
Great library! I do have a question related to #19. How easy would it be to extend CavalierContours with new curve kinds (e.g. rational B-splines). BTW: I'm looking for a 2D geometric kernel with boolean ops and offset functionality.
I think it's possible but definitely not trivial. One of the reasons the curve input type is simple is to greatly reduce the complexity and volume of code necessary for the algorithms, that's why many libraries and academic papers focus entirely on polylines that don't even support arcs.
The main supporting functions would first need to be implemented: offsetting, finding intersects, trimming (clipping) at intersects, shortest distance check (with other B splines, lines, and arcs), bounding box, winding number, etc. I don't think any of those are too challenging, but I'd have to do some research to know how to efficiently implement the shortest distance check and winding number check (those might actually be somewhat difficult).
Once those were implemented then it would be a matter of defining a new polyline structure that supports not just the vertexes with bulges but also B splines, and then adding all of the offset and boolean ops code to support the new polyline structure with B splines using the functions described above.
If for your application it's OK to lose fidelity on the B splines then the other approach is to convert the B splines into approximating lines and arcs and then use this library, I think there are quite a few papers on how to perform that approximation (and that approximation function could even be added to this library).
I'm very busy with other things in life right now so I wouldn't count on me to fully implement it. That said I am open to reviewing pull requests and collaborating to get it integrated, and I can answer questions about the existing code.
Thank you for your extensive reply! I was of course not suggesting that you implement this on your own, just trying to estimate how feasible it would be.
We tried the arc approximation with disappointing results. B-Splines or Nurbs curve radius of curvature vary in each curve point and therefore you end up quickly with 1000-5000 arcs that slow everything down. If you keep conversion tolerance loose instead, you get fewer arcs but big gaps between arcs themselves that are even more painful.
In my opinion, the best option is the 4 control points cubic Bezier segment. It is the basic constituent of more complex B-Spline and Nurbs curves and will prevent the prohibitive increase of items to be used by booleans/offsets.
With a detailed description of what methods I need to implement/change to support bezier segments, I can contribute to this project.
Can anybody provide the checklist? Thanks.
@timoria21 I'm not entirely sure what you mean by big gaps between arcs? I haven't looked into it too much so I'm not sure how challenging the approximation problem is but my initial searching returns a lot of papers on it and it definitely should be feasible without complete loss of fidelity or an explosion in the number of segments.
Unfortunately the functions/methods are not factored around the operations required by the segments (pretty much all the code operates on a pairs of vertexes representing a segment which can be just a line or arc, see plinesegment.hpp) so probably quite a bit of structural changes are required to get everything working.
For the algorithm in this library the following is necessary (but may not be exhaustive) for any new segment type:
- Parallel offset by some value (with collapse point if it collapses), note this will require some sort of heuristic/approximation because a cubic Bezier curve cannot be parallel offset exactly into another cubic Bezier curve (code)
- Finding intersects with any other type of segment (including handling of coincident/overlapping) (code)
- Trimming/clipping at points along the segment (done in trim/join code and when creating slices from intersects)
- Trim/join to other segments for raw offset (sometimes involves clipping to a point, other times just arbitrary line/arc join that will be pruned later) (code)
- Distance check against other segments to prune invalid slices from raw offset polyline, I think for a Bezier segment this would likely end up being a general minimum distance function but I'm not sure, with only lines and arcs the code only tests end points (which may be intersect points) and midpoints (possibly between intersect points) and it works (code here and here)
- AABB (axis aligned bounding box, may also be an approximate AABB that is always at least as large as the real AABB) (code)
- Test if point "is left" query used for winding number and in boolean operations (code here and here)
Given the amount of work required to implement the the above functionality integrate it and it not being entirely clear to me that it will be any more effective than approximating to lines/arcs (given the offset step will require a heuristic approximation anyways), I think implementing a high quality bezier curve to line/arc approximation algorithm would be much easier and more effective.
Thanks for the prompt and detailed reply @jbuckmccready.
Regarding this:
I'm not entirely sure what you mean by big gaps between arcs
When you approximate a freeform curve with arcs and lines you need to accept one of the following compromises:
- many, I mean thousands, of arcs.
- gaps between a more reasonable amount of arcs
- broken tangency between a more reasonable amount of arcs
After testing spline to arcs in real-world scenarios, I realized that it is not suitable for a reliable library like this one. This is the reason why I'm thinking to try cubic bezier segments.
After checking the points you listed, I'm scared about the spline partial overlapping test and wondering if generating more items from a single bezier segment offset could be an issue.
Thanks again.
@timoria21 Thanks for explaining more of what you mean, I think I am maybe underestimating the resolution/fidelity requirements of your application.
I am curious to learn more about what types of inputs you are working with and for what application - since you mention maintaining arc tangency my guess would be CNC/tool pathing?
How many cubic bezier segments do you start with and if they were approximated as lines how many line segments would be required so the error is within tolerance for the desired resolution?
If your application is tool pathing then you could maybe get the CNC controller to perform smoothing through non-tangent arcs (up to some angle/error), or after offsetting you could approximate the lines/arcs back into B-Splines if that's the only thing the CNC controller will perform smoothing on.
In a lot of applications inputs are given as B-Splines or NURBs for compact generality, but they are just approximations in the domain and therefore approximating transforms are effective for processing.
I am doing 2D region booleans and trim (cut with open contour) with any type of curve: line/circle/arc/ellipse/elliptical arc/spline. If the spline starts tangent to a line I cannot accept a broken tangency because of conversion to arcs/lines.
One of my spline curves can be decomposed into bezier segments in a number of 10-20 items. An ellipse or elliptical arc in 2-3 bezier segments.
I aim to merge bezier segments back to a single curve once processed without any precision loss.
Ah, if you're just interested in boolean operations then it should be easier to implement. You wont need to implement the distance or offsetting function, just need to implement the following:
- Finding intersects with any other type of segment (including handling of coincident/overlapping) (code)
- Trimming/clipping at points along the segment (done in trim/join code and when creating slices from intersects)
- AABB (axis aligned bounding box, may also be an approximate AABB that is always at least as large as the real AABB) (code)
- Test if point "is left" query used for winding number and in boolean operations (code here and here)
I think you're right that the spline partial overlap detection may be challenging to implement.
I think you're right that the spline partial overlap detection may be challenging to implement.
The first idea would be to check point coincidence and tangent coincidence, pretty loose but effective.
It will take some time, I will share my results on this issue.
@timoria21 If possible can you share the cubic bezier segments decomposed from the spline curve that was not working out when you approximated as arcs (or if not that exact curve then another example)? And can you share what algorithm you are using for the approximation?
I'd like to learn and look into this more and possibly work on implementing something.
Looking at your comments again it seems one of your key requirements is avoiding loss of tangency so you can merge things back into a single continuous curve without requiring adjustments that cause precision loss?
First I wanted to contribute on a couple of bugs I found you may already know (I test thousands of real-world cases).
Bug number 1:
case Circle2Circle2IntrType::Coincident:
// determine if arcs overlap along their sweep
// start and sweep angles
auto arc1StartAndSweep = startAndSweepAngle(v1.pos(), arc1.center, v1.bulge());
// we have the arcs go the same direction to simplify checks
auto arc2StartAndSweep = [&]
{
if (v1.bulgeIsNeg() == u1.bulgeIsNeg())
{
return startAndSweepAngle(u1.pos(), arc2.center, u1.bulge());
}
return startAndSweepAngle(u2.pos(), arc2.center, -u1.bulge());
}();
// end angles (start + sweep)
auto arc1End = arc1StartAndSweep.first + arc1StartAndSweep.second;
auto arc2End = arc2StartAndSweep.first + arc2StartAndSweep.second;
/* // Missing case
if (<Both Condition A and Condition B are true>)
{
result.intrType = PlineSegIntrType::OneIntersect;
result.point1 = v1.pos();
result.point2 = arcHasSameDirection ? u1.pos() : u2.pos()
}
*/
if (std::abs(utils::deltaAngle(arc1StartAndSweep.first, arc2End)) < utils::realThreshold<Real>()) // <-- Condition A
{
// only end points touch at start of arc1
result.intrType = PlineSegIntrType::OneIntersect;
result.point1 = v1.pos();
}
else if (std::abs(utils::deltaAngle(arc2StartAndSweep.first, arc1End)) < utils::realThreshold<Real>()) // <-- Condition B
{
// only end points touch at start of arc2
result.intrType = PlineSegIntrType::OneIntersect;
result.point1 = u1.pos(); // arcHasSameDirection ? u1.pos() : u2.pos()
}
else
{
const bool arc2StartsInArc1Sweep = utils::angleIsWithinSweep(arc1StartAndSweep.first, arc1StartAndSweep.second, arc2StartAndSweep.first);
const bool arc2EndsInArc1Sweep = utils::angleIsWithinSweep(arc1StartAndSweep.first, arc1StartAndSweep.second, arc2End);
if (arc2StartsInArc1Sweep && arc2EndsInArc1Sweep)
{
// arc2 is fully overlapped by arc1
result.intrType = PlineSegIntrType::ArcOverlap;
result.point1 = u1.pos();
result.point2 = u2.pos();
}
else if (arc2StartsInArc1Sweep)
{
// overlap from arc2 start to arc1 end
result.intrType = PlineSegIntrType::ArcOverlap;
result.point1 = u1.pos(); // arcHasSameDirection ? u1.pos() : u2.pos()
result.point2 = v2.pos();
}
else if (arc2EndsInArc1Sweep)
{
// overlap from arc1 start to arc2 end
result.intrType = PlineSegIntrType::ArcOverlap;
result.point1 = v1.pos();
result.point2 = u2.pos(); // arcHasSameDirection ? u1.pos() : u2.pos()
}
else
{
const bool arc1StartsInArc2Sweep = utils::angleIsWithinSweep(arc2StartAndSweep.first, arc2StartAndSweep.second, arc1StartAndSweep.first);
if (arc1StartsInArc2Sweep)
{
result.intrType = PlineSegIntrType::ArcOverlap;
result.point1 = v1.pos();
result.point2 = v2.pos();
}
else
{
result.intrType = PlineSegIntrType::NoIntersect;
}
}
}
break;
Bug number 2:
/// Type to hold all the slices collected after slicing at intersects, slicesRemaining holds the
/// following: non-coincident slices from plineA, followed by non-coincident slices from plineB,
/// followed by coincident slices from plineA, followed by coincident slices from plineB
/// other fields hold the index markers
template <typename Real> struct CollectedSlices {
std::vector<Polyline<Real>> slicesRemaining;
std::size_t startOfPlineBSlicesIdx;
std::size_t startOfPlineACoincidentSlicesIdx;
std::size_t startOfPlineBCoincidentSlicesIdx;
};
template <typename Real, typename PlineAPointOnSlicePred, typename PlineBPointOnSlicePred>
CollectedSlices<Real> collectSlices(Polyline<Real> const &plineA, Polyline<Real> const &plineB,
ProcessForCombineResult<Real> const &combineInfo,
PlineAPointOnSlicePred &&plineAPointOnSlicePred,
PlineBPointOnSlicePred &&plineBPointOnSlicePred,
bool setOpposingOrientation) {
CollectedSlices<Real> result;
auto &slicesRemaining = result.slicesRemaining;
// slice plineA
sliceAtIntersects(plineA, combineInfo, false, plineAPointOnSlicePred, slicesRemaining);
// slice plineB
result.startOfPlineBCoincidentSlicesIdx = slicesRemaining.size(); // <-- Shouldn't it be result.startOfPlineBSlicesIdx?
sliceAtIntersects(plineB, combineInfo, true, plineBPointOnSlicePred, slicesRemaining);
// add plineA coincident slices
result.startOfPlineACoincidentSlicesIdx = slicesRemaining.size();
slicesRemaining.insert(slicesRemaining.end(), combineInfo.coincidentSlices.begin(),
combineInfo.coincidentSlices.end());
// add plineB coincident slices
std::size_t startOfPlineBCoincidentSlices = slicesRemaining.size(); // <-- Why not put into the result
slicesRemaining.insert(slicesRemaining.end(), combineInfo.coincidentSlices.begin(),
combineInfo.coincidentSlices.end());
// invert direction of plineB coincident slices to match original orientation
std::size_t coincidentSliceIdx = 0;
for (std::size_t i = startOfPlineBCoincidentSlices; i < slicesRemaining.size(); ++i) {
if (combineInfo.coincidentIsOpposingDirection[coincidentSliceIdx]) {
invertDirection(slicesRemaining[i]);
}
++coincidentSliceIdx;
}
if (setOpposingOrientation != combineInfo.plineOpposingDirections()) {
// invert plineB slice directions to match request of setOpposingOrientation
for (std::size_t i = result.startOfPlineBSlicesIdx; i < result.startOfPlineACoincidentSlicesIdx;
++i) {
invertDirection(slicesRemaining[i]);
}
for (std::size_t i = result.startOfPlineBCoincidentSlicesIdx; i < slicesRemaining.size(); ++i) {
invertDirection(slicesRemaining[i]);
}
}
return result;
}
Then I wanted to share the algorithm I am using for arc/line conversion: BsplineToArcs.pdf. Sorry, it's written in Italian.
Here is the pseudocode, I am not using C++.
private ConvertToArcApproxType Subdivide(out Curve approximation, double epsDs, double sqrEpsLc, double epsAt)
{
approximation = null;
double ls01 = length();
Point startPoint = StartPoint;
Point endPoint = EndPoint;
Line candidateSegment = new Line(startPoint, endPoint);
if (ls01 < epsDs)
{
approximation = candidateSegment;
return ConvertToArcApproxType.Segment;
}
Point Pm = PointAt(Domain.Mid);
// Trying to approximate with a segment
Vector startTangent = StartTangent;
Vector endTangent = EndTangent;
if (AreParallel(startTangent, endTangent))
{
if (ls01 - Distance(startPoint, endPoint) < epsDs)
{
if (DistanceSquared(Pm, candidateSegment.MidPoint) < sqrEpsLc)
{
approximation = candidateSegment;
return ConvertToArcApproxType.Segment;
}
}
return ConvertToArcApproxType.None;
}
// Trying to approximate with an arc
Vector t0 = (Vector) startTangent.Clone();
Vector t1 = (Vector) endTangent.Clone();
t0.Normalize();
t1.Normalize();
Plane pArc = new Plane(startPoint, Pm, endPoint);
Line tmpPt0 = new Line(Point.Origin, Vector.Cross(t0, pArc.AxisZ).AsPoint);
Line tmpPt1 = new Line(Point.Origin, Vector.Cross(t1, pArc.AxisZ).AsPoint);
tmpPt0.TransformBy(new Translation(startPoint.AsVector));
tmpPt1.TransformBy(new Translation(endPoint.AsVector));
Segment pT0 = new Segment(tmpPt0.StartPoint, tmpPt0.EndPoint);
Segment pT1 = new Segment(tmpPt1.StartPoint, tmpPt1.EndPoint);
if (Intersection(pT0, pT1, true, out Point arcCenter, out _))
{
if (arcCenter != startPoint && arcCenter != endPoint)
{
Arc candidateArc = new Arc(arcCenter, startPoint, endPoint);
if (DistanceSquared(Pm, candidateArc.MidPoint) < sqrEpsLc && (ls01 - candidateArc.Length()) < epsDs)
{
if (Vector.AreCoincident(candidateArc.EndTangent, endTangent, epsAt))
{
approximation = candidateArc;
return ConvertToArcApproxType.Arc;
}
}
}
}
return ConvertToArcApproxType.None;
}
Last and not least, some interesting articles on the bezier to arcs topic I found online:
- https://dlacko.org/blog/2016/10/19/approximating-bezier-curves-by-biarcs/
- https://github.com/domoszlai/bezier2biarc
Regarding this question:
Looking at your comments again it seems one of your key requirements is avoiding loss of tangency so you can merge things back into a single continuous curve without adjustments causing precision loss?
Yes, I aim to avoid converting splines to arc/lines do booleans, and convert back. I already have the proof that to get the required precision you need to tight the tolerance and deal with 10K items for a few spline curves. This. obviously. badly affects the library speed.
@timoria21 Thanks for the bug reports and information.
Both of those bugs are fixed in my continued work on this project in Rust here: https://github.com/jbuckmccready/cavalier_contours. The Rust code has the same performance or better and is written in 100% safe code. The Rust code has a lot more test coverage and code documentation, and I think a few other bugs were fixed along the way. Going forward I'm only going to be working on the Rust code, but the algorithm development is useful regardless of language. I will archive the C++ project to make it clear beyond just the readme that it is not being actively maintained (I will unarchive it if others want to contribute code changes to keep it alive).
I don't know what language you're working in but I created C bindings in the Rust project (which I use from C#) and python bindings could also be easy to create using PyO3.
Rust also compiles to wasm nicely and I use a web app for visual testing and demo:
- Interactive hosted page here: https://jbuckmccready.github.io/cavalier_contours_web_demo_page
- Interactive page source code here: https://github.com/jbuckmccready/cavalier_contours_web_demo
I haven't yet added multipolyline/island offsetting algorithm since I had some ideas for creating an abstract shape type for it but, there are some additional functions I added that are not in the C++ as well.
I have been somewhat absent from this project for the last few months but I am still actively interested in it, I will update the Rust project readme to better convey the state of it.
I also found this resource for bezier curve operations (including approximating as arcs): http://pomax.github.io/bezierjs/
EDIT: I will hold off on archiving this repository since that makes all the issues read only as well...
I think a few other bugs were fixed along the way
Did you remember which ones?
I also found this resource for bezier curve operations (including approximating as arcs)
Cool, please remember that the way to go is to support bezier, not to accept gaps or thousands of items with arcs approximations. I did a lot of tests that confirm this.
Now the question is: will you address bezier any time soon or should I try by myself first?
Did you remember which ones?
I remember finding issues for boolean operations between alternating winding/direction of polylines that mostly overlap (where only 1 or 2 segments deviated, the vertexes overlapped, and then alternating winding direction). I created a large set of tests to verify everything and that's how I found the issues, the Rust code has all the tests here.
Cool, please remember that the way to go is to support bezier, not to accept gaps or thousands of items with arcs approximations. I did a lot of tests that confirm this.
Now the question is: will you address bezier any time soon or should I try by myself first?
I wont be working on the bezier segments any time soon.
If I can find some time I'd like to add the bezier to arc approximation to the library for convenience for cases where the tradeoff is acceptable.
I also want to look into ways of efficiently approximating to arcs with a possible loss in continuity (maybe drops to just G0 positional) to perform a parallel offset efficiently, and then morphs the result back into a B-spline with higher continuity (G1 tangency or G2 curvature) in a controlled approximating way that minimizes the number of segments - this is maybe a really hard problem and maybe you've already looked into it but I was going to do some research.
In the past I briefly looked into how to robustly parallel offset B-Splines (including cases of open curves, self intersecting curves, and offsets that result in collapsed global sub regions) it seemed very difficult to do robustly and if possible seemed computationally inefficient.
I remember finding issues for boolean operations between alternating winding/direction of polylines that mostly overlap (where only 1 or 2 segments deviated, the vertexes overlapped, and then alternating winding direction). I created a large set of tests to verify everything and that's how I found the issues, the Rust code has all the tests here.
It would be nice to have the changeset numbers.
I wont be working on the bezier segments any time soon.
Ok, I will be the first trying. I will share my results here. I may need more help when I dig into the code you referenced here:
Ah, if you're just interested in boolean operations then it should be easier to implement. You wont need to implement the distance or offsetting function, just need to implement the following:
- Finding intersects with any other type of segment (including handling of coincident/overlapping) (code)
- Trimming/clipping at points along the segment (done in trim/join code and when creating slices from intersects)
- AABB (axis aligned bounding box, may also be an approximate AABB that is always at least as large as the real AABB) (code)
- Test if point "is left" query used for winding number and in boolean operations (code here and here)
Hello,
I'm getting in trouble with the sliceAtIntersects() function and closed splines. I define closed splines as a single vertex in polylines but when for some reason they are split into one or multiple points I probably need to add one vertex at the end (when you split a closed spline define by one vertex you need two segments defined by three vertices).
Can you help me to understand where I should add this third vertex inside the sliceAtIntersect function?
Thanks.
Hello,
I'm getting in trouble with the
sliceAtIntersects()function and closed splines. I define closed splines as a single vertex in polylines but when for some reason they are split into one or multiple points I probably need to add one vertex at the end (when you split a closed spline define by one vertex you need two segments defined by three vertices).Can you help me to understand where I should add this third vertex inside the
sliceAtIntersectfunction?Thanks.
The sliceAtIntersects() function takes in all the intersects and cuts the polyline at the intersect points to return a new set of polylines that all start and end at sequential pairs of intersect points. I'm not as familiar with cutting/splitting splines but it seems like you would iterate through the intersects in order (the same way the function works now) and when a slice polyline is to be formed between two intersect points you would add any additional vertexes required to form the geometric "slice" as a sequence of vertexes. This is done now (additional vertexes added) in the case where there are multiple intersects along the same line or arc segment.
Described in another way:
The main operation is the following: as inputs you have two intersect points (and the associated segment indexes they lie on) and the source polyline, from that you need to form a new polyline that represents all the segments (or sub segment if both intersects lie on the same segment) between the two intersect points along the source polyline.
The above operation is performed on all sequential pairs (ordered along the length of the polyline) of intersect points to return all the "slices" as new polylines.
Hopefully that explanation helps.
you would iterate through the intersects in order
Yes, this was one of the problems. Thank you again for the detailed explanation.
I am doing good progress on splines. I'm currently dealing with more than two intersections (spline can intersect line/arc many times) and the following "is point inside" bugs, do you know them? Did you already fix them? Can you remember the changeset? I can provide the exact point coords if necessary.
I do not recall exact changes - and they may not even be fixed in the Rust code. If you post the coordinates for point and polylines I will look into it. It looks like it has to do with the point being inline with multiple segment's directional vectors.
Here you go, thanks.
Polyline p5 = new Polyline();
Polyline p6 = new Polyline();
p5.addVertex(new PlineVertex(-10, 0, 1));
p5.addVertex(new PlineVertex(10, 0, 0));
p5.addVertex(new PlineVertex(20, 0, 0));
p5.addVertex(new PlineVertex(20, -10, 0));
p5.addVertex(new PlineVertex(-20, -10, 0));
p5.addVertex(new PlineVertex(-20, 0, 0));
p5.iSClosed = true;
Vector v5 = new Vector(0, 0);
p6.addVertex(-5.51073e-15, -30, 0.269712);
p6.addVertex(26.0788, -14.8288, 0);
p6.addVertex(76.0788, 73.104, 0.12998);
p6.addVertex(80, 87.9329, 0);
p6.addVertex(80, 130, 0);
p6.addVertex(50, 130, 0);
p6.addVertex(50, 95, -0.414214);
p6.addVertex(40, 85, 0);
p6.addVertex(0, 85, 0);
p6.iSClosed = true;
Vector v6 = new Vector(-20, 85);
I just tested now in both Rust code and C++, the winding number returned is 0 (point is outside polyline) which is correct. Is there something I am missing?
I'm terribly sorry, my fault. As always happens with deep modifications I left something behind. Please consider removing my request and your answer (two posts above) as they don't add anything to the discussion.
I will soon paste here all the data structures I'm using to handle Nurbs curves. I have not finished yet but some pieces of code are good enough to be posted.
Thanks again.
I was finally able to support not only Bezier segments but also Nurbs curves in both booleans and offset. In the following days, I will start sharing some code/screenshot and explain all the issues/limitations I faced. The good news is that extending this library for spline/bezier/nurbs is possible and it works correctly on a large number of cases.
Let's start from here. I will split my discussion into several posts.
I simplify linear and circular Nurbs curves before passing them to the CavContour (CC) library. Linear Nurbs curves can be identified by degree 1 and 2 control points or by checking the control point distance from the line connecting the first and last vertex. Circular Nurbs curves can be simplified if the control points are in some special configurations like 90, 180, 270, 360 degrees. Actually, this can be further improved, see "Computing offsets of NURBS curves and surfaces" Les A. Piegl, Wayne Tiller.
Here you see how I have changed the PlineVertex class, adding: SegmentType, Curve and BezSegs properties. Basically, what I wanted to do, was to split complex Nurbs curves in bezier segments and with a simple "belly" (no inflection point) to be able to check the winding number in the same way CC does for arcs.
class PlineVertex
{
...
/// <summary>
/// Actual nurbs curve data.
/// </summary>
public Curve Curve { get; set; }
private Curve[] _bezSegs;
public Curve[] BezSegs
{
get
{
if (_bezSegs == null && Curve != null)
_bezSegs = DecomposeSegments();
return _bezSegs;
}
}
/// <summary>
/// Gets the type of the segment. Should replace the bulgeIsZero, bulgeIsPos and bulgeIsNeg
/// </summary>
public int SegmentType
{
get
{
if (Curve != null)
{
return 2;
}
if (bulgeIsZero())
{
return 0;
}
return 1;
}
}
private Curve[] DecomposeSegments()
{
Curve[] decomposedSegments = Curve.Decompose();
List<ICurve> result = new List<ICurve>(decomposedSegments.Length);
foreach (Curve seg in decomposedSegments)
{
double tolLinear = 1e-6;
if(seg.IsLinear(tolLinear, out _))
{
result.Add(new Line(seg.ControlPoints.First(), seg.ControlPoints.Last()).ConvertToNurbs());
}
else if (seg.IsPoint)
{
continue;
}
else
{
double[] ts = GetInflectionParams(seg);
if (ts != null && ts.Length > 0)
{
List<ICurve> splits = new List<ICurve>();
if(Utility.SplitAtParamArray(seg, ts, splits))
result.AddRange(splits);
else
result.Add(seg);
}
else
{
result.Add(seg);
}
}
}
return result.Cast<Curve>().ToArray();
}
private static double[] GetInflectionParams(Curve bezier)
{
if (bezier.Degree == 2)
return null;
if (bezier.Degree != 3)
return null;
if (bezier.IsRational)
throw new NotImplementedException();
return compute_Bezier_degree_3_inflections(bezier.ControlPoints[0].X, bezier.ControlPoints[0].Y,
bezier.ControlPoints[1].X, bezier.ControlPoints[1].Y,
bezier.ControlPoints[2].X, bezier.ControlPoints[2].Y,
bezier.ControlPoints[3].X, bezier.ControlPoints[3].Y);
}
/// <summary>
/// From https://github.com/qnzhou/nanospline
/// </summary>
private static double[] compute_Bezier_degree_3_inflections(double cx0, double cy0, double cx1, double cy1, double cx2, double cy2, double cx3, double cy3, double t0 = 0, double t1 = 1)
{
List<double> result = new List<double>();
double tol = 1e-9;
find_real_roots_in_interval_2(new double[]{
-18 * cx0 * cy1 + 18 * cx0 * cy2 + 18 * cx1 * cy0 - 18 * cx1 * cy2 - 18 * cx2 * cy0 + 18 * cx2 * cy1,
36 * cx0 * cy1 - 54 * cx0 * cy2 + 18 * cx0 * cy3 - 36 * cx1 * cy0 + 54 * cx1 * cy2 - 18 * cx1 * cy3 + 54 * cx2 * cy0 - 54 * cx2 * cy1 - 18 * cx3 * cy0 + 18 * cx3 * cy1,
-18 * cx0 * cy1 + 36 * cx0 * cy2 - 18 * cx0 * cy3 + 18 * cx1 * cy0 - 54 * cx1 * cy2 + 36 * cx1 * cy3 - 36 * cx2 * cy0 + 54 * cx2 * cy1 - 18 * cx2 * cy3 + 18 * cx3 * cy0 - 36 * cx3 * cy1 + 18 * cx3 * cy2
}, result, t0, t1, tol);
return result.ToArray();
}
private static void find_real_roots_in_interval_2(double[] coeffs, List<double> roots, double t0, double t1, double eps)
{
Debug.Assert(coeffs.Length > 2);
// Largest degree is zero, the polynomial is one degree less
if (Math.Abs(coeffs[2]) < eps)
{
find_real_roots_in_interval_1(coeffs, roots, t0, t1, eps);
return;
}
double a = coeffs[2];
double b = coeffs[1];
double c = coeffs[0];
double discr = b * b - 4 * a * c;
// no real root
if (discr < 0) return;
double root;
// duplicate root
if (Math.Abs(discr) < eps)
{
root = -b / (2 * a);
if (root >= t0 && root <= t1) roots.Add(root);
return;
}
double sqrt_discr = Math.Sqrt(discr);
root = (-b - sqrt_discr) / (2 * a);
if (root >= t0 && root <= t1) roots.Add(root);
root = (-b + sqrt_discr) / (2 * a);
if (root >= t0 && root <= t1) roots.Add(root);
}
private static void find_real_roots_in_interval_1(double[] coeffs, List<double> roots, double t0, double t1, double eps)
{
Debug.Assert(coeffs.Length > 1);
// Largest degree is zero, the polynomial is one degree less
if (Math.Abs(coeffs[1]) < eps)
{
find_real_roots_in_interval_0(coeffs, roots, t0, t1, eps);
return;
}
double a = coeffs[1];
double b = coeffs[0];
double root = -b / a;
if (root >= t0 && root <= t1) roots.Add(root);
}
static void find_real_roots_in_interval_0(double[] coeffs, List<double> roots, double t0, double t1, double eps)
{
Debug.Assert(coeffs.Length > 0);
if (Math.Abs(coeffs[0]) < eps) throw new Exception("Infinite root error");
}
Below, the changes to support Curve overlap and more than two intersections:
/// Split the segment defined by v1 to v2 at some point defined along it.
enum PlineSegIntrType
{
NoIntersect,
TangentIntersect,
OneIntersect,
TwoIntersects,
SegmentOverlap,
ArcOverlap,
CurveOverlap,
CurveMoreThan2Intersects
}
class IntrPlineSegsResult
{
public PlineSegIntrType intrType;
public Vector point1 = new Vector(2);
public Vector point2 = new Vector(2);
// Used to store more than 2 intersections on Curves.
public Point3D[] intersections;
}
These are the changes to splitAtPoint() method:
/// Result of splitting a segment v1 to v2.
static SplitResult splitAtPoint(PlineVertex v1, PlineVertex v2, Vector point, Vector limitSplitCurveToPoint = null)
{
SplitResult result = new SplitResult();
if (Vector.fuzzyEqual(v1.Pos, point))
{
result.updatedStart = new PlineVertex(point, 0);
result.splitVertex = new PlineVertex(point, v1.Bulge) { Curve = v1.Curve };
}
else if ((v1.SegmentType != 2) && Vector.fuzzyEqual(v1.Pos, v2.Pos))
{
result.updatedStart = new PlineVertex(point, 0);
result.splitVertex = new PlineVertex(point, v1.Bulge) { Curve = v1.Curve };
}
else if (Vector.fuzzyEqual(v2.Pos, point))
{
result.updatedStart = v1;
result.splitVertex = new PlineVertex(v2.Pos, 0);
}
else if (v1.SegmentType == 0 && v1.Curve == null)
{
result.updatedStart = v1;
result.splitVertex = new PlineVertex(point, 0);
}
else if (v1.SegmentType == 1)
{
ArcRadiusAndCenter radiusAndCenter = arcRadiusAndCenter(v1, v2);
Vector arcCenter = radiusAndCenter.center;
double a = Vector2.angle(arcCenter, point);
double arcStartAngle = Vector2.angle(arcCenter, v1.Pos);
double theta1 = utils.mathutils.deltaAngle(arcStartAngle, a);
double bulge1 = Math.Tan(theta1 / 4);
double arcEndAngle = Vector2.angle(arcCenter, v2.Pos);
double theta2 = utils.mathutils.deltaAngle(a, arcEndAngle);
double bulge2 = Math.Tan(theta2 / 4);
result.updatedStart = new PlineVertex(v1.Pos, bulge1);
result.splitVertex = new PlineVertex(point, bulge2);
}
else
{
Point3D splitPt = new Point3D(point.X, point.Y);
Point4D[] ctrlPts = v1.Curve.ControlPoints;
// A larger tolerance is used to check start/end coincidence for nurbs
double thrCoincidence = mathutils.realThreshold;
if (Point3D.DistanceSquared((Point3D)ctrlPts[0], splitPt) < thrCoincidence)
{
result.updatedStart = new PlineVertex(v1.Pos, 0);
result.splitVertex = new PlineVertex(point, 0) { Curve = (Curve)v1.Curve.Clone() };
}
else if (Point3D.DistanceSquared((Point3D)ctrlPts[ctrlPts.Length - 1], splitPt) < thrCoincidence)
{
result.updatedStart = new PlineVertex(v1.Pos, 0) { Curve = (Curve)v1.Curve.Clone() };
result.splitVertex = new PlineVertex(point, v2.Bulge);
}
else if (v1.Curve.SplitBy(splitPt, out ICurve lower, out ICurve upper))
{
result.updatedStart = new PlineVertex(v1.Pos, 0) { Curve = (Curve)lower };
result.splitVertex = new PlineVertex(point, 0) { Curve = (Curve)upper };
}
}
// To manage the situation of overlapping curves, that must be trimmed: lines and arcs doesn't have this problem, since every successor in the polyline ends the predecessor
if (limitSplitCurveToPoint != null)
{
LimitPlineVertexCurve(limitSplitCurveToPoint, result.splitVertex);
}
return result;
}
Then the limitPlineVertexCurve() method definition:
/// <summary>
/// Limits Curve of pline vertex to given vector.
/// </summary>
/// <param name="limitPt">The limiting point.</param>
/// <param name="pV">The pline vertex to limit.</param>
static void limitPlineVertexCurve(Vector limitPt, PlineVertex pV)
{
if (pV.Curve == null)
throw new EyeshotException("Limiting can be done only on a splitVertex with Curve.");
if (pV.Curve.SplitBy(new Point3D(limitPt.X, limitPt.Y), out ICurve left, out _))
{
pV.Curve = (Curve) left;
}
}