Integrating new curve types

[adam-urbanczyk]

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.

[jbuckmccready]

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.

[adam-urbanczyk]

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.

[timoria21]

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.

[jbuckmccready]

@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.

[timoria21]

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:

1. many, I mean thousands, of arcs.
2. gaps between a more reasonable amount of arcs
3. 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.

[jbuckmccready]

@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.

[timoria21]

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.

[jbuckmccready]

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.

[timoria21]

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.

[jbuckmccready]

@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?

[timoria21]

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

[timoria21]

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

[timoria21]

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

[timoria21]

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.

[jbuckmccready]

@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...

[timoria21]

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?

[jbuckmccready]

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.

[timoria21]

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)

[timoria21]

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.

[jbuckmccready]

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.

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.

[timoria21]

you would iterate through the intersects in order

Yes, this was one of the problems. Thank you again for the detailed explanation.

[timoria21]

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.

[jbuckmccready]

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.

[timoria21]

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

[jbuckmccready]

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?

[timoria21]

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.

[timoria21]

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.

[jbuckmccready]

@timoria21 That sounds awesome, I look forward to checking it out.

[timoria21]

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

[timoria21]

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

[timoria21]

Then changes to segTangentVector() method:

        static Vector segTangentVector(PlineVertex v1, PlineVertex v2, Vector pointOnSeg)
        {
            if (v1.SegmentType == 2)
            {
                Vector3D tan = v1.Curve.StartTangent;

                return new Vector(tan.X, tan.Y);
            }

            if (v1.SegmentType == 0)
            {
                return Vector.operatorMinus(v2.Pos, v1.Pos);
            }

            ArcRadiusAndCenter arc = arcRadiusAndCenter(v1, v2);
            if (v1.bulgeIsPos())
            {
                // ccw, rotate vector from center to pointOnCurve 90 degrees
                return new Vector(-(pointOnSeg.Y - arc.center.Y), pointOnSeg.X - arc.center.X);
            }

            // cw, rotate vector from center to pointOnCurve -90 degrees
            return new Vector(pointOnSeg.Y - arc.center.Y, -(pointOnSeg.X - arc.center.X));
        }

Then changed to closestPointOnSeg() method:

        /// Compute the closest point on a segment defined by v1 to v2 to the point given.
        static Vector closestPointOnSeg(PlineVertex v1, PlineVertex v2, Vector point)
        {
            if (v1.SegmentType == 2)
            {
                v1.Curve.ClosestPointTo(new Point3D(point.X, point.Y), out double t);

                Point3D closest = v1.Curve.PointAt(t);

                return new Vector(closest.X, closest.Y);
            }

            if (v1.SegmentType == 0)
            {
                return Vector2.closestPointOnLineSeg(v1.Pos, v2.Pos, point);
            }

            ArcRadiusAndCenter arc = plineSegmentGlobalMembers.arcRadiusAndCenter(v1, v2);

            if (Vector.fuzzyEqual(point, arc.center))
            {
                // avoid normalizing zero length vector (point is at center, just return start point)
                return v1.Pos;
            }

            if (Vector2.pointWithinArcSweepAngle(arc.center, v1.Pos, v2.Pos, v1.Bulge, point))
            {
                // closest point is on the arc
                Vector vToPoint = Vector.operatorMinus(point, arc.center);
                Vector.normalize(vToPoint);
                return Vector.operatorPlus(Vector.operatorInto(vToPoint, arc.radius), arc.center);//return arc.radius * vToPoint + arc.center;
            }

            // else closest point is one of the ends
            double dist1 = Vector2.distSquared(v1.Pos, point);
            double dist2 = Vector2.distSquared(v2.Pos, point);
            if (dist1 < dist2)
            {
                return v1.Pos;
            }

            return v2.Pos;
        }

[timoria21]

Then changed to createFastApproxBoundingBox() method:

        /// Computes a fast approximate AABB of a segment described by v1 to v2, bounding box may be larger
        /// than the true bounding box for the segment
        static AABB createFastApproxBoundingBox(PlineVertex v1, PlineVertex v2)
        {
            AABB result = new AABB();

            if (v1.SegmentType == 2)
            {
                v1.Curve.ControlBoundingBox(out Point3D min, out Point3D max);

                result.xMin = min.X;
                result.yMin = min.Y;
                result.xMax = max.X;
                result.yMax = max.Y;

                return result;
            }

            if (v1.SegmentType == 0)
            {
                if (v1.X < v2.X)
                {
                    result.xMin = v1.X;
                    result.xMax = v2.X;
                }
                else
                {
                    result.xMin = v2.X;
                    result.xMax = v1.X;
                }

                if (v1.Y < v2.Y)
                {
                    result.yMin = v1.Y;
                    result.yMax = v2.Y;
                }
                else
                {
                    result.yMin = v2.Y;
                    result.yMax = v1.Y;
                }

                return result;
            }

            // For arcs we don't compute the actual extents which is slower, instead we create an approximate
            // bounding box from the rectangle formed by extending the chord by the sagitta, NOTE: this
            // approximate bounding box is always equal to or bigger than the true bounding box
            double b = v1.Bulge;
            double offsX = b * (v2.Y - v1.Y) / 2;
            double offsY = -b * (v2.X - v1.X) / 2;

            double pt1X = v1.X + offsX;
            double pt2X = v2.X + offsX;
            double pt1Y = v1.Y + offsY;
            double pt2Y = v2.Y + offsY;

            double endPointXMin;
            double endPointXMax;
            if (v1.X < v2.X)
            {
                endPointXMin = v1.X;
                endPointXMax = v2.X;
            }
            else
            {
                endPointXMin = v2.X;
                endPointXMax = v1.X;
            }

            double ptXMin;
            double ptXMax;
            if (pt1X < pt2X)
            {
                ptXMin = pt1X;
                ptXMax = pt2X;
            }
            else
            {
                ptXMin = pt2X;
                ptXMax = pt1X;
            }

            double endPointYMin;
            double endPointYMax;
            if (v1.Y < v2.Y)
            {
                endPointYMin = v1.Y;
                endPointYMax = v2.Y;
            }
            else
            {
                endPointYMin = v2.Y;
                endPointYMax = v1.Y;
            }

            double ptYMin;
            double ptYMax;
            if (pt1Y < pt2Y)
            {
                ptYMin = pt1Y;
                ptYMax = pt2Y;
            }
            else
            {
                ptYMin = pt2Y;
                ptYMax = pt1Y;
            }

            result.xMin = Math.Min(endPointXMin, ptXMin);
            result.yMin = Math.Min(endPointYMin, ptYMin);
            result.xMax = Math.Max(endPointXMax, ptXMax);
            result.yMax = Math.Max(endPointYMax, ptYMax);

            return result;
        }

Then segLength() method changes:

        /// Calculate the path length for the segment defined from v1 to v2.
        static double segLength(PlineVertex v1, PlineVertex v2)
        {
            if (Vector.fuzzyEqual(v1.Pos, v2.Pos))
            {
                return 0;
            }

            if (v1.SegmentType == 2)
            {
                return v1.Curve.Length();
            }

            if (v1.SegmentType == 0)
            {
                return Math.Sqrt(Vector2.distSquared(v1.Pos, v2.Pos));
            }

            ArcRadiusAndCenter arc = plineSegmentGlobalMembers.arcRadiusAndCenter(v1, v2);
            double startAngle = Vector2.angle(arc.center, v1.Pos);
            double endAngle = Vector2.angle(arc.center, v2.Pos);
            return Math.Abs(arc.radius * utils.mathutils.deltaAngle(startAngle, endAngle));
        }

Then segMidpoint() method changes:

        /// Return the mid point along a segment path.
        static Vector segMidpoint(PlineVertex v1, PlineVertex v2)
        {
            if (v1.SegmentType == 2)
            {
                double t = v1.Curve.Domain.ParameterAt(0.5);

                Point3D pt = v1.Curve.PointAt(t);

                return new Vector(pt.X, pt.Y);
            }

            if (v1.SegmentType == 0)
            {
                return Vector2.midpoint(v1.Pos, v2.Pos);
            }

            ArcRadiusAndCenter arc = arcRadiusAndCenter(v1, v2);
            double a1 = Vector2.angle(arc.center, v1.Pos);
            double a2 = Vector2.angle(arc.center, v2.Pos);
            double angleOffset = Math.Abs(utils.mathutils.deltaAngle(a1, a2) / 2);
            // use arc direction to determine offset sign to robustly handle half circles
            double midAngle = v1.bulgeIsPos() ? a1 + angleOffset : a1 - angleOffset;
            return Vector2.pointOnCircle(arc.radius, arc.center, midAngle);
        }

[timoria21]

Then intrPlineSegs() method changes:

        static IntrPlineSegsResult intrPlineSegs(PlineVertex v1, PlineVertex v2, PlineVertex u1, PlineVertex u2)
        {
            IntrPlineSegsResult result = new IntrPlineSegsResult();
            bool vIsLine = v1.SegmentType == 0;
            bool uIsLine = u1.SegmentType == 0;
            bool vIsArc = v1.SegmentType == 1;
            bool uIsArc = u1.SegmentType == 1;

            if (vIsLine && uIsLine)
            {
                IntrLineSeg2LineSeg2Result intrResult = IntrLineSeg2LineSeg2.intrLineSeg2LineSeg2(v1.Pos, v2.Pos, u1.Pos, u2.Pos);
                switch (intrResult.intrType)
                {
                    case LineSeg2LineSeg2IntrType.None:
                        result.intrType = PlineSegIntrType.NoIntersect;
                        break;
                    case LineSeg2LineSeg2IntrType.True:
                        result.intrType = PlineSegIntrType.OneIntersect;
                        result.point1 = intrResult.point;
                        break;
                    case LineSeg2LineSeg2IntrType.Coincident:
                        result.intrType = PlineSegIntrType.SegmentOverlap;
                        // build points from parametric parameters (using second segment as defined by the function)
                        result.point1 = Vector2.pointFromParametric(u1.Pos, u2.Pos, intrResult.t0);
                        result.point2 = Vector2.pointFromParametric(u1.Pos, u2.Pos, intrResult.t1);
                        break;
                    case LineSeg2LineSeg2IntrType.False:
                        result.intrType = PlineSegIntrType.NoIntersect;
                        break;
                }
            }

            else if (vIsLine && uIsArc)
            {
                processLineArcIntr(v1.Pos, v2.Pos, u1, u2, result);
            }

            else if (uIsLine && vIsArc)
            {
                processLineArcIntr(u1.Pos, u2.Pos, v1, v2, result);
            }

            else if (uIsArc && vIsArc)
            {
                ArcRadiusAndCenter arc1 = arcRadiusAndCenter(v1, v2);
                ArcRadiusAndCenter arc2 = arcRadiusAndCenter(u1, u2);

                IntrCircle2Circle2Result intrResult = intrcircle2circle2.intrCircle2Circle2(arc1.radius, arc1.center, arc2.radius, arc2.center);

                switch (intrResult.intrType)
                {
                    case Circle2Circle2IntrType.NoIntersect:
                        result.intrType = PlineSegIntrType.NoIntersect;
                        break;
                    case Circle2Circle2IntrType.OneIntersect:
                        if (bothArcsSweepPoint(intrResult.point1, v1, v2, u1, u2, arc1, arc2))
                        {
                            result.intrType = PlineSegIntrType.OneIntersect;
                            result.point1 = intrResult.point1;
                        }
                        else
                        {
                            result.intrType = PlineSegIntrType.NoIntersect;
                        }
                        break;
                    case Circle2Circle2IntrType.TwoIntersects:
                        {
                            bool pt1InSweep = bothArcsSweepPoint(intrResult.point1, v1, v2, u1, u2, arc1, arc2);
                            bool pt2InSweep = bothArcsSweepPoint(intrResult.point2, v1, v2, u1, u2, arc1, arc2);
                            if (pt1InSweep && pt2InSweep)
                            {
                                result.intrType = PlineSegIntrType.TwoIntersects;
                                result.point1 = intrResult.point1;
                                result.point2 = intrResult.point2;
                            }
                            else if (pt1InSweep)
                            {
                                result.intrType = PlineSegIntrType.OneIntersect;
                                result.point1 = intrResult.point1;
                            }
                            else if (pt2InSweep)
                            {
                                result.intrType = PlineSegIntrType.OneIntersect;
                                result.point1 = intrResult.point2;
                            }
                            else
                            {
                                result.intrType = PlineSegIntrType.NoIntersect;
                            }
                        }
                        break;
                    case Circle2Circle2IntrType.Coincident:
                        // determine if arcs overlap along their sweep
                        // start and sweep angles
                        Tuple<double, double> arc1StartAndSweep = startAndSweepAngle(v1.Pos, arc1.center, v1.Bulge);
                        Tuple<double, double> arc2StartSweep = arc2StartAndSweep(v1, u1, u2, arc2, out bool sameDirection);

                        // end angles (start + sweep)
                        double arc1End = arc1StartAndSweep.Item1 + arc1StartAndSweep.Item2;
                        double arc2End = arc2StartSweep.Item1 + arc2StartSweep.Item2;//check here

                        bool endPoint2TouchAtStartArc1 = Math.Abs(utils.mathutils.deltaAngle(arc1StartAndSweep.Item1, arc2End)) < utils.mathutils.realThreshold;
                        bool endPoint1TouchAtStartArc2 = Math.Abs(utils.mathutils.deltaAngle(arc2StartSweep.Item1, arc1End)) < utils.mathutils.realThreshold;
                        if (endPoint2TouchAtStartArc1 && endPoint1TouchAtStartArc2)
                        {
                            // only end points touch at start of arc1
                            result.intrType = PlineSegIntrType.TwoIntersects;
                            result.point1 = v1.Pos;
                            result.point2 = sameDirection ? u1.Pos : u2.Pos; // TE4_CylSeam_HalfSlot
                        }
                        else if (endPoint2TouchAtStartArc1)
                        {
                            // only end points touch at start of arc1
                            result.intrType = PlineSegIntrType.OneIntersect;
                            result.point1 = v1.Pos;
                        }
                        else if (endPoint1TouchAtStartArc2)
                        {
                            // only end points touch at start of arc2
                            result.intrType = PlineSegIntrType.OneIntersect;
                            result.point1 = sameDirection ? u1.Pos : u2.Pos; // TE4_CylSeam_HalfSlot
                        }
                        else
                        {
                            bool arc2StartsInArc1Sweep = utils.mathutils.angleIsWithinSweep(arc1StartAndSweep.Item1, arc1StartAndSweep.Item2, arc2StartSweep.Item1);
                            bool arc2EndsInArc1Sweep = utils.mathutils.angleIsWithinSweep(arc1StartAndSweep.Item1, arc1StartAndSweep.Item2, 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 = sameDirection ? 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 = sameDirection ? u2.Pos : u1.Pos;
                            }
                            else
                            {
                                bool arc1StartsInArc2Sweep = utils.mathutils.angleIsWithinSweep(arc2StartSweep.Item1, arc2StartSweep.Item2, arc1StartAndSweep.Item1);
                                if (arc1StartsInArc2Sweep)
                                {
                                    result.intrType = PlineSegIntrType.ArcOverlap;
                                    result.point1 = v1.Pos;
                                    result.point2 = v2.Pos;
                                }
                                else
                                {
                                    result.intrType = PlineSegIntrType.NoIntersect;
                                }
                            }
                        }

                        break;
                }
            }
            else if (uIsArc) // v1 is curve
            {
                Curve crv1 = v1.Curve;
                Arc crv2 = ConvertSegmentToArc(u1, u2);

                FindIntersection(crv1, crv2, result, v1, v2, u1, u2);
            }
            else if (vIsArc) // u1 is curve
            {
                Arc crv1 = ConvertSegmentToArc(v1, v2);
                Curve crv2 = u1.Curve;

                FindIntersection(crv1, crv2, result, v1, v2, u1, u2);
            }
            else if (uIsLine) // v1 is curve
            {
                Curve crv1 = v1.Curve;
                Line crv2 = new Line(u1.X, u1.Y, u2.X, u2.Y);

                FindIntersection(crv1, crv2, result, v1, v2, u1, u2);
            }
            else if (vIsLine) // u1 is curve
            {
                Line crv1 = new Line(v1.X, v1.Y, v2.X, v2.Y);
                Curve crv2 = u1.Curve;

                FindIntersection(crv1, crv2, result, v1, v2, u1, u2);
            }
            else // both u1 and v1 are curve
            {
                FindIntersection(v1.Curve, u1.Curve, result, v1, v2, u1, u2);
            }

            return result;
        }

[timoria21]

That's enough for today, I will continue to post changes in the following days.

Meanwhile, I attach one of my initial test cases. It includes one ellipse curve, that can be easily converted to Nurbs.

[jbuckmccready]

Nice seeing it all come together. The code looks pretty clean due to how you were able to tackle each sub operation to get the boolean operation working. I like the approach of decomposing into bezier segments without inflection points for the winding number check.

It looks like you're working in C#? I'm curious to know what your use case for the operations are, or if it's purely academic/exploratory? If you are willing/can you should make a public git repository for the code, it's a very cool project!

Thanks for sharing your work and progress.

I have made some minor API improvements, added a few tests, and fixed a few edge case bugs in the Rust repository but unfortunately have not found the time to build out more on the algorithm side, it's inspiring to see what you're doing.

[timoria21]

No. I cannot make a public repository for the code, sorry. I cannot share the function from the book above. Yes, It would be nice to know what improvements (changesets) you did on the Rust project only.

[timoria21]

And now the Nurbs curve AreOverlapped() method. In my opinion, the most complex part of the code for supporting Nurbs. Now I regret to have changed the original CC behavior with closed Nurbs curves. I am now handling them with one single PlineVertex only. You'll see in the relevant code that I will post. Probably, it would have been much easier to split closed Nurbs curves in two (midpoint) as CC does for arcs.

Basically, I project and check curves endpoints and endtangents. In one case, also curves length.

        /// <summary>
        /// Returns true if the curves are overlapping.
        /// </summary>
        /// <param name="cu1">First curve</param>
        /// <param name="cu2">Second curve</param>
        /// <param name="p0">First overlapping point</param>
        /// <param name="p1">Second overlapping point</param>
        /// <returns>True, if the two curves are overlapping.</returns>
        static bool AreOverlapped(GCurve cu1, GCurve cu2, out Point3D p0, out Point3D p1)
        {
            p0 = p1 = null;

            int succesfulProjCount = 0;

            const double TANTOL = 1e-3;
            const double COINCTOL = 1e-4;
            const double LENGTHTOL = 1e-4;

            if (AreEqual(cu1, cu2, COINCTOL))
            {
                p0 = cu1.StartPoint;
                p1 = cu1.EndPoint;

                return true;
            }

            if (Math.Abs(cu2.Length() - cu1.Length()) < LENGTHTOL) // precisely overlapped curves with same exact length
            {
                if (cu2.StartPoint.DistanceTo(cu1.StartPoint) < COINCTOL && // the same point
                    cu2.EndPoint.DistanceTo(cu1.EndPoint) < COINCTOL && // the same point
                    Vector3D.AreCoincident(cu2.StartTangent, cu1.StartTangent, TANTOL) &&
                    Vector3D.AreCoincident(cu2.StartTangent, cu1.StartTangent, TANTOL))
                {
                    p0 = cu1.StartPoint;
                    p1 = cu1.EndPoint;

                    return true;
                }

                if (cu2.StartPoint.DistanceTo(cu1.EndPoint) < COINCTOL && // the same point
                    cu2.EndPoint.DistanceTo(cu1.StartPoint) < COINCTOL && // the same point
                    Vector3D.AreOpposite(cu2.StartTangent, cu1.EndTangent, TANTOL) &&
                    Vector3D.AreOpposite(cu2.EndTangent, cu1.StartTangent, TANTOL))
                {
                    p0 = cu1.StartPoint;
                    p1 = cu1.EndPoint;
                    return true;
                }
            }

            Vector3D proCu1Start = null, proCu1End = null; // projection
            Vector3D tanCu1Start = null, tanCu1End = null; // tangent

            double t;

            // first curve start point on second curve
            if (cu2.StartPoint.DistanceTo(cu1.StartPoint) > COINCTOL && // not the same point
                cu2.EndPoint.DistanceTo(cu1.StartPoint) > COINCTOL &&   // not the same point

                cu2.Project(cu1.StartPoint, COINCTOL, false, out t))
            {
                Vector3D[] CK = cu2.Evaluate(t, 1);
                proCu1Start = CK[0];

                if (proCu1Start.AsPoint.DistanceTo(cu1.StartPoint) < COINCTOL) // TODO: mettere una tol relativa
                {
                    succesfulProjCount++;

                    tanCu1Start = CK[1];
                    tanCu1Start.Normalize();
                }
            }

            // first curve end point on second curve
            if (cu2.StartPoint.DistanceTo(cu1.EndPoint) > COINCTOL && // not the same point
                cu2.EndPoint.DistanceTo(cu1.EndPoint) > COINCTOL &&   // not the same point
                cu2.Project(cu1.EndPoint, COINCTOL, false, out t))
            {
                Vector3D[] CK = cu2.Evaluate(t, 1);
                proCu1End = CK[0];

                if (proCu1End.AsPoint.DistanceTo(cu1.EndPoint) < COINCTOL) // TODO: mettere una tol relativa
                {
                    succesfulProjCount++;

                    tanCu1End = CK[1];
                    tanCu1End.Normalize();
                }
            }

            Vector3D proCu2Start = null, proCu2End = null; // projection
            Vector3D tanCu2Start = null, tanCu2End = null; // tangent

            // second curve start point on first curve
            if (cu1.StartPoint.DistanceTo(cu2.StartPoint) > COINCTOL && // not the same point
                cu1.EndPoint.DistanceTo(cu2.StartPoint) > COINCTOL &&   // not the same point
                cu1.Project(cu2.StartPoint, COINCTOL, false, out t))
            {
                Vector3D[] CK = cu1.Evaluate(t, 1);
                proCu2Start = CK[0];

                if (proCu2Start.AsPoint.DistanceTo(cu2.StartPoint) < COINCTOL) // TODO: mettere una tol relativa
                {
                    succesfulProjCount++;

                    tanCu2Start = CK[1];
                    tanCu2Start.Normalize();
                }
            }

            if (cu1.StartPoint.DistanceTo(cu2.EndPoint) > COINCTOL && // not the same point
                cu1.EndPoint.DistanceTo(cu2.EndPoint) > COINCTOL &&   // not the same point
                cu1.Project(cu2.EndPoint, COINCTOL, false, out t))
            {
                Vector3D[] CK = cu1.Evaluate(t, 1);
                proCu2End = CK[0];

                if (proCu2End.AsPoint.DistanceTo(cu2.EndPoint) < COINCTOL) // TODO: mettere una tol relativa
                {
                    succesfulProjCount++;

                    tanCu2End = CK[1];
                    tanCu2End.Normalize();
                }
            }

            if (succesfulProjCount < 1)
            {
                return false;
            }

            if (proCu1End != null && proCu2Start != null && proCu1End != proCu2Start)
            {
                if (tanCu1End != null && Vector3D.AreParallel(tanCu1End, cu1.EndTangent, TANTOL) &&
                    tanCu2Start != null && Vector3D.AreParallel(tanCu2Start, cu2.StartTangent, TANTOL))
                {
                    p0 = proCu1End.AsPoint;
                    p1 = proCu2Start.AsPoint;

                    return true;
                }
            }

#region Start/End same curve 
            if (tanCu1Start != null && Vector3D.AreParallel(tanCu1Start, cu1.StartTangent, TANTOL) &&
                tanCu1End != null && Vector3D.AreParallel(tanCu1End, cu1.EndTangent, TANTOL))
            {
                p0 = cu1.StartPoint;
                p1 = cu1.EndPoint;

                return true;
            }

            if (tanCu2Start != null && Vector3D.AreParallel(tanCu2Start, cu2.StartTangent, TANTOL) &&
                tanCu2End != null && Vector3D.AreParallel(tanCu2End, cu2.EndTangent, TANTOL))
            {
                p0 = cu2.StartPoint;
                p1 = cu2.EndPoint;

                return true;
            }
#endregion

#region Start/End mixed curves
            if (tanCu2End != null && Vector3D.AreParallel(tanCu2End, cu2.EndTangent, TANTOL) &&
                tanCu1Start != null && Vector3D.AreParallel(tanCu1Start, cu1.StartTangent, TANTOL))
            {
                p0 = proCu2End.AsPoint;
                p1 = proCu1Start.AsPoint;

                return true;
            }

            if (tanCu2Start != null && Vector3D.AreParallel(tanCu2Start, cu2.StartTangent, TANTOL) &&
                tanCu1End != null && Vector3D.AreParallel(tanCu1End, cu1.EndTangent, TANTOL))
            {
                p0 = proCu2Start.AsPoint;
                p1 = proCu1End.AsPoint;

                return true;
            }
#endregion

#region Start/Start or End/End same curve
            if (tanCu2Start != null && Vector3D.AreParallel(tanCu2Start, cu2.StartTangent, TANTOL) &&
                tanCu1Start != null && Vector3D.AreParallel(tanCu1Start, cu1.StartTangent, TANTOL))
            {
                p0 = proCu2Start.AsPoint;
                p1 = proCu1Start.AsPoint;

                return true;
            }

            if (tanCu2End != null && Vector3D.AreParallel(tanCu2End, cu2.EndTangent, TANTOL) &&
                tanCu1End != null && Vector3D.AreParallel(tanCu1End, cu1.EndTangent, TANTOL))
            {
                p0 = proCu2End.AsPoint;
                p1 = proCu1End.AsPoint;

                return true;
            }
#endregion

#region The two curves have one coincident point at start/end
            if (proCu1Start != null)
            {
                if ((cu2.StartPoint.DistanceTo(cu1.EndPoint) < COINCTOL || cu2.EndPoint.DistanceTo(cu1.EndPoint) < COINCTOL) &&
                    Vector3D.AreParallel(cu1.EndTangent, cu2.EndTangent, TANTOL) &&
                    tanCu1Start != null && Vector3D.AreParallel(tanCu1Start, cu1.StartTangent, TANTOL))
                {
                    p0 = cu1.StartPoint;
                    p1 = cu1.EndPoint;

                    return true;
                }
            }

            if (proCu1End != null)
            {
                if ((cu2.StartPoint.DistanceTo(cu1.StartPoint) < COINCTOL || cu2.EndPoint.DistanceTo(cu1.StartPoint) < COINCTOL) &&
                    Vector3D.AreParallel(cu1.StartTangent, cu2.StartTangent, TANTOL) &&
                    tanCu1End != null && Vector3D.AreParallel(tanCu1End, cu1.EndTangent, TANTOL))
                {
                    p0 = cu1.StartPoint;
                    p1 = cu1.EndPoint;
                    return true;
                }
            }

            if (proCu2Start != null)
            {
                if ((cu1.StartPoint.DistanceTo(cu2.EndPoint) < COINCTOL || cu1.EndPoint.DistanceTo(cu2.EndPoint) < COINCTOL) &&
                    Vector3D.AreParallel(cu1.EndTangent, cu2.EndTangent, TANTOL) &&
                    tanCu2Start != null && Vector3D.AreParallel(tanCu2Start, cu2.StartTangent, TANTOL))
                {
                    p0 = cu2.StartPoint;
                    p1 = cu2.EndPoint;

                    return true;
                }
            }

            if (proCu2End != null)
            {
                if ((cu1.StartPoint.DistanceTo(cu2.StartPoint) < COINCTOL || cu1.EndPoint.DistanceTo(cu2.StartPoint) < COINCTOL) &&
                    Vector3D.AreParallel(cu1.StartTangent, cu2.StartTangent, TANTOL) &&
                    tanCu2End != null && Vector3D.AreParallel(tanCu2End, cu2.EndTangent, TANTOL))
                {
                    p0 = cu2.StartPoint;
                    p1 = cu2.EndPoint;

                    return true;
                }
            }
#endregion

            return false;
        }

Here the AreEqual() method:

        /// <summary>
        /// Returns true if the curves are identical (orientation can be opposite)
        /// </summary>
        private static bool AreEqual(Curve cu1, Curve cu2, double tol)
        {
            if (cu1.Degree != cu2.Degree ||
                cu1.KnotVector.Length != cu2.KnotVector.Length ||
                cu1.ControlPoints.Length != cu2.ControlPoints.Length)
            {
                return false;
            }

            int numU = cu1.ControlPoints.Length;

            for (int k = 0; k < numU; k++)
            {
                if (Point3D.DistanceSquared(cu1.ControlPoints[k].Euclid, cu2.ControlPoints[k].Euclid) > tol * tol &&
                    Point3D.DistanceSquared(cu1.ControlPoints[k].Euclid, cu2.ControlPoints[numU - k - 1].Euclid) > tol * tol)
                {
                    return false;
                }
            }

            // pick the smaller tolerance (curve domain * tolFac)
            double domainTol = cu1.Domain.Length < cu2.Domain.Length
                ? cu1.Domain.Length * Utility.TOL6
                : cu2.Domain.Length * Utility.TOL6;

            if (Math.Abs(cu1.Domain.Length - cu2.Domain.Length) > domainTol)
            {
                return false;
            }

            int len = cu1.KnotVector.Length;

            for (int k = 0; k < len; k++)
            {
                if (Math.Abs(cu1.KnotVector[k] - cu2.KnotVector[k]) > domainTol &&
                    Math.Abs(cu1.KnotVector[k] + cu2.KnotVector[len - k - 1]) > domainTol) // normally curve knot values are not reversed but negated (knot values always grow)
                {
                    return false;
                }
            }

            return true;
        }

[timoria21]

Here are some test cases that validate the methods above. Big black points are overlapping points. Small black points are Nurbs curves control points.

[timoria21]

addVertex() method:

        void addVertex(double x, double y, double bulge, Curve crv = null)
        {
            m_vertexes.Add(new PlineVertex(x, y, bulge) { Curve = crv });
        }

        void addVertex(PlineVertex vertex)
        {
            addVertex(vertex.X, vertex.Y, vertex.Bulge, vertex.Curve);
        }

I had to comment this check on all visitSegIndices() overloads:

        void visitSegIndices(Polyline pline, Vector point, ref int result)
        {
            /*if (m_vertexes.Count < 2) // closed Nurbs curve can have one vertex only 
            {
                return;
            }*/

            uint i;
            uint j;
            if (iSClosed)
            {
                i = 0;
                j = (uint)(m_vertexes.Count - 1);
            }
            else
            {
                i = 1;
                j = 0;
            }

            while (i < m_vertexes.Count && cavc.Polyline.visitor(pline, j, i, point, ref result))
            {
                j = i;
                i = i + 1;
            }
        }

visitorGetExtents() method now using SegmentType:

        static bool visitorGetExtents(uint i, uint j, Polyline pline, AABB result)
        {
            PlineVertex v1 = pline.m_vertexes[(int)i];
            PlineVertex v2 = pline.m_vertexes[(int)j];

            if (v1.SegmentType == 0)
            {
                if (v2.X < result.xMin)
                {
                    result.xMin = v2.X;
                }
                else if (v2.X > result.xMax)
                {
                    result.xMax = v2.X;
                }

                if (v2.Y < result.yMin)
                {
                    result.yMin = v2.Y;
                }
                else if (v2.Y > result.yMax)
                {
                    result.yMax = v2.Y;
                }

            }
            else if (v1.SegmentType == 1)
            {
                AABB arcBB = new AABB();
                //bounds of chord
                if (v1.X < v2.X)
                {
                    arcBB.xMin = v1.X;
                    arcBB.xMax = v2.X;
                }
                else
                {
                    arcBB.xMin = v2.X;
                    arcBB.xMax = v1.X;
                }

                if (v1.Y < v2.Y)
                {
                    arcBB.yMin = v1.Y;
                    arcBB.yMax = v2.Y;
                }
                else
                {
                    arcBB.yMin = v2.Y;
                    arcBB.yMax = v1.Y;
                }

                // it's an arc segment, find axes crossings (with origin at arc center) to determine full
                // extents
                ArcRadiusAndCenter arc = plineSegmentGlobalMembers.arcRadiusAndCenter(v1, v2);
                Vector circleXMaxPt = new Vector(arc.center.X + arc.radius, arc.center.Y);
                Vector circleYMaxPt = new Vector(arc.center.X, arc.center.Y + arc.radius);

                // must check if arc is an axis aligned half circle, in which case the isLeft checks will fail
                // at the boundaries so we must handle as special case
                if (utils.mathutils.fuzzyEqual(v1.X, arc.center.X))
                {
                    //y axis aligned half circle
                    bool updateXMin = (v1.Y > v2.Y) ^ (v1.bulgeIsNeg());
                    if (updateXMin)
                    {
                        arcBB.xMin = arc.center.X - arc.radius;
                    }
                    else
                    {
                        arcBB.xMax = circleXMaxPt.X;
                    }
                }
                else if (utils.mathutils.fuzzyEqual(v1.Y, arc.center.Y))
                {
                    //x axis aligned half circle
                    bool updateYMin = (v1.X > v2.X) ^ (v1.bulgeIsPos());
                    if (updateYMin)
                    {
                        arcBB.yMin = arc.center.Y - arc.radius;
                    }
                    else
                    {
                        arcBB.yMax = circleYMaxPt.Y;
                    }
                }
                else
                {
                    //not axis aligned, use isLeft checks to find quadrants and crossings
                    int startPtQuad = getQuadrant(v1.Pos, arc, circleXMaxPt, circleYMaxPt);
                    int endPtQuad = getQuadrant(v2.Pos, arc, circleXMaxPt, circleYMaxPt);
                    if (startPtQuad == 1)
                    {
                        if (endPtQuad == 2)
                        {
                            //crosses YMax
                            arcBB.yMax = circleYMaxPt.Y;
                        }
                        else if (endPtQuad == 3)
                        {
                            if (v1.bulgeIsNeg())
                            {
                                //crosses xMax then yMin
                                arcBB.xMax = circleXMaxPt.X;
                                arcBB.yMin = arc.center.Y - arc.radius;
                            }
                            else
                            {
                                //crosses yMax then xMin
                                arcBB.yMax = circleYMaxPt.Y;
                                arcBB.xMin = arc.center.X - arc.radius;
                            }
                        }
                        else if (endPtQuad == 4)
                        {
                            //crosses XMax
                            arcBB.xMax = circleXMaxPt.X;
                        }
                    }
                    else if (startPtQuad == 2)
                    {
                        if (endPtQuad == 1)
                        {
                            //crosses yMax
                            arcBB.yMax = circleYMaxPt.Y;
                        }
                        else if (endPtQuad == 3)
                        {
                            //crosses xMin
                            arcBB.xMin = arc.center.X - arc.radius;
                        }
                        else if (endPtQuad == 4)
                        {
                            if (v1.bulgeIsNeg())
                            {
                                //crosses yMax then xMax
                                arcBB.yMax = circleYMaxPt.Y;
                                arcBB.xMax = circleXMaxPt.X;
                            }
                            else
                            {
                                //crosses xMin then yMin
                                arcBB.xMin = arc.center.X - arc.radius;
                                arcBB.yMin = arc.center.Y - arc.radius;
                            }
                        }
                    }
                    else if (startPtQuad == 3)
                    {
                        if (endPtQuad == 1)
                        {
                            if (v1.bulgeIsNeg())
                            {
                                //crosses xMin then yMax
                                arcBB.xMin = arc.center.X - arc.radius;
                                arcBB.yMax = circleYMaxPt.Y;
                            }
                            else
                            {
                                //crosses yMin then xMax
                                arcBB.yMin = arc.center.Y - arc.radius;
                                arcBB.xMax = circleXMaxPt.X;
                            }
                        }
                        else if (endPtQuad == 2)
                        {
                            //crosses xMin
                            arcBB.xMin = arc.center.X - arc.radius;
                        }
                        else if (endPtQuad == 4)
                        {
                            //crosses yMin
                            arcBB.yMin = arc.center.Y - arc.radius;
                        }
                    }
                    else
                    {
                        Debug.Assert(startPtQuad == 4, "shouldn't get here without start pt in quadrant 4");

                        if (endPtQuad == 1)
                        {
                            //crosses xMax
                            arcBB.xMax = circleXMaxPt.X;
                        }
                        else if (endPtQuad == 2)
                        {
                            if (v1.bulgeIsNeg())
                            {
                                //crosses yMin then xMin
                                arcBB.yMin = arc.center.Y - arc.radius;
                                arcBB.xMin = arc.center.X - arc.radius;
                            }
                            else
                            {
                                //crossses xMax then yMax
                                arcBB.xMax = circleXMaxPt.X;
                                arcBB.yMax = circleYMaxPt.Y;
                            }
                        }
                        else if (endPtQuad == 3)
                        {
                            //crosses yMin
                            arcBB.yMin = arc.center.Y - arc.radius;
                        }
                    }
                }

                result.xMin = Math.Min(result.xMin, arcBB.xMin);
                result.yMin = Math.Min(result.yMin, arcBB.yMin);
                result.xMax = Math.Max(result.xMax, arcBB.xMax);
                result.yMax = Math.Max(result.yMax, arcBB.yMax);
            }
            else
            {
                v1.Curve.GetBBox(out Point3D min, out Point3D max); 

                result.xMin = min.X;
                result.yMin = min.Y;
                result.xMax = max.X;
                result.yMax = max.Y;
            }

            //return true to iterate all segments
            return true;
        }

visitor() method:

        static bool visitor(Polyline pline, uint i, uint j, Vector point, ref int windingNumber)
        {
            PlineVertex v1 = pline.m_vertexes[(int)i];
            PlineVertex v2 = pline.m_vertexes[(int)j];

            if (v1.SegmentType == 0)
            {
                lineVisitor(point, v1, v2, ref windingNumber);
            }
            else if (v1.SegmentType == 1)
            {
                arcVisitor(v1, v2, point, ref windingNumber);
            }
            else
            {
                curveVisitor(v1, v2, point, ref windingNumber);
            }

            return true;
        }

curveVisitor() method:

        static void curveVisitor(PlineVertex v1, PlineVertex v2, Vector point, ref int windingNumber)
        {
            double max = double.MinValue;

            foreach (Point4D hp in v1.Curve.ControlPoints)
            {
                Point3D euclid = hp.Euclid;

                if (euclid.X > max)
                {
                    max = euclid.X;
                }
            }

            if (point.X > max)
                return;

            Curve[] decomposedSegments = v1.BezSegs;
            Line ln = new Line(point.X, point.Y, max + (0.1 * Math.Abs(max)), point.Y);

            foreach (Curve curve in decomposedSegments)
            {
                Point3D[] res = curve.IntersectWith(ln);

                double tol = 1e-6;
                bool alignedWithEnd = false;
                bool alignedWithStart = false;

                Vector3D[] tgAtPoint = new Vector3D[res.Length];
                bool[] ptToSkip = new bool[res.Length];

                Point4D[] crvCtrlPts = curve.ControlPoints;
                int lastCtrlPtIx = crvCtrlPts.Length - 1;

                bool mustComputeBox = true;
                Point3D bMin = null;
                Point3D bMax = null;

                for (var i = 0; i < res.Length; i++)
                {
                    curve.Project(res[i], out double t);
                    if (Point3D.DistanceSquared(res[i], curve.PointAt(t)) > tol)
                    {
                        ptToSkip[i] = true;
                        continue;
                    }

                    InterPoint ip = (InterPoint)res[i];
                    Vector3D tg = curve.TangentAt(ip.u);
                    tgAtPoint[i] = tg;

                    alignedWithStart |= Math.Abs(ip.Y - crvCtrlPts[0].Euclid.Y) < tol;
                    alignedWithEnd |= Math.Abs(ip.Y - crvCtrlPts[lastCtrlPtIx].Euclid.Y) < tol;

                    // Check for tangency on closed curves: skip them
                    if (Vector3D.AngleBetween(tg, ln.StartTangent) < 0.001)
                    {
                            if (mustComputeBox)
                            {
                                mustComputeBox = false;
                                curve.GetBBox(out bMin, out bMax);
                            }

                            // If the curve is completely on one side of the intersecting line ignore it
                            bool completelyAbove = bMin.Y - ip.Y >= -Utility.TOL9;
                            bool completelyBelow = ip.Y - bMax.Y >= -Utility.TOL9;
                            if (completelyAbove || completelyBelow)
                                ptToSkip[i] = true;
                    }
                }

                if (alignedWithStart || alignedWithEnd)
                {
                    // Assuming the decomposed segment has a single curvature, we create a dummy arc with the same start and end points and call arcVisitor
                    curve.IsPlanar(1e-4, out Plane pCurve);
                    bool isCW = pCurve.AxisZ.Z < 0;

                    double tolConcidenceV1V2 = mathutils.realPrecision * mathutils.realPrecision;

                    bool startOfBezSegCoincideWithSlice = Point3D.DistanceSquared(new Point3D(v1.Pos.X, v1.Pos.Y), gCrvCtrlPts[0].Euclid) < tolConcidenceV1V2;
                    bool endOfBezSegCoincideWithSlice = Point3D.DistanceSquared(new Point3D(v2.Pos.X, v2.Pos.Y), gCrvCtrlPts[lastCtrlPtIx].Euclid) < tolConcidenceV1V2;

                    // If the segment's start or end coincide with original v1 and v2, use them to prevent collateral effects of arcVisitor due to rounding errors in splitting
                    Vector pos1 = startOfBezSegCoincideWithSlice ? v1.Pos : new Vector(gCrvCtrlPts[0].X, gCrvCtrlPts[0].Y);
                    Vector pos2 = endOfBezSegCoincideWithSlice ? v2.Pos : new Vector(gCrvCtrlPts[lastCtrlPtIx].X, gCrvCtrlPts[lastCtrlPtIx].Y);

                    double smallBulge = 0.2;
                    PlineVertex dummyV1 = new PlineVertex(pos1, isCW ? -smallBulge : smallBulge);
                    PlineVertex dummyV2 = new PlineVertex(pos2, 0);

                    arcVisitor(dummyV1, dummyV2, point, ref windingNumber);
                }
                else
                {
                    for (var i = 0; i < res.Length; i++)
                    {
                        if (ptToSkip[i])
                            continue;

                        Vector3D tg = tgAtPoint[i];

                        if (tg.Y > 0) // left and upward crossing
                        {
                            windingNumber += 1;
                        }
                        else
                        {
                            windingNumber -= 1;
                        }
                    }
                }
            }
        }