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

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

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

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


Возиться

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

Две пересекающиеся кубические кривые Безье

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

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


person zneak    schedule 28.10.2010    source источник
comment
Пожалуйста, не вычисляйте детерминанты, если вам это не нужно. Если вы хотите проверить сингулярность, вычислите наименьшее и наибольшее сингулярное значение. И если вам по какой-то причине нужен определитель, не используйте расширение Лапласа! Он имеет экспоненциальную временную сложность. Вы можете сделать это за O (n ^ 3)!   -  person sellibitze    schedule 28.10.2010
comment
Подключив матрицу Сильвестра к калькулятору матриц по адресу bluebit.gr/matrix-calculator, дало -2916 для определителя. Возможно, вам потребуется исправить вашу детерминантную функцию.   -  person Kyle Lutz    schedule 28.10.2010
comment
@Kyle Lutz: Да, я обнаружил это примерно через 5 минут после публикации и исправил свою детерминантную функцию.   -  person zneak    schedule 28.10.2010
comment
@sellibitze Я с радостью откажусь от него, как только кто-нибудь объяснит другой способ поиска пересечений кривых Безье.   -  person zneak    schedule 28.10.2010
comment
Как вы получили уравнения x = 9t ^ 3 - 9t ^ 2 - 3t + 2; y = 9t ^ 3 - 9t ^ 2 - 6t + 4 от контрольных точек кривых?   -  person user200783    schedule 28.10.2010
comment
@ Пол Бейкер. Я не знал, это сделал мой калькулятор. Он преобразовал известное нам кубическое уравнение Безье в (-P0 + 3*P1 - 3*P2 + P4) * t^3 + 3*(P0 - 2*P1 + P2) * t^2 - 3*(P0 - P1) * t + P0; и два уравнения, которые я показал, - это просто функции для X и Y двух вычитаемых кривых. WolframAlpha подтверждает, что они идентичны.   -  person zneak    schedule 28.10.2010
comment
Это идеальный алгоритм, но он обнаруживает, если две параметрические кривые пересекаются при одном и том же значении параметра (это то, на что указывает ответ @PaulBaker). Настоящая проблема (пересекаются ли кривые вообще?) - это двухвариативный полином, для которого вы хотите найти корни, проблема, для которой я не знаю, существует ли аналитическое решение. Я думаю, что вопрос следует отредактировать, чтобы включить это замечание.   -  person Ad N    schedule 19.09.2013


Ответы (8)


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

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

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

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

bezInt (B 1, B 2):

  1. Does bbox(B1) intersect bbox(B2)?
    • No: Return false.
    • Да: продолжить.
  2. Is area(bbox(B1)) + area(bbox(B2)) < threshold?
    • Yes: Return 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] Похоже, алгоритм разделения кривой Безье на две называется алгоритм де Кастельжау.

person j_random_hacker    schedule 28.10.2010
comment
Я подумывал об использовании этого решения. Это называется вырезкой по Безье. - person zneak; 28.10.2010
comment
@zneak: на самом деле, из статьи, на которую несколько раз указывалось в этой ветке (cagd .cs.byu.edu / ~ 557 / text / ch7.pdf), это метод подразделения Безье. - person Ad N; 19.09.2013

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

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

person tfinniga    schedule 03.12.2012
comment
Да, в настоящее время принятый ответ - вырезка Безье. Но здорово, что у вас есть аргументы, доказывающие, что это лучшее решение. - person zneak; 03.12.2012
comment
Я видел, что ... использование жирной линии вместо ограничивающих рамок даст гораздо более быстрое схождение, и это не намного больше работы. - person tfinniga; 03.12.2012

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

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

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

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

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

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

person user200783    schedule 28.10.2010

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

person Andrew    schedule 28.10.2010
comment
Я недостаточно хорошо разбираюсь в математике, чтобы знать, как решать параметрические уравнения. Не хотите ли привести пример с двумя кривыми из моего вопроса? - person zneak; 28.10.2010
comment
Я не знаю, как именно это сделать (какой метод использовать). Я просто знаю, что это такие методы. Возможно, вам удастся найти библиотеку численных методов C ++. - person Andrew; 29.10.2010

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

введите описание изображения здесь

foreach SampledLineSegment line1 in Bezier1
    foreach SampledLineSegment line2 in Bezier2
        if( line1 intersects line2 )
            then Bezier1 intersects Bezier2
person bobobobo    schedule 03.04.2013
comment
Вероятно, это разумное приближение, но мне больше нравится метод отсечения Безье, потому что он обеспечивает измеримую точность, и при этом он все еще довольно быстр. - person zneak; 03.04.2013
comment
Согласованный. Сначала я неправильно понял жалкую точность, поэтому был обеспокоен. - person bobobobo; 03.04.2013
comment
Существует ограниченный шанс пропустить пересечение, если две кривые имеют одну точку, касательную к двум кривым. Сегменты линии в конечном итоге окажутся немного ближе к среднему значению, поэтому кривые могут соприкасаться, а аппроксимация сегмента линии - нет. - person Tatarize; 10.07.2016
comment
Также существует ограниченная вероятность того, что у вас могут быть две разные кривые Безье, которые в основном параллельны, которые в конечном итоге пересекаются по ошибке. Но, с другой стороны, это было бы очень легко вычислить, и на самом деле вы могли бы разделить два пересекающихся сегмента на еще пять сегментов или еще много чего и получить лучшее представление о пересечении. Или просто подберите его, приблизив к минимальной ошибке, которую вы более или менее примете. - person Tatarize; 10.07.2016
comment
И, учитывая метод быстрого вычисления ограничивающей рамки для пересекающихся линий, позволяющий очень быстрое чистое сравнение не удается, вы могли бы лучше в основном решить это приближение за доли секунды. Потому что на самом деле сортировка кривой на монотонные сегменты и поиск одного или двух линейных сегментов, которые вам действительно нужно проверить, можно было бы уговорить в основном с постоянным временем. - person Tatarize; 10.07.2016

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

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

person GWW    schedule 28.10.2010
comment
Спасибо за ссылки. Хотя второй касается квадратичных кривых, которые, как я считаю, на порядок легче решить, поскольку я думаю, вам просто нужно вычесть их уравнения и использовать квадратную формулу. Я все же прочту вашу первую ссылку. - person zneak; 28.10.2010
comment
Ах, извините за вторую ссылку, надеюсь, первая поможет. - person GWW; 28.10.2010

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

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

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

person Tatarize    schedule 11.07.2016

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

Прежде всего, я хотел бы поблагодарить вас, 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);

В этот момент 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)

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

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

person Biły    schedule 05.11.2016
comment
Редкий случай, когда кривые принадлежат к одному семейству K, не так уж неожидан для программ для 2D-рисования. Пользователь может создать кривую Безье, а затем разделить ее на части. - person Malcolm McLean; 05.07.2017