diff --git a/22_segment_intersect/Bentley_Ottmann_tests.py b/22_segment_intersect/Bentley_Ottmann_tests.py new file mode 100644 index 0000000..9646637 --- /dev/null +++ b/22_segment_intersect/Bentley_Ottmann_tests.py @@ -0,0 +1,104 @@ +from entities import Point + + +def test(id) : + return { + 1: [2, + [Point(0, 0), Point(0, 24), + Point(-5, -25), Point(40, 40)], + []], + 2: [2, + [Point(0, 0), Point(3, 6), + Point(6, 6), Point(37, 37)], + []], + 3: [2, + [Point(0, 0), Point(2, 27), + Point(27, 27), Point(34, 34)], + []], + 4: [2, + [Point(0, 0), Point(2, 30), + Point(30, 30), Point(28, 28)], + []], + 5: [4, + [Point(10, 20), Point(50, 60), + Point(40, 30), Point(30, 60), + Point(50, 60), Point(50, 40), + Point(30, 20), Point(50, 40)], + [Point(35, 45), + Point(50, 60), + Point(50, 40), + Point(40, 30)]], + 6: [7, + [Point(4, 5), Point(1, 3), + Point(3, 3), Point(4, 5), + Point(4, 3), Point(4, 5), + Point(4, 5), Point(5, 3), + Point(4, 1), Point(3, 3), + Point(5, 3), Point(4, 1), + Point(4, 1), Point(4, 3)], + [Point(4, 5), + Point(3, 3), + Point(4, 3), + Point(5, 3), + Point(4, 1)]], + 7: [14, + [Point(-300, 0), Point(-300, 100), + Point(-300, 100), Point(-200, 200), + Point(-200, 200), Point(-100, 200), + Point(-100, 200), Point(0, 100), + Point(0, 100), Point(100, 200), + Point(100, 200), Point(200, 200), + Point(200, 200), Point(300, 100), + Point(300, 100), Point(300, 0), + Point(300, 0), Point(200, -100), + Point(200, -100), Point(100, -200), + Point(100, -200), Point(0, -300), + Point(0, -300), Point(-100, -200), + Point(-100, -200), Point(-200, -100), + Point(-200, -100), Point(-300, 0)], + [Point(-300, 100), + Point(-200, 200), + Point(-100, 200), + Point(0, 100), + Point(100, 200), + Point(200, 200), + Point(300, 100), + Point(300, 0), + Point(200, -100), + Point(100, -200), + Point(0, -300), + Point(-100, -200), + Point(-200, -100), + Point(-300, 0)]], + 8: [8, + [Point(-20, 0), Point(-20, 39), + Point(-20, 40), Point(29, 40), + Point(30, 40), Point(30, 11), + Point(30, 10), Point(-9, 10), + Point(-10, 10), Point(-10, 29), + Point(-10, 30), Point(19, 30), + Point(20, 30), Point(20, 21), + Point(20, 20), Point(0, 20)], + []], + 9: [4, + [Point(12, 455), Point(35, -61), + Point(123, 46), Point(-166, -181), + Point(-100, -13), Point(56, 700), + Point(-100, 1), Point(80, -1)], + [Point(63397847386456, -815531637627, 1000000000000), + Point(-48472157447999, 483023971644, 500000000000), + Point(32301958334230, -470021759269, 1000000000000), + Point(6673742589653, -4880485924390, 200000000000)]], + 10: [9, + [Point(-230, 48), Point(-84, 116), + Point(-137, 18), Point(-94, 342), + Point(-3000000000000000, 2704532564500349, 10000000000000), Point(-31, 0), + Point(-200, 200), Point(600, 500), + Point(300, 200), Point(500, 500), + Point(400, 200), Point(400, 500), + Point(500, 200), Point(300, 500), + Point(600, 200), Point(200, 500), + Point(200, 350), Point(600, 350)], + [Point(400, 350), + Point(-126628841820639, 96145470932852, 1000000000000)]] + }.get(id) diff --git a/22_segment_intersect/Line_segment_intersection.ipynb b/22_segment_intersect/Line_segment_intersection.ipynb new file mode 100644 index 0000000..4b8f3cd --- /dev/null +++ b/22_segment_intersect/Line_segment_intersection.ipynb @@ -0,0 +1,419 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import bintrees\n", + "\n", + "import entities\n", + "import solutions\n", + "import checker" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": false + }, + "source": [ + "# Пересечение множества отрезков\n", + "Перед нами стоит задача: дано множество $S$, состоящее из $n$ отрезков на плоскости. Нам нужно найти все точки их пересечений. \n", + "\n", + "Очевидно, что наивный алгоритм, перебирающий все пары отрезков, имеет сложность $O(n^2)$. В некотором смысле это оптимальная оценка, так как в случае, когда каждая пара отрезков из множества пересекается, любой алгоритм займет $\\Omega(n^2)$ времени. Однако на практике такие ситуации встречаются редко, обычно пересечений гораздо меньше. Поэтому хотелось бы уметь решать данную задачу быстрее в таких ситуациях. Другими словами, нужен алгоритм, время работы которого будет зависеть не только от $n$, но и от количества пересечений. Такие алгоритмы называются *output-sensitive*. \n", + "\n", + "В этом конспекте мы рассмотрим алгоритм Бентли-Оттмана (*англ. Bentley-Ottmann*), который позволяет решать поставленную задачу за время $O((n + I)\\log{n})$, где $I$ - количество персечений." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Перед тем, как приступить к изучению этого материала, настоятельно рекомендуется ознакомиться с конспектом про [афинное пространство](https://github.com/CT-18/cg/tree/5_affine_space), а именно с предикатом \"левый поворот\" и задачей о пересечении двух отрезков." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "## Алгоритм Бентли-Оттмана\n", + "\n", + "\n", + "Посмотрим на множество $S$. Нам нужно найти такое свойство или правило, по которому мы смогли бы отметать отрезки, которые точно не имеют общих точек, и не проверять их на пересечение.\n", + "\n", + "Начнем с достаточно простого наблюдения. Спроецируем все отрезки на ось $X$. Заметим, что отрезки, проекции которых не пересекаются, не могут иметь общих точек. А для отрезков, проекции которых накладываются друг на друга, существует параллельная оси $Y$ прямая, пересекающая их. \n", + "\n", + "Представим прямую $l$, скользящую по плоскости параллельно оси $Y$, начинающую свое движение из $-\\infty$. Будем называть её заметающей прямой (*sweep line*). Тогда, принимая во внимание наше наблюдение, мы можем утверждать, что пересекаются только такие отрезки, для которых существует положение прямой $l$, при котором она одновременно их пересекает. Таким образом, во время движения $l$ вдоль плоскости мы сможем отслеживать все отрезки, пересекающие её в один момент времени, и находить их пересечения.\n", + "\n", + "Будем называть статусом (*status*) заметающей прямой множество отрезков, пересекающих $l$ в данный момент времени, и событиями (*event points*) - точки, в которых статус меняется. Очевидно, что в нашем случае это точки начал и концов отрезков. \n", + "\n", + "Заметим, что алгоритм совершает какие-то дейтвия только при достижении $l$ события: обновляется статус прямой и происходят проверки отрезков на наличие пересечений. Например, если событие - \"левый конец отрезка\" $s_i$, то это значит, что $s_i$ начал пересекать заметающую прямую и должен быть добавен в статус. После добавления $s_i$ проверяется на пересечение с отрезками, уже имеющимися в статусе. При достижении события \"правый конец\", отрезок перестает пересекать $l$ и должен быть удален из статуса. Таким образом, мы проверяем на пересечение только такие отрезки, для которых существует вертикальная прямая, пересекающая их.\n", + "\n", + "Однако такой алгоритм все еще не output-sensitive, поскольку есть ситуации, при которых происходит проверка $O(n^2)$ пар отрезков, хотя пересечений на самом деле много меньше: простейший пример - множество горизонтальных отрезков, пересекающих ось $Y$. Проблема в том, что два отрезка, пересекающих заметающую прямую, могут находиться далеко друг от друга в вертикальном направлении. \n", + "\n", + "\n", + "Упорядочим отрезки, пересекающие $l$, снизу вверх: это поможет нам понять, какие отрезки находятся близко. Теперь будем проверять на пересечение только ближайших соседей, то есть каждый новодобавленный в статус отрезок будем проверять только с его непосредственными соседями снизу и сверху, если они есть. Позже у этого отрезка могут измениться соседи, и тогда они снова будут проверены на пересечение. Такой подход отразится на статусе заметающей прямой: теперь статус - это упорядоченная последовательность отрезков, пересекающих $l$. Соответственно, теперь он меняется не только при достижении концов отрезков, но и в точках пересечения, так как в них меняется порядок пересекающихся отрезков. Следовательно, мы должны проверить эти два отрезка на пересечение с их новыми соседями. Таким образом, теперь у нас появился третий тип событий - точки пересечения отрезков.\n", + "\n", + "Итак, мы разработали так называемый *plane sweep algorithm*. Для полного понимания принципа его работы, ниже представлена визуализация движения заметающей прямой:\n", + "" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "Теперь давайте подробнее разберем, что именно должно происходить при достижении заметающей прямой каждого типа события. Так как статус, как ранее было сказано, это упорядоченный набор пересекающих заметающую прямую отрезков в конкретный момент времени, то два отрезка являются соседними в статусе, если они являются ближайшими вдоль $l$.\n", + "1. Если событие - начало отрезка, то это значит, что новый отрезок начал пересекать заметающую прямую. Поэтому он должен быть добавлен в статус и проверен на пересечение с двумя своими непосредственными соседями сверху и снизу, если они есть. Нас интересуют только пересечения правее заметающей прямой, так как если пересение находится левее, то оно уже было обнаружено (этот достаточно очевидный факт будет доказан чуть позже). Например, если отрезки $s_k$ и $s_l$ были соседними в статусе, и в какой-то момент между ними появился новый отрезок $s_i$, то мы должны проверить, не пересекается ли $s_i$ с $s_k$ и $s_i$ с $s_l$. Если мы находим точку пересечения левее $l$, то мы должны запомнить это событие (позже мы разберем подробнее, что это значит) и впоследствии обработать его.

\n", + "2. Если событие - точка пересечения, то это значит, что эти отрезки поменяли свой порядок вдоль $l$. У каждого из них появился как максимум один новый сосед, с которым нужно проверить соответствующий отрезок на пересечение. Опять же нас интересуют только пересечения правее заметающей прямой (по той же причине). Приведём пример: предположим, что четыре отрезка идут в статусе в порядке $s_i$, $s_k$, $s_j$, $s_m$ в момент, когда $l$ достигает точки пересечения $s_k$ и $s_j$ (как на рисунке справа). В следующий момент эти два отрезка меняются в статусе местами, и мы должны проверить на перечечение $s_i$ с $s_j$ и $s_k$ с $s_m$. Если мы находим новые точки пересечений, то запоминаем их для дальнейшей обрабртки. Однако здесь стоит не забывать о том, что эти новые пересечения мы могли уже обнаружить ранее.

\n", + "3. Если событие - конец отрезка, то это значит, что он перестал пересекать заметающую прямую. Следовательно, его соседи сверху и снизу, если они оба есть, теперь являются ближайшими вдоль $l$ и должны быть проверены на пересечение. Если мы обнаружили пересечение правее заметающей прямой, то, опять же, запоминаем это событие для дальнейшей обработки. Предположим, что в статусе отрезки $s_k$, $s_j$ и $s_m$ идут в указанном порядке, когда заметающая прямая достигает канца отрезка $s_j$. Тогда отрезки $s_k$ и $s_m$ должны быть проверены на пересечение. Здесь снова стоит помнить, что новое событие \"пересечение отрезков\" могло быть уже обраружено ранее." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Теперь давайте обсудим, какие структуры данных нам понадобятся для реализации алгоритма. На самом деле их всего две: очередь событий $\\Omega$ и статус $T$. \n", + "\n", + "Как раньше было показано, алгоритм совершает какие-то действия только при достижении заметающей прямой какого-либо события. Следовательно, программно сэмитировать движение $l$ мы можем просто рассматривая события в том порядке, в котором прямая настигла бы их по ходу своего движения справа налево. Определим отношение \"меньше\" $\\prec$ для двух событий $p$ и $q$: будем считать, что событие \"левый конец\" меньше события \"пересечение\", а последнее меньше собятия \"правый конец\". Если события одного типа, то \n", + "$$ p < q \\Longleftrightarrow \\left[ \n", + " \\begin{align} \n", + " &p_x < q_x \\\\ \n", + " &p_x = q_x \\wedge p_y < q_y \\\\ \n", + " \\end{align} \n", + "\\right. $$\n", + "Таким образом, для очереди нам нужна такая структура данных, в которой мы могли бы хранить события в соответствии с определённым выше порядком. Кроме того надо помнить, что события \"пересечение отрезков\" появляются динамически, то есть мы должны уметь вставлять в очередь новые элементы, а также некоторые события могут быть вычислены больше одного раза, но в очереди при этом они должны встречаться лишь единожды. И ещё нам нужна операция поиска следующего события в очереди. Мы будем использовать сбалансированное бинарное дерево поиска, в котором обе операции - вставка нового элемента и поиск и удаление следующего элемента - занимают $O(\\log{m})$ времени, где $m$ - количество событий в $\\Omega$. \n", + "\n", + "Статус заметающей прямой нужен нам для поиска соседей данного отрезка. Структура данных статуса должна быть динамической, так как отрезки добавляются и удаляются из него в ходе работы алгоритма. Отрезки в статусе упорядочены, поэтому мы снова будем использовать сбалансированное бинарное дерево поиска. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Итак, мы придумали эффективный алгоритм, сильно сократив количество проверок (строгая оценка будет даначуть позже). Теперь нужно убедиться в том, что он по-прежнему находит все точки пересечений. Для начала уберем из рассмотрения перекрывающиеся (имеющие больше одной общей точки) и вертикальные отрезки, а также случай, когда больше двух отрезков пересекаются в одной точке. Эти три случая обрабатываются несложно, не выбиваясь из общего алгоритма. Однако для простоты рассуждений временно будем считать, что таких ситуаций быть не может, и покажем, что в таком предположении алгоритм работает корректно." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Теорема\n", + "> Пусть $s_i$ и $s_j$ - два невертикальных пересекающихся в точке $p$ отрезка, причем $p$ - единственная их общая точка. Также предположим, что нет третьего отрезка, проходящего через $p$. Тогда существует событие, которое случится раньше, чем заметающая прямая $l$ достигнет $p$, при котором $s_i$ и $s_j$ станут соседями и будут проверены на пересечение.\n", + "\n", + "$\\triangleright$
\n", + "Пусть $l$ находится строго левее $p$. Если $l$ достаточно близка к $p$, то $s_i$ и $s_j$ должны быть соседями (более формально мы должны взять такое положение $l$, при котором между ней и вертикальной прямой, параллельной оси $Y$ и проходящей через точку $p$, нет событий). Другими словами, существует такое положение заметающей прямой, в котором отрезки будут соседями. С другой стороны $s_i$ и $s_j$ не были соседями, когда $l$ начинала свое движение, поскольку она стартует левее всех точек. Следовательно, по ходу ее движение должно было произойти событие $q$, при котором $s_i$ и $s_j$ стали соседями и были проверены на пересечение.\n", + "
\n", + "$\\triangleleft$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Отметим, что началом отрезка будем считать тот его конец, который меньше второго в лексикографическом порядке.\n", + "\n", + "Теперь вернемся к трем случаям, о которых мы временно забывали, и рассмотрим, как обрабатывается каждый из них:\n", + "1. Если в множестве $S$ есть вертикальный отрезок $s_i$, то в очереди $\\Omega$ у нас будут лежать два события с одинаковыми $x$-координатами, отвечающие за начало и конец $s_i$. Благодаря вышеопределенному порядку $\\prec$, сначала будет рассмотрено событие \"начало отрезка\", так как его $y$-координата меньше, и отрезок будет проверен на пересечение с соседями. Через какое-то время будет рассмотрено событие \"конец отрезка\" и будут проведены соответствующие проверки. То есть случай вертикального отрезка ничем не отличается от случая невертикального, и нам не нужно совершать никаких дополнительных действий для его обработки. Для лучшего понимания можно просто представить, что заметающая прямая не вертикальная, а наклонена от вертикальной оси под очень небольшим углом влево, как на рисунке справа (такой прием называется *skew*). Тогда очевидно, что нижний конец вертикального отрезка всегда будет достигаться первым.

\n", + "2. Для случая, когда два отрезка накладываются друг на друга, мы просто заранее должны решить, что мы хотим от алгоритма в таких ситуациях. То есть для корректной обработки в этом случае нужно четкое определение того, что ожидается на выходе. Это может быть любая точка из общих точек двух отрезков, может быть конкретно обозначенная. Кроме того может быть кому-то понадобится возвращать список пар пересекающихся отрезков, в таком случае это опять же не будет выбиваться из общей картины.

\n", + "3. Теперь разберем случай, когда в одной точке пересекается больше двух отрезков. Так как у нас события отсортированы в порядке \"начало отрезка\" < \"пересечение отрезков\" < \"конец отрезка\", то случаи, когда в точке пересечения отрезков начинаются или заканчиваются другие отрезки, вписывается в общую картину и не требуют отдельного рассмотрения. Теперь вспомним, что в точке пересечения мы меняем порядок пересекающихся отрезков в статусе и проверяем их на пересечение с новыми соседями. Поэтому сколько бы отрезков не пересекалось в точке, благодаря таким перестановкам все пары пересекающихся отрезков будут найдены. Рассмотрим пример справа. Кратко опишем действия алгоритма: сначала в пустой статус добавляется отрезок $s_1$, далее - отрезок $s_2$, после проверки находится пересечение $s_2$ с $s_1$, далее в статус добавляется $s_3$ и находится пересечение $s_3$ с $s_1$. После этого добавляются $s_4$, $s_5$ и $s_6$, при добавлении последнего находятся пересечения $s_6$ с $s_1$ и $s_6$ с $s_3$. Если бы мы не рассматривали события \"пересечение отрезков\", то далее алгоритм прошел бы по событиям \"конец отрезка\", не нашел бы больше пересечений (поскольку мы берем только пересечения правее) и завершил свою работу. Но так как у нас три типа событий, то после добавления $s_6$ алгоритм перейдет к обработке события \"пересечение отрезков $s_2$ и $s_1$\", и проверит на пересечение $s_2$ и $s_3$ (найдет новое пересечение) и $s_1$ и $s_4$, далее будут обработаны события \"пересечение $s_1$ и $s_3$\" и \"пересечение $s_6$ и $s_1$\". Во время проверки последнего $s_6$ и $s_1$ поменяются местами и на пересечение будут проверены пары $s_6$ с $s_2$ (найдет новое пересечение) и $s_1$ с $s_3$. Далее будут проверены события \"пересечение $s_6$ и $s_3$\" (указанные отрезки поменяются местами и будут проверены пары $s_3$ с $s_1$ и $s_6$ с $s_5$), \"пересечение $s_6$ и $s_2$\", после чего отрезки в порядке $s_2$, $s_4$, $s_3$, $s_5$, $s_1$, $s_6$ будут удалены из статуса. Таким образом, алгоритм найдет шесть пар пересекающихся отрезков." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Реализация алгоритма\n", + "Ранее мы подробно обсудили идею алгоритма, теперь давайте соберем все вместе и напишем практическую реализацию. Здесь будет дан каркас основных моментов и описание их ожидаемой функциональности, читателю в качестве тренировки предлагается самому дописать недостающие части. Входные данные содержатся в файле *segments_set.in* в следующем формате: в первой строке задано $n$ - количество отрезков в множестве $S$, в каждой из следующих $n$ строк указаны координаты начала и конца соответствующего отрезка. Все координаты являются целыми числами.\n", + "\n", + "**Input:** Множество отрезков на плоскости $S$.
\n", + "**Output:** Множество точек пересечения отрезков из $S$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Теперь давайте опишем алгоритм *in general*, пока не рассматривая обработку событий и другие детали: \n", + "1. Инициализируем $status$ пустым множеством.\n", + "2. Инициализируем множество событий $events$ пустым множеством.\n", + "3. Вставим в $events$ события начал и концов отрезков из $S$.\n", + "4. Пока очередь событий не пуста, будем делать следующее:\n", + " 1. Достаем из очереди следующее событие $p$\n", + " 2. Обрабатываем $p$\n", + " \n", + "Для простоты реализации можно просто удалять пройденные события и на каждом шаге брать минимальное. Реализация этой части может быть следующей:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "Point = entities.Point\n", + "Segment = entities.Segment\n", + "Event = entities.Event\n", + "Type = entities.Type\n", + "\n", + "def find_intersections(segments):\n", + " status = bintrees.RBTree()\n", + " events = bintrees.RBTree()\n", + " intersection_points = set()\n", + " for s in segments:\n", + " events.insert(Event(s, Type.left), Event(s, Type.left))\n", + " events.insert(Event(s, Type.right), Event(s, Type.right)) \n", + " while not events.is_empty():\n", + " event, _ = events.pop_min()\n", + " intersection_points.update(handle_event_point(events, status, event)) # declared below\n", + " return intersection_points\n", + "\n", + "segments = []\n", + "with open(\"segments_set.in\") as inp:\n", + " n = int(inp.readline()) # number of segments\n", + " for i in range(0, n):\n", + " x1, y1, x2, y2 = map(int, inp.readline().split())\n", + " s = Segment(Point(x1, y1), Point(x2, y2), i)\n", + " segments.append(s)\n", + " plt.plot([x1, x2], [y1, y2], c = 'black')\n", + " plt.scatter(x1, y1, c = 'black', s = 30)\n", + " plt.scatter(x2, y2, c = 'black', s = 30)\n", + "\n", + "print(\"Found intersections:\", find_intersections(segments))\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Как видно из вышеприведенного кода, за обработку событий у нас будет отвечать функция *handle_event_point(event)*. В ней нужно будет обрабатывать каждый тип событий: проверять соседние отрезки на пересечение и искать новые события. Для реализации данной функции нам понадобится дополнительная, вычисляющая пересечение двух отрезков.\n", + "\n", + "Из курса линейной алгебры мы знаем, что точку пересечения $q$ двух непараллельных отрезков $a_1 b_1$ и $a_2 b_2$ можно найти следующим образом:\n", + "\n", + "\\begin{equation} \n", + "\\left( \\begin{array}{cc}\n", + " a_1.y - b_1.y & b_1.x - a_1.x\\\\\n", + " a_2.y - b_2.y & b_2.x - a_2.x\n", + "\\end{array} \\right) \n", + "\\cdot \n", + "\\left( \\begin{array}{c}\n", + " x\\\\\n", + " y\n", + "\\end{array} \\right) =\n", + "\\left( \\begin{array}{c}\n", + " (a_1.y - b_1.y) \\cdot a_1.x + (b_1.x - a_1.x) \\cdot a_1.y\\\\\n", + " (a_2.y - b_2.y) \\cdot a_2.x + (b_2.x - a_2.x) \\cdot a_2.y\n", + "\\end{array} \\right) \\end{equation}\n", + "\n", + "\\begin{equation} \\left( \\begin{array}{c}\n", + " x\\\\\n", + " y\n", + "\\end{array} \\right) = \n", + "\\frac{\n", + "\\left( \\begin{array}{cc}\n", + "b_2.x - a_2.x & a_1.x - b_1.x\\\\\n", + "b_2.y - a_2.y & a_1.y - b_1.y\n", + "\\end{array} \\right) \n", + "\\cdot\n", + "\\left( \\begin{array}{c}\n", + "b_1.x \\cdot a_1.y - b_1.y \\cdot a_1.x\\\\\n", + "b_2.x \\cdot a_2.y - b_2.y \\cdot a_2.x\n", + "\\end{array} \\right) \n", + "}{\n", + "\\left| \\begin{array}{cc}\n", + "a_1.y - b_1.y & b_1.x - a_1.x\\\\\n", + "a_2.y - b_2.y & b_2.x - a_2.x\n", + "\\end{array} \\right| } \\end{equation}\n", + "\n", + "Предположим, что два отрезка пересекаются (это, как мы знаем из конспекта про афинное пространство, легко проверяется с помощью предиката поворота и *bounding-box'а*). Сейчас Вам предлагается написать функцию, возвращающую точку пересечения этих отрезков. Для начала давайте будем работать в целых числах и все точки будем хранить в однородных координатах. Пусть на вход функции подаются четыре целых числа в $R^2$ - точки концов отрезков. Если отрезки имеют больше одной точки пересечения, то вывести можно любую." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def get_intersection_point(a, b, c, d): # концы отрезков в формате Point(x, y, 1)\n", + " # Заполните тело функции. Она должна возвращать точку пересечения\n", + " # отрезков ab и cd в однородных координатах.\n", + " return solutions.get_intersection_point(a, b, c, d) # точка пересечения в формате Point(x, y, det)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Чтобы убедиться в корректности написанной Вами функции, запустите ее на наборе тестов:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "checker.check_first_exercise(plt, get_intersection_point)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Теперь вернемся к обработке событий. Мы уже не раз повторили, что обработка события состоит из (не всегда) вставки или удаления нового отрезка из статуса и (всегда) из тестирования соседей на пересечение. Логика этой функции описана ниже.\n", + "\n", + "Для начала введем несколько обозначений:\n", + "* $U(p)$ - множество отрезков, верхний конец которых есть $p$\n", + "* $C(p)$ - множество отрезков, внутри которых содержится $p$\n", + "* $L(p)$ - множество отрезков, нижний конец которых есть $p$\n", + "\n", + "Теперь опишем последовательность операций, необходимых для обработки очередного события:\n", + "1. Найдем в статусе все отрезки, содержащие $p$. \n", + "2. Если $U(p) \\cup C(p) \\cup L(p)$ содержит больше одного отрезка, то $p$ - точка пересечения всех отрезков (каждого с каждым) в этом объединении.\n", + "3. Удалим из статуса отрезки из объединения $C(p) \\cup L(p)$.\n", + "4. Вставим отрезки из объединения $U(p) \\cup C(p)$ в статус. Порядок этих отрезков в $T$ должен соответствовать порядку, в котором отрезки пересекут заметающую прямую в следующий момент времени (то есть после пересечения). В частности, мы удаляем и заново вставляем отрезки из множества $C(p)$, чтобы обновился их порядок.\n", + "5. Если объединение $U(p) \\cup C(p) = \\varnothing$, то делаем следующее:\n", + " 1. Обозначим за $s_{left}$ левого соседа $p$ в $T$.\n", + " 2. Обозначим за $s_{right}$ правого соседа $p$ в $T$.\n", + " 3. Если оба отрезка $s_{left}$ и $s_{right}$ существуют, то мы должны проверить, не породят ли они новые события. Эту часть можно вынести в отдельную функцию, так как такая проверка нам будет нужна еще не раз. Сигнатура у неё будет следующая: *find_new_event(s_left, s_right, p)*. Подробнее опишем её чуть ниже.\n", + "6. Если $U(p) \\cup C(p) \\ne \\varnothing$, то:\n", + " 1. Обозначим за $s_{bound}$ самый левый отрезок из $U(p) \\cup C(p)$ в $T$.\n", + " 2. Обозначим за $s_{left}$ левого соседа $s_{bound}$ в $T$.\n", + " 3. Если оба отрезка существуют, то проверим их на новые события: вызовем *find_new_event*, передав в качестве параметров последовательно $s_{left}$, $s_{bound}$, $p$.\n", + " 4. Обозначим за $s_{bound}$ самый правый отрезок из $U(p) \\cup C(p)$ в $T$.\n", + " 5. Обозначим за $s_{right}$ правого соседа $s_{bound}$ в $T$.\n", + " 6. Если оба отрезка существуют, то проверим их на новые события: вызовем *find_new_event*, передав в качестве параметров последовательно $s_{bound}$, $s_{right}$, $p$.\n", + "\n", + "Последный момент, на котором стоит подробнее остановиться, это функция *find_new_event(s_left, s_right, p)*. Она должна проверять отрезки $s_{left}$ и $s_{right}$ на пересечение, и если оно действительно есть, то смотреть, что точка пересечния находится справа от текущего положения заметающей прямой или на ней, но строго выше $p$. Если такое событие найдено и его еще нет в очереди (важно, что второй раз функция его добавлять не должна!), то оно должна быть добавлено.\n", + " \n", + "Сейчас Вам предлагается написать данную функцию по описанному шаблону." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [], + "source": [ + "def handle_event_point(events, status, event):\n", + " # Заполните тело функции. Она должна возвращать множество уникальных \n", + " # точек пересечения отрезков из входного набора отрезков.\n", + " return set();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Чтобы проверить правильность алгоритма в целом, используйте чекер:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "checker.check_Bentley_Ottmann_algorithm(plt, find_intersections)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Теорема (корректность алгоритма)\n", + ">Алгоритм находит все точки пересечения и отрезки, которые их содержат. \n", + "\n", + "$\\triangleright$
\n", + "Воспользуемся индукцией по событиям, отсортированным в вышеуказанном порядке. Предположим что все события $q : q < p$ были обработаны корректно. Докажем, что $p$ также будет корректно обработано и все отрезки, содержащие $p$, будут найдены. \n", + "1. Если событие $p$ - точка пересечения отрезков, то по предыдущей теореме найдется предшествующее $p$ событие $q$ такое, что при его обработке отрезки, пересекающиеся в $p$, будут соседними в статусе. Следовательно, в этом случае $p$ будет обнаружена.\n", + "2. Если $p$ - начало или конец отрезка, то она была добавлена в $events$ в самом начале работы алгоритма. Все содержащие её отрезки из статуса в момент обработки $p$ будут найдены. Для остальных отрезков, содержащих $p$, верно, что $p$ — их левый конец, и они будут найдены, так как мы храним их вместе с $p$ в $events$.\n", + "
\n", + "$\\triangleleft$\n", + "\n", + "Следовательно, после того, как заметающая прямая пройдет всю плоскость, а точнее после того, как будет рассмотрено последнее событие, мы вычислим все точки пересечений отрезков из множества $S$. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Теорема (оценка времени работы)\n", + ">Время работы алгоритма - $O((n + I)\\log{n})$, где $I$ - количество пересечений.\n", + "\n", + "$\\triangleright$
\n", + "На первом шаге алгорим строит множество событий $events$, при использовании сбалансированного бинарного дерева в качестве структуры данных это выполняется за $O(n\\log{n})$. Инициализация статуса занимает константное время. \n", + "\n", + "Далее на каждой итерации мы обрабатываем по одному событию. Обработка может включать в себя удаление отрезка из статуса за $O(\\log{n})$ и проверку на пересечение отрезка с не более, чем двумя соседями, что может вызвать вставку новых событий (снова $O(\\log{n})$). Всего при обработке события $p$ мы $m(p) = \\left| L(p) \\cup C(p) \\cup U(p) \\right|$ раз выполняем операции вставки, удаления и поиска отрезка в статусе, где $L(p)$ - множество отрезков, верхний конец которых есть $p$, $C(p)$ - множество отрезков, внутри которых содержится $p$, и $U(p)$ - множество отрезков, нижний конец которых есть $p$. Пусть $m = \\sum_{p}^{}m(p)$. Тогда время работы алгоритма - $O(n\\log{n} + m\\log{n})$. Осталось показать, что $m = n + I$. \n", + "\n", + "Рассмотрим планарный граф, вершинами которого являются концы отрезков, а также их точки пересечения, а ребрами - части отрезков, их соединяющие. Рассмотрим событие $p$. Ему соответствует вершина графа со степенью $O(m(p))$. Следовательно, $m$ не превосходит суммы степеней всех вершин графа: $m = O(\\sum_{p}^{}deg(p))$. Каждое ребро графа добавляет $+1$ к степени ровно двух вершин, из чего следует, что $m = O(2\\left| E \\right|) = O(\\left| E \\right|)$, где $\\left| E \\right|$ - количество ребер в графе. Пусть $\\left| V \\right|$ - количество вершин в графе, тогда $\\left| V \\right| \\leqslant 2n + I$. Мы знаем, что в планарном графе $\\left| E \\right| = O(\\left| V \\right|)$. Таким образом, $O(\\left| E \\right|) = O(\\left| V \\right|) = O(2n + I)$. Следовательно, время работы алгоритма - $O((n + I)\\log{n})$.\n", + "
\n", + "$\\triangleleft$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Требуемая память\n", + "Благодаря данному алгоритму можно сократить количество проверяемых отрезков с $n^2$ до $n + I$, что является отличным результатом. Теперь оценим затраты по памяти. Очевидно, что статус не может занимать больше $O(n)$ памяти, а вот очередь событий вполне может занимать $O(n + I)$, то есть если все пересечения окажутся правее заметающей прямой, то для очереди потребуется $O(n^2)$ памяти. Чтобы избежать этой оценки, поступают следующим образом: если два отрезка перестали быть соседями, то событие, отвечающее за их точку пересечения, удаляется из статуса. Таким образом, мы будем хранить не больше $O(n)$ событий, а на время работы алгоритма такая оптимизация не повлияет." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Заключение\n", + "Итак, мы научились решать задачу о нахождении всех точек пересечения множества отрезков за $O((n + I)\\log{n})$ времени и $O(n)$ памяти. Понятно, что в худшем случае данный алгоритм работает за $O(n^2\\log{n})$. Существуют алгоритмы, которые имеют лучшую временную оценку, например, в следующей статье описывается решение поставленной задачи за время $O(n\\log{n} + I)$ с сохранением оценки по памяти $O(n)$:\n", + "\n", + "* [Ivan J. Balaban - An optimal algorithm for finding segments intersections](https://pdfs.semanticscholar.org/f523/84b8cc2b44f91c1e4048c7bc0021f6e01b01.pdf)\n", + "\n", + "Однако алгоритмы Бентли-Оттмана и Ивана Балабана сложны в реализации, поэтому на практике используются достаточно редко. Чаще всего эта задача решается с использованием [квадродеревьев](https://github.com/CT-18/cg/tree/1_skip_quadtree), так как такой подход имеет ту же асимптотику и при этом более прост в реализации." + ] + } + ], + "metadata": { + "anaconda-cloud": {}, + "kernelspec": { + "display_name": "Python [default]", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} diff --git a/22_segment_intersect/checker.py b/22_segment_intersect/checker.py new file mode 100644 index 0000000..611964e --- /dev/null +++ b/22_segment_intersect/checker.py @@ -0,0 +1,88 @@ +import numpy as np +from entities import Point, Segment + +import exercise_1_tests +import Bentley_Ottmann_tests + + +def check_first_exercise(plt, get_intersection_point): + test = exercise_1_tests.test + + def calculate_det(a, b): + return a.x * b.y - b.x * a.y + + def bounding_box(p, a, b): + if a > b: + a, b = b, a + return a <= p <= b + + def check_overlap(p, a, b): + orientation = np.sign(calculate_det(b - a, p - a)) + if orientation != 0: + return False + return bounding_box(p.x, a.x, b.x) and bounding_box(p.y, a.y, b.y) + + def drow_points(axis, points, color): + for point in points: + axis.scatter(point.x, point.y, c=color, s=40) + + def drow_segments(axis, a, b, color): + axis.plot([a.x, b.x], [a.y, b.y], c=color) + + f, axes = plt.subplots(2, 4, figsize=(11, 6)) + print("Part 1. Intersection tests") + for i, axis in zip(range(1, 9), axes.reshape(8)): + axis.set_title("Test " + str(i)) + a, b, c, d, answ = test(i) + output = get_intersection_point(a, b, c, d) + drow_points(axis, [a, b], 'r') + drow_segments(axis, a, b, 'r') + drow_points(axis, [c, d], 'b') + drow_segments(axis, c, d, 'b') + if output == answ: + print("Test", i, "Ok") + else: + print("Test", i, "Failed:") + print("\tx={}, y={}, det={}".format(answ.x, answ.y, answ.det), "expected but\n", + "\tx={}, y={}, det={}".format(output.x, output.y, output.det), "found") + print("Part 2. Overlap tests") + for i in range(9, 21): + a, b, c, d = test(i) + output = get_intersection_point(a, b, c, d) + if check_overlap(output, a, b) and check_overlap(output, c, d): + print("Test", i, "Ok") + else: + print("Test", i, "Fail") + print("\ta={}, b={}, c={}, d={}".format(a, b, c, d)) + print("\tYour answer:", output) + + plt.show() + + +def check_Bentley_Ottmann_algorithm(plt, find_intersections): + test = Bentley_Ottmann_tests.test + + f, axes = plt.subplots(2, 5, figsize=(12, 6)) + for i, axis in zip(range(1, 11), axes.reshape(10)): + axis.set_title("Test " + str(i)) + input = test(i) + n = input[0] # number of segments + segments = [] + for j in range(0, 2 * n, 2): + p1, p2 = input[1][j], input[1][j + 1] + s = Segment(p1, p2, j) + segments.append(s) + axis.plot([p1.x, p2.x], [p1.y, p2.y], c='black') + axis.scatter(p1.x, p1.y, c='black', s=30) + axis.scatter(p2.x, p2.y, c='black', s=30) + intersection_points = find_intersections(segments) + intersection_points_answ = set() + intersection_points_answ.update(input[2]) + if intersection_points == intersection_points_answ: + print("Test", i, "Ok") + else: + print("Test", i, "Fail:") + print("\t", intersection_points_answ, "expected but") + print("\t", intersection_points if len(intersection_points) > 0 else "{}", "found") + + plt.show() diff --git a/22_segment_intersect/entities.py b/22_segment_intersect/entities.py new file mode 100644 index 0000000..5fb9000 --- /dev/null +++ b/22_segment_intersect/entities.py @@ -0,0 +1,113 @@ +from fractions import * +from enum import Enum +Type = Enum('Type', 'left right inter') + + +class Point: + __slots__ = ("x", "y", "det") + + def __init__(self, x=0, y=0, det=1): + self.x = x + self.y = y + self.det = det + + def __add__(self, p): + return Point(Fraction(self.x, self.det) + Fraction(p.x, p.det), + Fraction(self.y, self.det) + Fraction(p.y, p.det)) + + def __sub__(self, p): + return Point(Fraction(self.x, self.det) - Fraction(p.x, p.det), + Fraction(self.y, self.det) - Fraction(p.y, p.det)) + + def __neg__(self): + return Point(-self.x, -self.y, self.det) + + def __lt__(self, p): + return Fraction(self.x, self.det) < Fraction(p.x, p.det) or \ + Fraction(self.x, self.det) == Fraction(p.x, p.det) and Fraction(self.y, self.det) < Fraction(p.y, p.det) + + def __eq__(self, p): + return Fraction(self.x, self.det) == Fraction(p.x, p.det) and \ + Fraction(self.y, self.det) == Fraction(p.y, p.det) + + def __gt__(self, p): + return not (self < p or self == p) + + def __hash__(self): + return hash(tuple([Fraction(self.x, self.det), Fraction(self.y, self.det)])) + + def __repr__(self): + return "(%r, %r)" % (str(Fraction(self.x, self.det)), str(Fraction(self.y, self.det))) + + +class Segment: + __slots__ = ("a", "b", "id") + + def __init__(self, a, b, id): + self.a = min(a, b) + self.b = max(a, b) + self.id = id + + def __lt__(self, s): + if s is None: + return False + if self.id == s.id: + return False + if self.a == s.a: + return self.b < s.b + return self.a < s.a + + def __eq__(self, s): + return self.id == s.id + + def __gt__(self, p): + return not (self < p or self == p) + + def __repr__(self): + return "[%r, %r]" % (self.a, self.b) + + +class Type(Enum): + left = 0 + right = 1 + inter = 2 + + def __lt__(self, other): + return self.value < other.value + + +class Event: + __slots__ = ("segment", "t", "addit_segment", "inter_point") # addit_segment & inter_point are used only for type = inter + + def __init__(self, segment, t, addit_segment = None, inter_point = None): + self.segment = segment + self.t = t + self.addit_segment = addit_segment + self.inter_point = inter_point + + def __lt__(self, e): + a = self.segment.a if self.t == Type.left else \ + self.segment.b if self.t == Type.right else self.inter_point + b = e.segment.a if e.t == Type.left else \ + e.segment.b if e.t == Type.right else e.inter_point + if a == b: + if self.t == e.t: + if self.segment == e.segment: + return self.addit_segment < e.addit_segment + return self.segment < e.segment + return self.t < e.t + return a < b + + def __eq__(self, e): + return self.t == e.t and \ + self.segment.id == e.segment.id and \ + (self.addit_segment is None and e.addit_segment is None or + self.addit_segment is not None and e.addit_segment is None and self.addit_segment == e.addit_segment) + + def __gt__(self, e): + return not (self == e or self < e) + + def __repr__(self): + if type != Type.inter: + return "Event [%r, type = %r]" % (self.segment, self.t) + return "Event inter [%r, %r, point = %r]" % (self.segment, self.addit_segment, self.inter_point) \ No newline at end of file diff --git a/22_segment_intersect/exercise_1_tests.py b/22_segment_intersect/exercise_1_tests.py new file mode 100644 index 0000000..5b582fa --- /dev/null +++ b/22_segment_intersect/exercise_1_tests.py @@ -0,0 +1,27 @@ +from entities import Point + + +def test(id): + return { + 1: [Point(1, 1), Point(3, 3), Point(3, 3), Point(3, 2), Point(-6, -6, -2)], + 2: [Point(1, 1), Point(4, 5), Point(2, 1), Point(2, 5), Point(24, 28, 12)], + 3: [Point(-18, -34), Point(62, 15), Point(83, 7), Point(-8, 16), Point(277978, 51274, 5179)], + 4: [Point(0, 0), Point(1, 4), Point(2, 2), Point(0, 3), Point(6, 24, 9)], + 5: [Point(0, 0), Point(7, 1), Point(0, 1), Point(7, 0), Point(-49, -7, -14)], + 6: [Point(1, 2), Point(4, 2), Point(4, 2), Point(7, 2), Point(4, 2, 1)], + 7: [Point(0, 0), Point(5, 2), Point(3, 0), Point(2, 2), Point(30, 12, 12)], + 8: [Point(10, 1), Point(100, 80), Point(30, 2), Point(-60, 95), Point(330300, 169530, 15480)], + + 9: [Point(-1, 1), Point(3, 1), Point(1, 1), Point(5, 1)], + 10: [Point(-2, -2), Point(-2, 6), Point(-2, 3), Point(-2, 10)], + 11: [Point(1, 1), Point(5, 5), Point(3, 3), Point(7, 7)], + 12: [Point(10, -10), Point(1, -1), Point(-5, 5), Point(5, -5)], + 13: [Point(-80, -80), Point(30, -80), Point(30, -80), Point(-80, -80)], + 14: [Point(8, -70), Point(8, 70), Point(8, -70), Point(8, 70)], + 15: [Point(-15, -15), Point(32, 32), Point(32, 32), Point(-15, -15)], + 16: [Point(-15, 132), Point(14, 14), Point(14, 14), Point(-15, 132)], + 17: [Point(4, 18), Point(4, 46), Point(4, 46), Point(4, 43)], + 18: [Point(6, -17), Point(8, -17), Point(33, -17), Point(6, -17)], + 19: [Point(3, 18), Point(15, 1), Point(27, -16), Point(3, 18)], + 20: [Point(-1, -1), Point(8, 8), Point(-1, -1), Point(101, 101)] + }.get(id) diff --git a/22_segment_intersect/images/intersection_change_status.png b/22_segment_intersect/images/intersection_change_status.png new file mode 100644 index 0000000..814d61b Binary files /dev/null and b/22_segment_intersect/images/intersection_change_status.png differ diff --git a/22_segment_intersect/images/intersection_point.png b/22_segment_intersect/images/intersection_point.png new file mode 100644 index 0000000..5e75879 Binary files /dev/null and b/22_segment_intersect/images/intersection_point.png differ diff --git a/22_segment_intersect/images/lower_endpoint.png b/22_segment_intersect/images/lower_endpoint.png new file mode 100644 index 0000000..4dd98a5 Binary files /dev/null and b/22_segment_intersect/images/lower_endpoint.png differ diff --git a/22_segment_intersect/images/several_intersections_case.png b/22_segment_intersect/images/several_intersections_case.png new file mode 100644 index 0000000..b01051c Binary files /dev/null and b/22_segment_intersect/images/several_intersections_case.png differ diff --git a/22_segment_intersect/images/skew.png b/22_segment_intersect/images/skew.png new file mode 100644 index 0000000..d87561f Binary files /dev/null and b/22_segment_intersect/images/skew.png differ diff --git a/22_segment_intersect/images/sweep_line.png b/22_segment_intersect/images/sweep_line.png new file mode 100644 index 0000000..f3c3d4e Binary files /dev/null and b/22_segment_intersect/images/sweep_line.png differ diff --git a/22_segment_intersect/images/sweep_line_visualization.gif b/22_segment_intersect/images/sweep_line_visualization.gif new file mode 100644 index 0000000..e6fda0b Binary files /dev/null and b/22_segment_intersect/images/sweep_line_visualization.gif differ diff --git a/22_segment_intersect/images/upper_endpoint.png b/22_segment_intersect/images/upper_endpoint.png new file mode 100644 index 0000000..46de0d9 Binary files /dev/null and b/22_segment_intersect/images/upper_endpoint.png differ diff --git a/22_segment_intersect/segments_set.in b/22_segment_intersect/segments_set.in new file mode 100644 index 0000000..79a5a18 --- /dev/null +++ b/22_segment_intersect/segments_set.in @@ -0,0 +1,10 @@ +9 +1 3 1 -2 +6 1 1 0 +3 3 1 1 +5 -2 1 3 +8 3 4 9 +2 1 0 6 +2 3 3 2 +-1 3 5 6 +1 1 8 8 \ No newline at end of file diff --git a/22_segment_intersect/solutions.py b/22_segment_intersect/solutions.py new file mode 100644 index 0000000..a3a1eef --- /dev/null +++ b/22_segment_intersect/solutions.py @@ -0,0 +1,37 @@ +import numpy as np +from entities import Point + + +def point_belongs_to_segment(p, a, b): # as we know both segments lay on the same line + left, right, p_ = 0, 0, 0 + if a.x == b.x: + left, right, p_ = a.y, b.y, p.y + else: + left, right, p_ = a.x, b.x, p.x + if left > right: + left, right = right, left + return left <= p_ <= right + + +def get_intersection_point(a, b, c, d): + left = [[d.x - c.x, a.x - b.x], + [d.y - c.y, a.y - b.y]] + right = [[b.x * a.y - b.y * a.x], + [d.x * c.y - d.y * c.x]] + numerator = np.matmul(left, right) + denominator = (a.y - b.y) * (d.x - c.x) - (b.x - a.x) * (c.y - d.y) + if denominator != 0: + return Point(numerator[0][0], numerator[1][0], denominator) + p = None + if point_belongs_to_segment(a, c, d): + p = a + elif point_belongs_to_segment(b, c, d): + p = b + elif point_belongs_to_segment(c, a, b): + p = c + elif point_belongs_to_segment(d, a, b): + p = d + if p is None: + raise Exception('Incorrect test!') + return p +