Пересечение отрезков в последнее время стало основным моим занятием. Зафиксирую, почему именно так, а то потом забудется. Как найти пересечения прямых рассмотрено в статье Пересечение прямых, угол и координаты пересечения. Нас интересуют именно отрезки. Решение векторно-алгебраическое, поэтому будет интересно.
Теория
Итак, у нас есть отрезок AB, заданный точками A и B. И отрезок CD, заданный точками C и D соответственно. Необходимо найти точку их пересечения, либо сообщить, что отрезки не пересекаются.
Векторная часть
Представим отрезок AB, как вектор, направленный из A в B. И отрезок CD, тоже как вектор, из C в D. Сразу предположим хороший случай, когда они пересекаются в некоторой точке.
Введём понятие коэффициентов t и u, которые означают часть векторов AB (t) и CD (u) до точки пересечения.
Очевидно, что значения коэффициентов t и u лежат в диапазоне [0..1]. Таким образом, получаем вектора AS и CS, которые сходятся в точке пересечения S.
|
1 2 |
AS = t*AB CS = u*CD |
Введём новый вектор AC, из A в C.
Замечаем, что он является разностью AS и CS:
|
1 |
AC = AS - CS = t*AB - u*CD |
Таким образом, имеем систему уравнений:
|
1 2 |
t*AB.X-u*CD.X = AC.X t*AB.Y-u*CD.Y = AC.Y |
Которую будем решать алгебраически, методом Крамера.
Алгебраическая часть
Имеем следующую систему:
Где:
|
1 2 3 4 5 6 |
A1 = AB.X A2 = AB.Y B1 = CD.X B2 = CD.Y C1 = AC.X C2 = AC.Y |
Определители будут такими:
Убираем знак минуса в числителе и знаменателе, получаем решение:
Если определитель Δab = 0, значит отрезки параллельны друг другу.
Если параметр t=0, значит точка пересечения совпадает с A. При t=1, точка пересечения совпадает с B. Признаком «попадания» в отрезок AB является нахождение параметра t ∈ [0, 1]. Аналогично для параметра u.
Практическая часть
Напишем функцию, которая бы определяла факт пересечения отрезков, и, в случае пересечения, возвращала параметры t и u.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 |
// Проверка пересечения 2D-отрезков и возврат параметров пересечения function IsSegmentsIntersect(const A, B, C, D: TPointF; out t, u, denom: Double; Epsilon: Double): Boolean; var Dt, Du: Double; AB, CD, AC: TPointF; begin // Инициализируем параметры, даже если результат будет False Result := False; t := 0; u := 0; // Находим векторы AB := B-A; // Вектор AB CD := D-C; // Вектор CD AC := C-A; // Вектор AC // Находим определитель Δab denom := AB.X*CD.Y - AB.Y*CD.X; // Проверяем if IsZero(denom, Epsilon) then exit; // Параллельные сегменты, выходим // Находим определители Δt и Δu Dt := AC.X*CD.Y - AC.Y*CD.X; Du := AC.X*AB.Y - AC.Y*AB.X; // Решаем систему уравнений методом Крамера t := Dt/denom; u := Du/denom; // Проверка на попадание в отрезки Result := (t >= 0) and (t <= 1) and (u >= 0) and (u <= 1); end; |
Пересечение отрезков в 2D
Получив факт пересечения, нам достаточно одного параметра t, чтобы вычислить точку пересечения. Если отрезки пересекаются, то использование u для отрезка CD даст нам ту же самую точку. Короче, функция:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
// Поиск 2D-точки пересечения отрезков AB и СВ function Find2DIntersection(const A, B, C, D: TPointF; out Intersection: TPointF; Epsilon: Double=1e-6): Boolean; var t, u, denom: Double; begin Intersection := TPointF.Create(MaxSingle, MaxSingle); Result := IsSegmentsIntersect(A, B, C, D, t, u, denom, Epsilon); // Если пересечение зафиксировано, параметр u не понадобится if Result then begin // По большому счёту, u тут вообще не нужен // Либо нет пересечения вообще, либо достаточно t Intersection.X := A.X + t*(B.X - A.X); Intersection.Y := A.Y + t*(B.Y - A.Y); end; end; |
Пересечение отрезков в 3D
Если отрезки пересекаются в плоскости XY, далеко не факт, что они пересекутся в координате Z.
Вот тут нам и понадобится параметр u. После того, как мы удостоверились, что отрезки пересекаются в плоскости XY, находим Z-отметки в точках пересечения и сравниваем их:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
// Поиск 3D-точек пересечения отрезков AB и СВ // Ориентир только на плоскость XY // Если там не пересекается, решаем, что и в 3D не пересекается function Find3DIntersection(const A, B, C, D: TPoint3D; out AB, CD: TPoint3D; Epsilon: Double=1e-6): Boolean; overload; var t, u, denom: Double; begin // Инициализируем пустыми значениями AB := TPoint3D.Create(MaxSingle, MaxSingle, 0); CD := TPoint3D.Create(MaxSingle, MaxSingle, 0); // Пытаемся найти параметры пересечения в плоскости XY Result := IsSegmentsIntersect(PPointF(@A)^, PPointF(@B)^, PPointF(@C)^, PPointF(@D)^, t, u, denom, Epsilon); // Нет пересечений в XY, выходим if not Result then exit; // Находим точки пересечения в плоскости XY AB.X := A.X + t*(B.X - A.X); AB.Y := A.Y + t*(B.Y - A.Y); CD := AB; // Вычисляем Z-координаты для обоих отрезков AB.Z := A.Z + t*(B.Z - A.Z); CD.Z := C.Z + u*(D.Z - C.Z); // Проверяем совпадение Z (копланарность) Result := SameValue(AB.Z, CD.Z, Epsilon) end; |
Приведённая выше функция рассматривает вариант с поиском пересечений только в XY. Но думать, что если нет пересечения в плоскости XY, то и в 3D нет — это неправильно (Данил, спасибо, что обратил внимание!).
Например, в случае, когда отрезки в плоскости XY параллельны, они вполне могут пересечься в плоскости XZ.
Чтобы однозначно найти точку пересечения отрезков в 3D, надо проверить все плоскости: XY, YZ и XZ. В функции ниже мы находим факт пересечения и определитель для каждой плоскости, затем выбираем максимальный по модулю определитель из допустимых. Просто при делении на большее число, получаем более точное значение (деление происходит при вычислении параметров t и u).
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 |
type TPlaneMode = (pmXY, pmYZ, pmXZ); // Поиск 3D-точек пересечения отрезков AB и СВ function Find3DIntersection(const A, B, C, D: TPoint3D; out AB, CD: TPoint3D; out Mode: TPlaneMode; Epsilon: Double=1e-6): Boolean; overload; var t, u: array[TPlaneMode] of Double; function Calculate(const P1, P2, P3, P4: TPoint3D; out R1, R2, R3, R4: Single): Boolean; begin R1 := P1.X + t[Mode]*(P2.X - P1.X); R2 := P1.Y + t[Mode]*(P2.Y - P1.Y); CD := AB; R3 := P1.Z + t[Mode]*(P2.Z - P1.Z); R4 := P3.Z + u[Mode]*(P4.Z - P3.Z); Result := SameValue(R3, R4, Epsilon) end; var i: TPlaneMode; Dmax: Double; denom: array[TPlaneMode] of Double; Res: array[TPlaneMode] of Boolean; begin // Инициализируем пустыми значениями AB := TPoint3D.Create(MaxSingle, MaxSingle, 0); CD := TPoint3D.Create(MaxSingle, MaxSingle, 0); // Пытаемся найти параметры пересечения в плоскости XY Res[pmXY] := IsSegmentsIntersect(PointF(A.X, A.Y), PointF(B.X, B.Y), PointF(C.X, C.Y), PointF(D.X, D.Y), t[pmXY], u[pmXY], denom[pmXY], Epsilon); // Пытаемся найти параметры пересечения в плоскости YZ Res[pmYZ] := IsSegmentsIntersect(PointF(A.Y, A.Z), PointF(B.Y, B.Z), PointF(C.Y, C.Z), PointF(D.Y, D.Z), t[pmYZ], u[pmYZ], denom[pmYZ], Epsilon); // Пытаемся найти параметры пересечения в плоскости XZ Res[pmXZ] := IsSegmentsIntersect(PointF(A.X, A.Z), PointF(B.X, B.Z), PointF(C.X, C.Z), PointF(D.X, D.Z), t[pmXZ], u[pmXZ], denom[pmXZ], Epsilon); // Ищем максимальный по модулю успешный определитель Dmax := 0; for i := Low(TPlaneMode) to High(TPlaneMode) do if Res[i] and (Abs(denom[i])>Dmax) then begin Mode := i; Dmax := Abs(denom[i]); end; // Проверяем, есть ли пересечения Result := not IsZero(Dmax, Epsilon); // Нет пересечений, выходим if not Result then exit; // Вычисляем координаты точек пересечения case Mode of pmXY: // Для плоскости XY Result := Calculate(A, B, C, D, AB.X, AB.Y, AB.Z, CD.Z); pmYZ: // Для плоскости YZ Result := Calculate( Point3D(A.Y, A.Z, A.X), Point3D(B.Y, B.Z, B.X), Point3D(C.Y, C.Z, C.X), Point3D(D.Y, D.Z, D.X), AB.Y, AB.Z, AB.X, CD.X); pmXZ: // Для плоскости XZ Result := Calculate( Point3D(A.X, A.Z, A.Y), Point3D(B.X, B.Z, B.Y), Point3D(C.X, C.Z, C.Y), Point3D(D.X, D.Z, D.Y), AB.X, AB.Z, AB.Y, CD.Y); else // Непонятное значение Result := False; end; end; |
Добавим ещё функцию, которая просто возвращает факт пересечения отрезков в 3D, и, в случае успеха, вернёт 3D-точку пересечения:
|
1 2 3 4 5 6 7 8 |
// Поиск 3D-точки пересечения отрезков AB и СВ function Find3DIntersection(const A, B, C, D: TPoint3D; out Intersection: TPoint3D; Epsilon: Double=1e-6): Boolean; overload; var Dummy: TPoint3D; Mode: TPlaneMode; begin Result := Find3DIntersection(A, B, C, D, Intersection, Dummy, Mode, Epsilon); end; |
Для копипаста
Для копипаста в спойлере модуль.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 |
//****************************************************************************** // Project: IP76.RU // Created: 2025-07-12 // Описание: Пересечение отрезков в 2D и 3D // https://ip76.ru/segments-intersection //****************************************************************************** unit IP76.Segments; interface uses System.Types, System.Math.Vectors; // Поиск 2D-точки пересечения отрезков AB и СВ function Find2DIntersection(const A, B, C, D: TPointF; out Intersection: TPointF; Epsilon: Double=1e-6): Boolean; // Поиск 3D-точек пересечения отрезков AB и СВ. // Ориентир только на плоскость XY, // Если там не пересекается, значит и в 3D не пересекается function Find3DIntersection(const A, B, C, D: TPoint3D; out AB, CD: TPoint3D; Epsilon: Double=1e-6): Boolean; overload; // Поиск 3D-точек пересечения отрезков AB и СВ type TPlaneMode = (pmXY, pmYZ, pmXZ); function Find3DIntersection(const A, B, C, D: TPoint3D; out AB, CD: TPoint3D; out Mode: TPlaneMode; Epsilon: Double=1e-6): Boolean; overload; // Поиск 3D-точки пересечения отрезков AB и СВ function Find3DIntersection(const A, B, C, D: TPoint3D; out Intersection: TPoint3D; Epsilon: Double=1e-6): Boolean; overload; implementation uses System.Math; // Проверка пересечения 2D-отрезков и возврат параметров пересечения function IsSegmentsIntersect(const A, B, C, D: TPointF; out t, u, denom: Double; Epsilon: Double): Boolean; var Dt, Du: Double; AB, CD, AC: TPointF; begin // Инициализируем параметры, даже если результат будет False Result := False; t := 0; u := 0; // Находим векторы AB := B-A; // Вектор AB CD := D-C; // Вектор CD AC := C-A; // Вектор AC // Находим определитель Δab denom := AB.X*CD.Y - AB.Y*CD.X; // Проверяем if IsZero(denom, Epsilon) then exit; // Параллельные сегменты, выходим // Находим определители Δt и Δu Dt := AC.X*CD.Y - AC.Y*CD.X; Du := AC.X*AB.Y - AC.Y*AB.X; // Решаем систему уравнений методом Крамера t := Dt/denom; u := Du/denom; // Проверка на попадание в отрезки Result := (t >= 0) and (t <= 1) and (u >= 0) and (u <= 1); end; // Поиск 2D-точки пересечения отрезков AB и СВ function Find2DIntersection(const A, B, C, D: TPointF; out Intersection: TPointF; Epsilon: Double=1e-6): Boolean; var t, u, denom: Double; begin Intersection := TPointF.Create(MaxSingle, MaxSingle); Result := IsSegmentsIntersect(A, B, C, D, t, u, denom, Epsilon); // Если пересечение зафиксировано, параметр u не понадобится if Result then begin // По большому счёту, u тут вообще не нужен // Либо нет пересечения вообще, либо достаточно t Intersection.X := A.X + t*(B.X - A.X); Intersection.Y := A.Y + t*(B.Y - A.Y); end; end; // Поиск 3D-точек пересечения отрезков AB и СВ. // Ориентир только на плоскость XY, // Если там не пересекается, решаем, что и в 3D не пересекается function Find3DIntersection(const A, B, C, D: TPoint3D; out AB, CD: TPoint3D; Epsilon: Double=1e-6): Boolean; overload; var t, u, denom: Double; begin // Инициализируем пустыми значениями AB := TPoint3D.Create(MaxSingle, MaxSingle, 0); CD := TPoint3D.Create(MaxSingle, MaxSingle, 0); // Пытаемся найти параметры пересечения в плоскости XY Result := IsSegmentsIntersect(PPointF(@A)^, PPointF(@B)^, PPointF(@C)^, PPointF(@D)^, t, u, denom, Epsilon); // Нет пересечений в XY, выходим if not Result then exit; // Находим точки пересечения в плоскости XY AB.X := A.X + t*(B.X - A.X); AB.Y := A.Y + t*(B.Y - A.Y); CD := AB; // Вычисляем Z-координаты для обоих отрезков AB.Z := A.Z + t*(B.Z - A.Z); CD.Z := C.Z + u*(D.Z - C.Z); // Проверяем совпадение Z (копланарность) Result := SameValue(AB.Z, CD.Z, Epsilon) end; // Поиск 3D-точек пересечения отрезков AB и СВ function Find3DIntersection(const A, B, C, D: TPoint3D; out AB, CD: TPoint3D; out Mode: TPlaneMode; Epsilon: Double=1e-6): Boolean; overload; var t, u: array[TPlaneMode] of Double; function Calculate(const P1, P2, P3, P4: TPoint3D; out R1, R2, R3, R4: Single): Boolean; begin R1 := P1.X + t[Mode]*(P2.X - P1.X); R2 := P1.Y + t[Mode]*(P2.Y - P1.Y); CD := AB; R3 := P1.Z + t[Mode]*(P2.Z - P1.Z); R4 := P3.Z + u[Mode]*(P4.Z - P3.Z); Result := SameValue(R3, R4, Epsilon) end; var i: TPlaneMode; Dmax: Double; denom: array[TPlaneMode] of Double; res: array[TPlaneMode] of Boolean; begin // Инициализируем пустыми значениями AB := TPoint3D.Create(MaxSingle, MaxSingle, 0); CD := TPoint3D.Create(MaxSingle, MaxSingle, 0); // Пытаемся найти параметры пересечения в плоскости XY Res[pmXY] := IsSegmentsIntersect(PointF(A.X, A.Y), PointF(B.X, B.Y), PointF(C.X, C.Y), PointF(D.X, D.Y), t[pmXY], u[pmXY], denom[pmXY], Epsilon); // Пытаемся найти параметры пересечения в плоскости YZ Res[pmYZ] := IsSegmentsIntersect(PointF(A.Y, A.Z), PointF(B.Y, B.Z), PointF(C.Y, C.Z), PointF(D.Y, D.Z), t[pmYZ], u[pmYZ], denom[pmYZ], Epsilon); // Пытаемся найти параметры пересечения в плоскости XZ Res[pmXZ] := IsSegmentsIntersect(PointF(A.X, A.Z), PointF(B.X, B.Z), PointF(C.X, C.Z), PointF(D.X, D.Z), t[pmXZ], u[pmXZ], denom[pmXZ], Epsilon); // Ищем максимальный по модулю успешный определитель Dmax := 0; for i := Low(TPlaneMode) to High(TPlaneMode) do if Res[i] and (Abs(denom[i])>Dmax) then begin Mode := i; Dmax := Abs(denom[i]); end; // Проверяем, есть ли пересечения Result := not IsZero(Dmax, Epsilon); // Нет пересечений, выходим if not Result then exit; // Вычисляем координаты точек пересечения case Mode of pmXY: // В плоскости XY Result := Calculate(A, B, C, D, AB.X, AB.Y, AB.Z, CD.Z); pmYZ: // В плоскости YZ Result := Calculate( Point3D(A.Y, A.Z, A.X), Point3D(B.Y, B.Z, B.X), Point3D(C.Y, C.Z, C.X), Point3D(D.Y, D.Z, D.X), AB.Y, AB.Z, AB.X, CD.X); pmXZ: // В плоскости XZ Result := Calculate( Point3D(A.X, A.Z, A.Y), Point3D(B.X, B.Z, B.Y), Point3D(C.X, C.Z, C.Y), Point3D(D.X, D.Z, D.Y), AB.X, AB.Z, AB.Y, CD.Y); else // Непонятное значение Result := False; end; end; // Поиск 3D-точки пересечения отрезков AB и СВ function Find3DIntersection(const A, B, C, D: TPoint3D; out Intersection: TPoint3D; Epsilon: Double=1e-6): Boolean; overload; var Dummy: TPoint3D; Mode: TPlaneMode; begin Result := Find3DIntersection(A, B, C, D, Intersection, Dummy, Mode, Epsilon); end; end. |
Демо приложение
Кнопками 2D и 3D меняем вид отображения.
В 2D можем таскать концы отрезков. В 3D крутим сцену (крутилка так себе, но представление о расположении даёт).
Кнопки 0..4 — это для иллюстрации статьи. Показывает векторы.
Если поставить галку на Signal a problem при отсутствии пересечения будет появляться вредный кот. Если в 2D кота нет, это не значит, что в 3D он не появится: Z-отметки отрезков могут не совпадать. В этом случае они будут красными.
Кнопка Set Z for AB подвинет отрезок AB по Z на величину разности Z-отметок в точке пересечения в плоскости XY. Таким образом можно установить полное пересечение по всем трём координатам и изгнать кота.
Есть одна особенность — 2D вид это на самом деле 3D, но с осью Z, направленной строго вертикально. Отображение отрезков в 2D учитывает Z координату и отрисовывается с учётом перспективы. Поэтому, когда захотите расположить их параллельно, и точки пересечения найти не удалось, визуально может казаться, что отрезки на самом деле пересекаются. Имеет смысл перейти в 3D вид и посмотреть, как там дела.
Скачать
Друзья, спасибо за внимание!
Исходник (zip) 90 Кб. Delphi XE 7
Исполняемый файл (zip) 1.24 Мб (Скомпилирован в XE 7 x64)
Касперский снова истерит на 32