Проверка, пересекаются ли две кубические кривые Безье - PullRequest
30 голосов
/ 28 октября 2010

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

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

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


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

Я использовал свой графический калькулятор для рисования двух сплайнов Безье (которые мы назовем 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

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

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.

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

Ответы [ 8 ]

15 голосов
/ 28 октября 2010

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

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

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

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

bezInt (B 1 , B 2 ):

Пересекает ли bbox (B 1 ) bbox (B 2 )?
  • Нет: вернуть false.
  • Да: продолжить.
Is area (bbox (B 1 )) + area(bbox (B 2 )) <порог?<ul> Да: Вернуть true. Нет: Продолжить. Разделить B 1 на B 1a и B 1b при t = 0,5 Разделить B 2 на B 2a и B 2b при t = 0,5 Return bezInt (B 1a , B 2a ) ||bezInt (B 1a , B 2b ) ||bezInt (B 1b , B 2a ) ||bezInt (B 1b , B 2b ).

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

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

7 голосов
/ 03 декабря 2012

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

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

3 голосов
/ 28 октября 2010

Это ошибка с моей стороны или с их стороны?

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

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

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

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

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

2 голосов
/ 03 апреля 2013

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

enter image description here

foreach SampledLineSegment line1 in Bezier1
    foreach SampledLineSegment line2 in Bezier2
        if( line1 intersects line2 )
            then Bezier1 intersects Bezier2
2 голосов
/ 28 октября 2010

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

2 голосов
/ 28 октября 2010

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

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

0 голосов
/ 05 ноября 2016

Да, я знаю, эта тема уже давно принята и закрыта, но ...

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

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

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

Теперь код 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);

На данный момент Maxima ответила:

(%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)

Пусть Максима решит это уравнение:

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 случаев, здесь квадратик притворяется кубическим Безьекривые.

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

0 голосов
/ 11 июля 2016

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

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 список сегментов, когда мы сталкиваемся с ведущим сегментом, мы добавляем его в активный список.Когда мы нажимаем на значение, помеченное маленьким-y, мы удаляем этот сегмент из активного списка.

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

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

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

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

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

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


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

Часто цитируемое Глава 7документ отмечает, для алгоритма подразделения:

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

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


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

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

...