[C++, Алгоритмы, Математика] Кривые Безье. Немного о пересечениях и как можно проще

Автор Сообщение
news_bot ®

Стаж: 6 лет 9 месяцев
Сообщений: 27286

Создавать темы news_bot ® написал(а)
11-Окт-2020 19:30

Article
Вы сталкивались когда-нибудь с построением (непрерывного) пути обхода кривой на плоскости, заданной отрезками и кривыми Безье?
Вроде бы не сильно сложная задача: состыковать отрезки кривых в один путь и обойти его "не отрывая пера". Замкнутая кривая обходится в одном направлении, ответвления — в прямом и обратном, начало и конец в одном узле.
Всё было хорошо, пока из-под рук дизайнеров не стали вылезать монструозные пути, где отдельные кривые могли пересекаться или не точно состыковываться. Объяснение было предельно простым — визуально они все лежат как надо, а для станка, который этот путь будет обходить, такие отклонения незаметны.
Вооружившись знанием о величине максимально допустимого отклонения, я приступил к исследованию, результатами которого хочу поделиться.
Первое что я сделал — это разобрался как на сегодняшний день (октябрь 2020) обстоят дела с поиском точек пересечения кривых. То ли я не там искал, то ли не то спрашивал, но найти простого решения не получилось. Хотя, идея с результантом пары полиномов довольно занимательна. Много разных алгоритмов связанных с кривыми Безье собрано здесь.
Что не понравилось в известных способах и что точно не хочется делать, так это численно искать корни полиномов, или даже решать квадратичные уравнения. Очень не хочется исследовать кривые на экстремумы. Да и вообще, хотелось бы избежать деления, возведения в степень и всего того, что может привести к неопределённому поведению.

Пример

SPL
Если пытаться продолжить кривую, то не факт что она вообще пересечётся с другой кривой, хотя они находятся достаточно близко

Итак, с чем придётся работать.
  • Точки задаются типом 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)
{

1. Проверить, что эти кривые ещё не отмечены как подобные.

SPL
if (m_isSimilarCurve)
        return;

2. Проверить, что ограничивающие прямоугольники кривых пересекаются

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

3. Если отрезки кривых меньше R/2, то можно считать, что пересечение найдено

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

4. Разделить кривые

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

5. Продолжить поиск для каждого отрезка кривой

SPL
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) не подходит для вычисления множества точек на кривой (например для растрезации), для этого лучше подойдёт вычисление с помощью треугольной матрицы.

===========
Источник:
habr.com
===========

Похожие новости: Теги для поиска: #_c++, #_algoritmy (Алгоритмы), #_matematika (Математика), #_krivye_beze (кривые безье), #_polinomy_bernshtejna (полиномы бернштейна), #_rekursivnyj_algoritm (рекурсивный алгоритм), #_c++, #_c++, #_algoritmy (
Алгоритмы
)
, #_matematika (
Математика
)
Профиль  ЛС 
Показать сообщения:     

Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете голосовать в опросах
Вы не можете прикреплять файлы к сообщениям
Вы не можете скачивать файлы

Текущее время: 23-Ноя 04:40
Часовой пояс: UTC + 5