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