Подтвердить что ты не робот

Проверка пересечения двух кубических кривых Безье

Для личного проекта мне нужно выяснить, пересекаются ли две кривые Безье. Мне не нужно знать, где: мне просто нужно знать, делают ли они это. Тем не менее, мне нужно будет сделать это быстро.

Я убираю место, и я нашел несколько ресурсов. В основном, этот вопрос здесь, который дал многообещающий ответ.

Итак, после того, как я понял, что такое матрица Сильвестра, что является determinant, что означает resultant и почему это полезно, я думал, что понял, как работает решение. Однако реальность заставляет различать, и она не работает так хорошо.


Мессинг вокруг

Я использовал свой графический калькулятор, чтобы нарисовать два сплайна Безье (которые мы будем называть B 0 и B 1), которые пересекаются. Их координаты следуют (P 0, P 1, P 2, P 3):

(1, 1) (2, 4) (3, 4) (4, 3)
(3, 5) (3, 6) (0, 1) (3, 1)

Результат следующий: B 0 является "горизонтальной" кривой и B 1 другой:

Two cubic Bézier curves that intersect

В соответствии с вышеизложенным ответом на верхний голос я вычитал B 0 в B 1. Это оставило меня с двумя уравнениями (оси X и Y), которые, согласно моему калькулятору, следующие:

x = 9t^3 - 9t^2 - 3t + 2
y = 9t^3 - 9t^2 - 6t + 4

Матрица Сильвестра

И из этого я построил следующую матрицу Сильвестра:

9  -9  -3   2   0   0
0   9  -9  -3   2   0
0   0   9  -9  -3   2
9  -9  -6   4   0   0
0   9  -9  -6   4   0
0   0   9  -9  -6   4

После этого я создал функцию С++ для вычисления детерминант матриц с помощью Расширение Лапласа:

template<int size>
float determinant(float* matrix)
{
    float total = 0;
    float sign = 1;
    float temporaryMatrix[(size - 1) * (size - 1)];
    for (int i = 0; i < size; i++)
    {
        if (matrix[i] != 0)
        {
            for (int j = 1; j < size; j++)
            {
                float* targetOffset = temporaryMatrix + (j - 1) * (size - 1);
                float* sourceOffset = matrix + j * size;
                int firstCopySize = i * sizeof *matrix;
                int secondCopySize = (size - i - 1) * sizeof *matrix;
                memcpy(targetOffset, sourceOffset, firstCopySize);
                memcpy(targetOffset + i, sourceOffset + i + 1, secondCopySize);
            }
            float subdeterminant = determinant<size - 1>(temporaryMatrix);
            total += matrix[i] * subdeterminant * sign;
        }
        sign *= -1;
    }
    return total;
}

template<>
float determinant<1>(float* matrix)
{
    return matrix[0];
}

Кажется, он работает довольно неплохо на относительно маленьких матрицах (2x2, 3x3 и 4x4), поэтому я ожидаю, что он будет работать и с матрицами 6x6. Однако я не проводил обширных тестов, и есть вероятность, что он сломался.


Проблема

Если я правильно понял ответ от другого вопроса, определитель должен быть 0, так как кривые пересекаются. Однако, подавая мою программу матрицу Сильвестра, которую я сделал выше, она -2916.

Это ошибка на моем конце или на их конце? Каков правильный способ выяснить, пересекаются ли две кубические кривые Безье?

4b9b3361

Ответ 1

Пересечение кривых Безье выполняется (очень круто) Asymptote векторный графический язык: искать intersect() здесь.

Хотя они не объясняют алгоритм, который они фактически используют там, кроме как сказать, что он из p. 137 из "Книги метафонов", похоже, что ключ к ней - два важных свойства кривых Безье (которые объясняются в другом месте на этом сайте, хотя я не могу найти страницу прямо сейчас):

  • Кривая Безье всегда содержится в ограничивающей рамке, определяемой четырьмя контрольными точками
  • Кривая Безье всегда может быть подразделена при произвольном значении t на две кривые безрезерца

С этими двумя свойствами и алгоритмом пересечения многоугольников вы можете выполнить произвольную точность:

bezInt (B 1, B 2):

  1. Либо bbox (B 1) пересекает bbox (B 2)?
    • Нет: возвращает false.
    • Да: Продолжить.
  2. Область (bbox (B 1)) + область (bbox (B 2)) < порог?
    • Да: Вернуть true.
    • Нет: Продолжить.
  3. Разделить B 1 на B 1a и B 1b при t = 0.5
  4. Разделить B 2 на B 2a и B 2b при t = 0.5
  5. Возврат bezInt (B 1a, B 2a) || bezInt (B 1a, B 2b) || bezInt (B 1b, B 2a) || bezInt (B 1b, B 2b).

Это будет быстро, если кривые не пересекаются, - это то, что обычный случай?

[EDIT] Похоже, что алгоритм разделения кривой Безье в двух называется алгоритмом де Кастеляу.

Ответ 2

Если вы делаете это для производственного кода, я бы предложил алгоритм отсечения Безье. Он хорошо объяснил в разделе 7.7 этого бесплатного онлайн-текста CAGD (pdf), работает для любой степени кривой Безье и является быстрым и надежным.

Хотя использование стандартных корневых кодов или матриц может быть более простым с математической точки зрения, обрезание Безье относительно легко реализовать и отладить, и на самом деле будет иметь меньше ошибок с плавающей запятой. Это происходит потому, что всякий раз, когда он создает новые числа, он делает взвешенные средние (выпуклые комбинации), поэтому нет возможности экстраполяции на основе шумных входов.

Ответ 3

Это ошибка на моем конце или на их конце?

Вы основываете свою интерпретацию определителя на 4-м комментарии к этому ответу? Если это так, я считаю, что там, где ошибка. Воспроизведение комментария здесь:

Если детерминант равен нулю, существует корень в X и Y в * точно такой же значение t, поэтому существует пересечение двух кривых. (t может не находиться в интервале 0..1 хотя).

Я не вижу никаких проблем с этой частью, но автор продолжает:

Если определитель равен < > 0, вы можете убедитесь, что кривые не пересекаются где угодно.

Я не думаю, что это правильно. Вполне возможно, что две кривые пересекаются в месте, где значения t отличаются, и в этом случае будет пересечение, даже если матрица имеет ненулевой детерминант. Я считаю, что это то, что происходит в вашем случае.

Ответ 4

Я никоим образом не специалист по подобным вещам, но я следую за хорошим блоком в котором много говорится о кривых. У него есть ссылка на две хорошие статьи, рассказывающие о вашей проблеме (вторая ссылка имеет интерактивную демонстрацию и некоторый исходный код). Другие люди могут иметь гораздо лучшее представление о проблеме, но я надеюсь, что это поможет!

http://cagd.cs.byu.edu/~557/text/ch7.pdf

Ответ 5

Я не знаю, как быстро это будет, но если у вас есть две кривые C1 (t) и C2 (k), они пересекаются, если C1 (t) == C2 (k). Таким образом, для двух переменных (t, k) у вас есть два уравнения (на x и на y). Вы можете решить систему с помощью численных методов с достаточной для вас точностью. Когда вы найдете t, k параметров, вы должны проверить, есть ли параметр на [0, 1]. Если они пересекаются на [0, 1].

Ответ 6

Это сложная проблема. Я бы разделил каждую из двух кривых Безье на 5-10 сегментов дискретной линии и просто сделал пересечения линий.

enter image description here

foreach SampledLineSegment line1 in Bezier1
    foreach SampledLineSegment line2 in Bezier2
        if( line1 intersects line2 )
            then Bezier1 intersects Bezier2

Ответ 7

Я бы сказал, что самый простой и, скорее всего, самый быстрый ответ - разбить их на очень мелкие линии и найти точки, где кривые пересекаются, если они действительно делают.

public static void towardsCubic(double[] xy, double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3, double t) {
    double x, y;
    x = (1 - t) * (1 - t) * (1 - t) * x0 + 3 * (1 - t) * (1 - t) * t * x1 + 3 * (1 - t) * t * t * x2 + t * t * t * x3;
    y = (1 - t) * (1 - t) * (1 - t) * y0 + 3 * (1 - t) * (1 - t) * t * y1 + 3 * (1 - t) * t * t * y2 + t * t * t * y3;
    xy[0] = x;
    xy[1] = y;
}

public static void towardsQuad(double[] xy, double x0, double y0, double x1, double y1, double x2, double y2, double t) {
    double x, y;
    x = (1 - t) * (1 - t) * x0 + 2 * (1 - t) * t * x1 + t * t * x2;
    y = (1 - t) * (1 - t) * y0 + 2 * (1 - t) * t * y1 + t * t * y2;
    xy[0] = x;
    xy[1] = y;
}

Очевидно, что ответ грубой силы плохо. Ответьте на bo {4}, и там много линейной геометрии и обнаружения столкновений, которые действительно помогут довольно много.


Выберите количество разрезов для кривых. Что-то вроде 100 должно быть здорово.

Мы берем все сегменты и сортируем их на основе наибольшего значения y, которое они имеют. Затем добавим второй флаг в список для наименьшего значения y для этого сегмента.

Мы сохраняем список активных ребер.

Мы перебираем через отсортированный по y список сегментов, когда мы сталкиваемся с ведущим сегментом, мы добавляем его в активный список. Когда мы нажимаем значение отмеченного значения small-y, мы удаляем этот сегмент из активного списка.

Затем мы можем просто перебирать весь набор сегментов, что соответствует строке сканирования, увеличивая нашу у монотонно, поскольку мы просто перебираем список. Мы перебираем значения в нашем отсортированном списке, который обычно удаляет один сегмент и добавляет новый сегмент (или для разделенных и объединенных узлов, добавляет два сегмента или удаляет два сегмента). Таким образом, сохраняя активный список релевантных сегментов.

Мы запускаем проверку пересечения с быстрым сбоем, добавляя новый активный сегмент в список активных сегментов, против только этого сегмента и активных сегментов.

Таким образом, мы всегда знаем точно, какие сегменты линии имеют значение, поскольку мы итерируем по отборным сегментам наших кривых. Мы знаем, что эти сегменты разделяют совпадение в y-координатах. Мы можем быстро вывести из строя любой новый сегмент, который не перекрывается относительно его х-коордов. В редком случае, когда они пересекаются в х-координатах, мы затем проверяем, пересекаются ли эти сегменты.

Это, скорее всего, уменьшит количество проверок пересечения линий в основном на количество пересечений.

foreach(segment in sortedSegmentList) {
    if (segment.isLeading()) {
        checkAgainstActives(segment);
        actives.add(segment);
    }
    else actives.remove(segment)
}

checkAgainstActive() будет просто проверять, перекрывает ли этот сегмент и любой сегмент активного списка их х-координаты, а затем выполняет проверку пересечения строк и предпринимает соответствующие действия.


Также обратите внимание, что это будет работать для любого числа кривых, любого типа кривых, любого порядка кривых в любой смеси. И если мы перебираем весь список сегментов, он найдет каждое пересечение. Он найдет каждую точку, в которой Безье пересекает круг или любое пересечение, которое дюжина квадратичных кривых Безье имеет друг с другом (или сами) и все в одну и ту же долю секунды.

Часто цитируемый документ главы 7 для алгоритма разбиения:

"Как только пара кривых была разделена настолько, что каждая из них может быть аппроксимирована сегментом линии с точностью до допустимого значения

Мы можем буквально просто пропустить посредника. Мы можем сделать это достаточно быстро, чтобы просто сравнить сегменты линии с допустимым количеством ошибок от кривой. В конце концов, это типичный ответ.


Во-вторых, обратите внимание, что основная часть увеличения скорости для обнаружения столкновения здесь, а именно упорядоченный список сегментов, отсортированный по их наивысшему y для добавления в активный список, и самый низкий y для удаления из активного списка, может Аналогично делается и для полигонов корпуса кривой Безье. Наш сегмент линии - это просто многоугольник порядка 2, но мы можем делать любое количество точек так же тривиально, и проверять ограничивающий прямоугольник всех контрольных точек в любом порядке кривой, с которой мы имеем дело. Поэтому вместо списка активных сегментов у нас будет список активных точек полигонов корпуса. Мы могли бы просто использовать алгоритм Де Кастеляу, чтобы разбить кривые, таким образом, выборку, как разделенные кривые, а не сегменты линии. Таким образом, вместо 2 пунктов у нас будет 4 или 7 или еще что-то еще, и запустите ту же процедуру, которая будет весьма склонна к быстрому сбою.

Добавление соответствующей группы точек в max y, удаление ее при min y и сравнение только активного списка. Таким образом, мы можем быстро и лучше реализовать алгоритм подразделения Безье. Просто найдя перекрывающееся прямоугольное перекрытие, а затем разделив наложенные кривые и удалив те, которые этого не сделали. Поскольку, опять же, мы можем делать любое количество кривых, даже те, которые подразделяются на кривые предыдущей итерации. И, подобно нашему приближению сегмента, быстро решаются для каждой позиции пересечения между сотнями различных кривых (даже разных порядков) очень быстро. Просто проверьте все кривые, чтобы увидеть, перекрываются ли ограничивающие поля, если они это сделают, мы разделим их, пока наши кривые не станут достаточно маленькими или не выйдут из них.

Ответ 8

Да, я знаю, эта тема принята и закрыта надолго, но...

Во-первых, я хотел бы поблагодарить вас, zneak, за вдохновение. Ваши усилия позволяют найти правильный путь.

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

Отсутствующий трюк, который позволяет решить проблему способом, предложенным zneak: достаточно укоротить одну из кривых с коэффициентом k 0, тогда детерминант матрицы Сильвестра будет равен нулю. Очевидно, что если укороченная кривая будет пересекаться, тогда будет и оригинал. Теперь задача изменяется на поиск подходящего значения для коэффициента k. Это привело к проблеме решения одномерного многочлена 9 степени. Реальные и положительные корни этого многочлена отвечают за потенциальные точки пересечения. (Это не должно быть неожиданностью, две кубические кривые Безье могут пересекаться до 9 раз). Окончательный выбор выполняется, чтобы найти только те k факторы, которые обеспечивают 0 <= t <= 1 для обеих кривых.

Теперь код Maxima, где начальная точка установлена ​​в 8 точек, предоставлена ​​ zneak:

p0x:1; p0y:1;
p1x:2; p1y:4;
p2x:3; p2y:4;
p3x:4; p3y:3;

q0x:3; q0y:5;
q1x:3; q1y:6;
q2x:0; q2y:1;
q3x:3; q3y:1;

c0x:p0x;
c0y:p0y;
c1x:3*(p1x-p0x);
c1y:3*(p1y-p0y);
c2x:3*(p2x+p0x)-6*p1x;
c2y:3*(p2y+p0y)-6*p1y;
c3x:3*(p1x-p2x)+p3x-p0x;
c3y:3*(p1y-p2y)+p3y-p0y;

d0x:q0x;
d0y:q0y;
d1x:3*(q1x-q0x);
d1y:3*(q1y-q0y);
d2x:3*(q2x+q0x)-6*q1x;
d2y:3*(q2y+q0y)-6*q1y;
d3x:3*(q1x-q2x)+q3x-q0x;
d3y:3*(q1y-q2y)+q3y-q0y;

x:c0x-d0x + (c1x-d1x*k)*t+ (c2x-d2x*k^2)*t^2+ (c3x-d3x*k^3)*t^3;
y:c0y-d0y + (c1y-d1y*k)*t+ (c2y-d2y*k^2)*t^2+ (c3y-d3y*k^3)*t^3;

z:resultant(x,y,t);

На этом этапе Максима ответил:

(%o35)−2*(1004*k^9−5049*k^8+5940*k^7−1689*k^6+10584*k^5−8235*k^4−2307*k^3+1026*k^2+108*k+76)

Пусть Maxima разрешает это уравнение:

rr: float( realroots(z,1e-20))  

Ответ:

(%o36) [k=−0.40256438624399,k=0.43261490325108,k=0.84718739982868,k=2.643321910825066,k=2.71772491293651]

Теперь код для выбора правильного значения k:

for item in rr do ( 
  evk:ev(k,item),
  if evk>0 then (  
    /*print("k=",evk),*/ 
    xx:ev(x,item),  rx:float( realroots(xx,1e-20)),/*print("x(t)=",xx,"   roots: ",rx),*/
    yy:ev(y,item),  ry:float( realroots(yy,1e-20)),/*print("y(t)=",yy,"   roots: ",ry),*/
    for it1 in rx do (  t1:ev(t,it1),
    for it2 in ry do (  t2:ev(t,it2),
       dt:abs(t1-t2),
       if dt<1e-10 then (
         /*print("Common root=",t1,"  delta t=",dt),*/
         if (t1>0) and (t1<=1) then ( t2:t1*evk,
         if (t2>0) and (t2<=1) then (                           
                 x1:c0x + c1x*t1+ c2x*t1^2+ c3x*t1^3,
                 y1:c0y + c1y*t1+ c2y*t1^2+ c3y*t1^3,
                 print("Intersection point:     x=",x1, "      y=",y1)
)))))/*,disp ("-----")*/
));

Максимальный ответ:

"Intersection point:     x="1.693201254437358"      y="2.62375005067273
(%o37) done

Тем не менее, это не только мед. Представленный метод непригоден, если:

  • P0 = Q0 или, более общо, если P0 лежит на второй кривой (или ее продолжении). Можно попытаться поменять кривые.
  • чрезвычайно редкие случаи, когда обе кривые принадлежат одному K-семейству (например, их бесконечные расширения одинаковы).
  • Следите за (sqr (c3x) + sqr (c3y)) = 0 или (sqr (d3x) + sqr (3y)) = 0 случаев, здесь квадратичные притворяются кубическими кривыми Безье.

Можно спросить, почему сокращение выполняется только один раз. Этого достаточно из-за обратного обратного закона, который был обнаружен сперва, но это еще одна история.