Bezier curves. A little about intersections and as simple as possible

Have you ever encountered building a (continuous) path to traverse a curve in a plane defined by line segments and Bezier curves?

It seems to be not a very difficult task: to join the segments of the curves in one path and go around it "without lifting the pen." A closed curve traverses in one direction, branches traverse forward and backward, and start and end at the same node.

Everything was fine until monstrous paths began to crawl out from under the hands of the designers, where individual curves could intersect or not exactly dock. The explanation was extremely simple - visually they all lie as they should, and for the machine that will bypass this path, such deviations are invisible.

Armed with the knowledge of the maximum permissible deviation, I started the research, the results of which I want to share.

The first thing I did was figure out how things are today (October 2020) with finding the intersection points of curves. Either I was looking in the wrong place, or maybe I was asking, but I couldn't find a simple solution. Although, the idea with the resultant of a pair of polynomials is quite interesting. Many different algorithms related to Bezier curves are collected here .

What I didn’t like in the known methods and what I don’t want to do is numerically search for the roots of polynomials, or even solve quadratic equations. I really don't want to explore curves for extremes. Anyway, I would like to avoid division, exponentiation, and everything that can lead to undefined behavior.

So what do you have to work with.

• Points are specified by type Point, like this:

  
  
  

  using Point = std::array<double, 2>;

  
  
  

  For Pointthe operators of addition, subtraction, multiplication by a scalar, scalar multiplication are defined.

  
  
  

• The value of the Radmissible deviation of points is set.

  
  
  

• Curves are defined by arrays of control (control) points std::vector<Point>.

  
  
  

• Almost coincident curves should be marked and, if possible, deleted, for example, if it is a forgotten duplicate (copy-paste is evil).

  
  
  

, , . ( ):

Point point(const std::vector<Point> &curve, double t, int n, int i)
{
    return n == 0 ? curve[i] : (point(curve, t, n - 1, i - 1) * (1 - t) + point(curve, t, n - 1, i) * t);
}

— , :

Point point(const std::vector<Point> &curve, double t)
{
    return point(curve, t, curve.size() - 1, curve.size() - 1);
}

, curve — : , .

— - , R:

template <>
struct std::less<Point>
{
    bool operator()(const Point &a, const Point &b, const double edge = R) const
    {
        for (auto i = a.size(); i-- > 0;) {
            if (a[i] + edge < b[i])
                return true;
            if (a[i] > b[i] + edge)
                return true;
        }
        return false;
    }
};

. , , , . — :

struct Rect
{
    Point topLeft, bottomRight;
    Rect(const Point &point);
    Rect(const std::vector<Point> &curve);
    bool isCross(const Rect &rect, const double edge) const
    {
        for (auto i = topLeft.size(); i-- > 0;) {
            if (topLeft[i] > rect.bottomRight[i] + edge
                || bottomRight[i] + edge < rect.topLeft[i])
                return false;
        }
        return true;
    }
};

void find(const std::vector<Point> &curveA, const std::vector<Point> &curveB,
          double tA, double tB, double dt)
{

    if (m_isSimilarCurve)
        return;

    Rect aRect(curveA);
    Rect bRect(curveB);
    if (!aRect.isCross(bRect, R))
        return;

    if (isNear(aRect.tl, aRect.br, R / 2) && isNear(bRect.tl, bRect.br, R / 2)) {
        // 3.1        
        addBest(curveA.front(), curveA.back(), curveB.front(), curveB.back(), tA, tB, dt);
        m_isSimilarCurve = (m_result.size() > curveA.size() * curveB.size());
        return;
    }

    const auto curveALeft = subCurve(curveA, 0, 0.5);
    const auto curveARight = subCurve(curveA, 0.5, 1.0);
    const auto curveBLeft = subCurve(curveB, 0, 0.5);
    const auto curveBRight = subCurve(curveB, 0.5, 1.0);

    const auto dtHalf = dt / 2;
    find(curveALeft, curveBLeft, tA, tB, dtHalf);
    find(curveALeft, curveBRight, tA, tB + dtHalf, dtHalf);
    find(curveARight, curveBLeft, tA + dtHalf, tB, dtHalf);
    find(curveARight, curveBRight, tA + dtHalf, tB + dtHalf, dtHalf);

}

- : , t t1 t2?

t = (t2 - t1) t' + t1 . t = t1, t = t2. , ( ) . : k t2 t1, k:

Point point(const std::vector<Point> &curve, double t1, int n, int i, double t2, int k)
{
    return n > k ? (point(curve, t1, n - 1, i - 1, t2, k) * (1 - t2) +
                    point(curve, t1, n - 1, i, t2, k) * t2)
                 : point(curve, t1, n, i);
}

, :

std::vector<Point> subCurve(const std::vector<Point> &curve, double t1, double t2)
{
    std::vector<Point> result(curve.size());
    for (auto k = result.size(); k-- > 0;) {
        result[k] = point(curve, t1, curve.size() - 1, curve.size() - 1, t2, curve.size() - 1 - k);
    }
    return result;
}

, , .

.

1. t1and t2can be any:
   
   
   • subCurve(curve, 1, 0)will give a curve that "moves" from the end point curveto the start point ,
   • subCurve(curve, 1, 2)extrapolates curvebeyond the last anchor point.
2. Some method implementations are intentionally omitted because they do not contain anything particularly interesting.
3. The function is point(curve, t)not suitable for calculating many points on a curve (for example, for rasterization); for this, calculation using a triangular matrix is ​​better suited.