?

Log in

No account? Create an account
nyaload

Журнал Пушыстого

Журнал Пушыстого

Previous Entry Share Next Entry
Геометрические Алгоритмы. Пересечение двух отрезков на плоскости.
nyaload
_winnie
Пусть есть два отрезка A,B и C,D. Надо понять, пересекаются ли они, и если да, то где именно. Во всех исходниках и книжках, которые я смотрел, это было сделано с кучей ненужного мусора, не смотря на то, что вроде бы очень простая и распространённая задача. Обычно как делают:
Параметризуют первый отрезок как A + (B-A)*u, u <= [0, 1], второй как C + (D-C)*v, v <= [0, 1]

Получается уравнение

A + (B-A)*u = C + (D-C)*v

Если переписать покомпонентно

Ax + (Bx-Ax)*u = Cx + (Dx-Cx)*v
Ay + (By-Ay)*u = Cy + (Dy-Cy)*v

То есть обычное уравнение из двух неизвестных.

Считают обычно по правилу крамера, потом смотрят, лежат ли корни в пределах от 0 до 1. Если определитель равен нулю, то считают что отрезки параллельны и не пересекаются.
Примечание: Или лежат на одной прямой, в этом случае лучше считать что не пересекаются. Так как любое дрожание начальных параметров приведёт к изменению решения, а дрожание такое есть всегда из-за неточности float, и часто бессмыслено выделять особый случай, который если и получится, то только если случайно float округлится до нужных значений. Если нам надо от алгоритма надо bool-значение "(Не)пересекаются", а не три энума "пересекаются в точке", "пересекаются по прямой", "не пересекаются", то проще считать второй случай сразу третьим.

Иногда сравнивают определитель сразу с 0, иногда пытаются ввести взятое с потолка EPSILON. Получается много if, взятое с потолка EPSILON (0.001 ? 0.0001 ?), непонятная работа если мы всё-таки начинаем делить большие числа на число чуть больше EPSILON, неправильная работа если не только числитель, но и знаменатель дроби меньше EPSILON.

Тут нужно не "разделять и властвовать", решая подзадачи (Ага, линейное уравнение. 1) решим 2) Ага, надо проверить что результат от 0 до 1), а решать задачу целиком.


Если у нас есть отрезок P + Q*t, то тогда нормаль к прямой проходящей через этот отрезок это perp(Q) = (-Qy, +Qx), а её уравнение это ((-Qy*x + Qx*y) - (-Qy*Px + Qx*Py) = 0). Переобозначив -Qy как A, и тп, получим привычное уравнение прямой (A*x + B*y + D == 0) Это же уравнение записанное как неравенство даст уравнение полуплоскости. Если подставить какую-то точку (x,y) в выражение A*x + B*y + С, то если оно получилось положительным - значит точка с одной стороны прямой, а если отрицательным - то с другой. Если ноль - значит, ровно на прямой. Степень отклонения от нуля этого выражения - это то, насколько точка далеко от прямой (это длина проекции отклонения на нормаль к прямой). Если оно почти ноль - значит, точка почти на прямой. Если выражение в два раза больше - значит точка в два раза дальше от прямой.



Мы не будем решать уравнение, а затем смотреть, лежат ли корни u и v там где нам надо, в промежутке 0-1. Определим сразу есть ли пересечение, сначала записав уравнение обоих отрезков, а затем посмотрев в каких полуплоскостях лежат концы отрезков. Если оба отрезка пересекают прямую своего отрезка-напарника (то есть концы лежат в разных полуплоскостях этой прямой), значит пересечение есть. Если концы какого-то отрезка лежат в одной полуплоскости - значит пересечения нет. Так как это выражение A*x + B*y + D характеризует насколько точка удалена от плоскости, то значит для двух концов отрезка два значения показывают на какие части делится отрезок пересечением. А это нам сразу даёт чему равно u или v в точке пересечения - u делит отрезок [0, 1] точно так же как точка пересечения делит отрезок AB, который задаётся как A + (B-A)*u.


Код:

    bool intersection(Point2f start1, Point2f end1, Point2f start2, Point2f end2, Point2f *out_intersection)
    {
        Point2f dir1 = end1 - start1;
        Point2f dir2 = end2 - start2;

        //считаем уравнения прямых проходящих через отрезки
        float a1 = -dir1.y;
        float b1 = +dir1.x;
        float d1 = -(a1*start1.x + b1*start1.y);

        float a2 = -dir2.y;
        float b2 = +dir2.x;
        float d2 = -(a2*start2.x + b2*start2.y);

        //подставляем концы отрезков, для выяснения в каких полуплоскотях они
        float seg1_line2_start = a2*start1.x + b2*start1.y + d2;
        float seg1_line2_end = a2*end1.x + b2*end1.y + d2;

        float seg2_line1_start = a1*start2.x + b1*start2.y + d1;
        float seg2_line1_end = a1*end2.x + b1*end2.y + d1;

        //если концы одного отрезка имеют один знак, значит он в одной полуплоскости и пересечения нет.
        if (seg1_line2_start * seg1_line2_end >= 0 || seg2_line1_start * seg2_line1_end >= 0) 
            return false;

        float u = seg1_line2_start / (seg1_line2_start - seg1_line2_end);
        *out_intersection =  start1 + u*dir1;

        return true;
    }
_Winnie C++ Colorizer



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


Смотри ещё: • Угол между двумя векторами
Смотри ещё: • Проверка точки внутри многоугольника

Симпатично.

Если подставить какую-то точку (x,y) в выражение A*x + B*y + С, то если оно получилось положительным - значит точка с одной стороны прямой, а если больше - то с другой.

напиши лучше отрицательным ^^ не каждый догадается что больше - значит меньше

а приятно с утра читать такие вот заметки. спасибо)

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

>а ты в какой точке нормаль проводишь?
В любой точке прямой :) в уравнении прямой после того как есть A и B - D можно получить подставив любую точку прямой, начало или конец отрезка, D при этом получится одинаковый.

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

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

Примечание: Или лежат на одной прямой, в этом случае лучше считать что не пересекаются.
Если я правильно понял, ты этот случай не выделяешь из соображений его неустойчивости. Но почему ты его относишь к случаю «не пересекаются»? Имхо, при малейших колебаниях начальных параметров изменится угол между отрезками. Тогда отрезки, возможно, будут пересекаться (быть может на конце).

Во всех исходниках и книжках, которые я смотрел, это было сделано с кучей ненужного мусора, не смотря на то, что вроде бы очень простая и распространённая задача... Считают обычно по правилу крамера, потом смотрят, лежат ли корни в пределах от 0 до 1.
Никогда с таким не сталкивался :-` Знаю такой способ: отрезки пересекаются тогда, и только тогда, когда одновременно:
1) Пересекаются их bounding box'ы.
2) Первый отрезок пересекает прямую, содержащую второй отрезок.
3) Второй отрезок пересекает прямую, содержащую первый отрезок.
Никаких делений не надо, нужно вычислять лишь векторные произведения, и сравнивать их с 0.

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

А для 3d?
Почему бы не поделить на ноль? В результате будет INF или NAN, а они не попадают в диапазон [0,1]. По операциям получается то же самое, только деление будет делаться всегда.

Хотя в 3d два отрезка скорее всего никогда не пересекаются.

А почему нет епсилонов? Надо бъ, тъ проверку делаеш на точку vs прямую, никак без епсилона...
Т.е. не понимаю, почему точка на прямой не будет 'дрожать'. Ето ведь не exact arithmetic, все сравнения у тебя зависят от размера отрезков, вот и float может подвести...

чем не годится [Graphics Gems III] pp. 199-202 "Faster Line Segment Intersection"?
в нем если не надо, то даже не умножают. правда, и точку пересечения не ищут.

Это что ли? http://tog.acm.org/GraphicsGems/gemsiii/insectc.c
Ну оно совершенно другое, куча if, почему-то fixed-point арифметика (автор пишет для мобильников?), почему-то больше кода. Как ты сказал, точку оно не ищет.
Пойнта в Bounding Box я не вижу - проще подставить x и y в Ax+By+D чем эта простыня из if.

Кстати, а нет ли под рукой ссылок на книги по аналитической геометрии вообще и приложениям её к комп.графике в частности?

К сожалению, концентрированного для компьютерной графики ничего не могу назвать. Разве что только исходники http://www.geometrictools.com//LibFoundation/Intersection/
Обычно у меня или в мозгу сразу выстраивается картинка, или я нахожу в гугле несколько решений, и выбираю из них наиболее "красивое" (т.е. более понятное и с меньшим количеством условных переходов).

<< Если определитель равен нулю, то считают что отрезки параллельны и не пересекаются.
Примечание: Или лежат на одной прямой, в этом случае лучше считать что не пересекаются. >>

Вообще-то, даже в случае, если отрезки лежат на одной прямой, то это еще не значит, что пересекаются:
либо не имеют общих точек, либо только одну, либо бесконечное множество(общий отрезок).
Следовательно, случай <<отрезки на одной прямой>> требует более тщательной проработки,или мы
ВЫНУЖДЕНЫ его отнести к варианту без пересечения.

Спасибо. Пригодилось =)

Здравствуйте!
Подскажите, как можно было бы решить такую задачу.
Есть отрезки x1-x2 и x3-x4, которые лежат на одной прямой. Нужно найти точки x5 и x6, которые являются концами отрезка, получившегося в результате наложения этих двух отрезков. Собственно, отрезка этого может и не быть, если наложения не происходит.

Предположил что x1 < x2 и x3 < x4, иначе надо сделать (x1', x2) = min(x1, x2), max(x1, x2) и тп.

(x5, x6) = (max(x1, x3), min(x2, x4)).

Если получилось что x5 > x6, пересечения нет.

http://segfault.kiev.ua/smart-questions-ru.html#homework

Извините за нескромный вопрос, но что делает оператор "-" класса Point2f? Хочу воспользоваться вашим алгоритмом, но нужно на Яву переделать...

эм. Если такие трудности с векторами - то что-либо с математикой будет очень тяжко делаться.
Point2f(x1, y1) - Point2f(x2, y2) это Point2f(x1 - x2, y1 - y2)
-Point2f(x1, y1) - это Point2f(-x1, -y1)


http://ru.wikipedia.org/wiki/Вектор_(геометрия)#.D0.9E.D0.BF.D0.B5.D1.80.D0.B0.D1.86.D0.B8.D0.B8_.D0.BD.D0.B0.D0.B4_.D0.B2.D0.B5.D0.BA.D1.82.D0.BE.D1.80.D0.B0.D0.BC.D0.B8
http://ru.wikipedia.org/wiki/Вектор_(алгебра)


Доброго Времени суток.
Решил для скорости не изобретать одноколесный велосипед, а взять Ваш код (переписав его на фортран).
Получилось вот как:

function Segm_intersection2D(start1, end1, start2, end2, out_intersection)
! Взято для скорости Winne http://users.livejournal.com/_winnie/152327.html
! ВАЖНО: Если Отрезки пересекаются по вершине, то возвращается .false.
! ВАЖНОЖ Если Отрезки не пересекаются, то возвращается (/HUGE(1.0d0), HUGE(1.0d0)/)
implicit none
logical :: Segm_intersection2D
real(8), intent(in) :: start1(2), end1(2), start2(2), end2(2)
real(8), intent(out) :: out_intersection(2)

real(8) :: dir1(2), dir2(2)
real(8) :: a1, b1, d1, a2, b2, d2, u, seg1_line2_start, seg1_line2_end, seg2_line1_start, seg2_line1_end

dir1 = end1 - start1;
dir2 = end2 - start2;

!считаем уравнения прямых проходящих через отрезки
a1 = -dir1(2)
b1 = +dir1(1)
d1 = -(a1*start1(1) + b1*start1(2))

a2 = -dir2(2)
b2 = +dir2(1)
d2 = -(a2*start2(1) + b2*start2(2))

!подставляем концы отрезков, для выяснения в каких полуплоскотях они
seg1_line2_start = a2*start1(1) + b2*start1(2) + d2;
seg1_line2_end = a2*end1(1) + b2*end1(2) + d2;

seg2_line1_start = a1*start2(1) + b1*start2(2) + d1;
seg2_line1_end = a1*end2(1) + b1*end2(2) + d1;

!если концы одного отрезка имеют один знак, значит он в одной полуплоскости и пересечения нет.
if ( (seg1_line2_start * seg1_line2_end >= 0) .or. (seg2_line1_start * seg2_line1_end >= 0)) then
Segm_intersection2D = .false.
out_intersection = (/HUGE(1.0d0), HUGE(1.0d0)/)
else

u = seg1_line2_start / (seg1_line2_start - seg1_line2_end);
out_intersection = start1 + u*dir1;
Segm_intersection2D = .true.
endif
endfunction Segm_intersection2D

Был закономерно, но приятно удивлен, что при следующих входных данных:
start1 = (/ 24050.0, 20100.0 /)
end1 = (/ 24050.0, 20050.0 /)
Start2 = (/ 0.0, 0.0/)
end2 = (/ 24050.0 20100.0/),
т.е. когда end2 равна start1 Ваша процедура выдает .false.

Другого и быть не могло, т.к. seg1_line2_start = 0 и seg2_line1_end = 0.
Подобные аспекты уже обсуждались выше в комментариях, но четко такого случая вроде бы не было.

Оно понятно: если поставить условие (seg1_line2_start * seg1_line2_end > 0) .or. (seg2_line1_start * seg2_line1_end > 0)
то нельзя различить счлучаю : оба отрезка лежат на прямой, но не пересекаются, оба отрезка лежат на одной прямой и перечсекаются(может быть по одной точке), оба отрезка не лежат на одной прямой, но они проходят через одну точку, объявленную на каждом из них, как концевая.


Алгоритм некорректен

Вы некорректно рассматриваете случай, когда отрезки на одной прямой. Действительно, если даны отрезки вида: (0, 0) - (5, 5) (4, 4) - (6, 6), то все определители равны нулю, и возвращается false, что неверно. marqueewinq@gmail.com (репост)

Edited at 2012-03-11 10:59 am (UTC)

Re: Алгоритм некорректен

Дело в том, что при дрожании входных параметров на любой эпсилон пересечение то будет, то не будет.

Для многих задач при таком дрожании - не принципиально, считать что пересечение есть или пересечения нет. Шум то ведь во входных значениях всегда есть, от неточности float и от погрешностей измерения. Данный алгоритм устойчив в том случае, если ответ устойчив к изменению параметров. Например, если нам надо решить задачу "пересекает ли игрок стену в компьютерной игре".

Действительно есть задачи где нужно больше чем true/false, знать так же "а может они соприкасаются", но там нужно очень хорошо понимать что такое соприкасаются. На расстоянии двух микрон - это уже соприкасаются, или ещё нет? Такой ответ кусок кода дать не может, это может решить только инженер который знает, микроны у него или нанометры.


Хорошая вещь) Я тут всю ночь на работе долбался как правильно пересечения искать... а тут такое... только один раз деление использовано))) да и то только чтобы найти точку) Автору респект. Помогло.

Пожалуйста :)