Вы сталкивались когда-нибудь с построением (непрерывного) пути обхода кривой на плоскости, заданной отрезками и кривыми Безье?

Вроде бы не сильно сложная задача: состыковать отрезки кривых в один путь и обойти его "не отрывая пера". Замкнутая кривая обходится в одном направлении, ответвления — в прямом и обратном, начало и конец в одном узле.

Всё было хорошо, пока из-под рук дизайнеров не стали вылезать монструозные пути, где отдельные кривые могли пересекаться или не точно состыковываться. Объяснение было предельно простым — визуально они все лежат как надо, а для станка, который этот путь будет обходить, такие отклонения незаметны.

Вооружившись знанием о величине максимально допустимого отклонения, я приступил к исследованию, результатами которого хочу поделиться.

Первое что я сделал — это разобрался как на сегодняшний день (октябрь 2020) обстоят дела с поиском точек пересечения кривых. То ли я не там искал, то ли не то спрашивал, но найти простого решения не получилось. Хотя, идея с результантом пары полиномов довольно занимательна. Много разных алгоритмов связанных с кривыми Безье собрано здесь.

Что не понравилось в известных способах и что точно не хочется делать, так это численно искать корни полиномов, или даже решать квадратичные уравнения. Очень не хочется исследовать кривые на экстремумы. Да и вообще, хотелось бы избежать деления, возведения в степень и всего того, что может привести к неопределённому поведению.

Итак, с чем придётся работать.

• Точки задаются типом Point, например так:

  
  

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

  
  

  Для Point определены операторы сложения, вычитания, умножения на скаляр, скалярного умножения.

  
  

• Задана величина R допустимого отклонения точек.

  
  

• Кривые заданы массивами опорных (контрольных) точек std::vector<Point>.

  
  

• Почти совпадающие кривые следует отмечать и, по возможности, удалять, например, если это забытый дубликат (копипаста — это зло).

  
  

Первое, что точно понадобиться, это простой способ вычисления значения параметрической кривой. (Простой в смысле реализации и читабельности):

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 false;
        }
        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. t1 и t2 могут быть любыми:
   
   • subCurve(curve, 1, 0) даст кривую, которая "движется" от конечной точки curve к начальной,
   • subCurve(curve, 1, 2) экстраполирует curve за пределами последней опорной точки.
2. Реализации некоторых методов опущены намеренно, так как не содержат ничего особенно интересного.
3. Функция point(curve, t) не подходит для вычисления множества точек на кривой (например для растрезации), для этого лучше подойдёт вычисление с помощью треугольной матрицы.