2D Geometric Constraint Solver.

Содержание

  1. Вступление
  2. Постановка задачи
  3. Математическая нотация
  4. Пример
  5. Оптимизация
  6. Параметризация геометрических примитивов
  7. Замена переменных
  8. Имплементация геометрических ограничений

Вступление

Я провожу кучу времени в MCAD программах, таких как SolidWorks, Onshape, Fusion 360 и т.д. Мне всегда было интересно, как они устроены “под капотом”, поэтому я решил создать свой собственный параметрический 2D-редактор чертежей с ограниченным набором геометрических примитивов (только отрезки и дуги окружностей) и базовыми геометрическими ограничениями. Главная проблема, которую нужно решить, это Geometric constraint solving, и для фана я решил разработать свой метод решения этой проблемы, не изучая предварительно существующие подходы.

Вот что в итоге у меня получилось:

Программа написана на Python, исходный код доступен на GitHub. Ниже я постараюсь объяснить, как все это устроено.

Постановка задачи

Основные функции разрабатываемого 2D редактора, которые должны поддерживаться:

  1. Добавление и удаление отрезков и дуг
  2. Передвижение и изменение с помощью мыши существующих отрезков и дуг, при этом программа должна автоматически следить за тем, чтобы существующие геометрические ограничения не нарушались
  3. Наложение новых геометрические ограничений на линии и дуги и удаление существующих геометрические ограничений

Геометрические ограничения, которые должны поддерживаться:

  1. COINCIDENCE: совпадение точек
  2. PARALLELITY: параллельность отрезков
  3. PERPENDICULARITY: перпендикулярность отрезков
  4. EQUAL_LENGTH_OR_RADIUS: равная длина отрезков или равный радиус дуг
  5. FIXED: неподвижность точек
  6. HORIZONTALITY: горизонтальность точек
  7. VERTICALITY: вертикальность точек
  8. TANGENCY: касательность дуги и отрезка или касательность двух дуг
  9. CONCENTRICITY: концентричность дуг

Для того, чтобы разобраться, как все это работает, рассмотрим пример:

Математическая нотация

Вот некоторые обозначения, которые будут использоваться далее:

  1. $x$ – скаляр
  2. $\vec{y}$ – вектор в математическом смысле
  3. $\vec{f}()$ – вектор-функция
  4. $\overline{z}$ – вектор в геометрическом смысле

Пример

Допустим, у нас есть два отрезка, $AB$ и $BC$. Также у нас есть вот такие геометрические ограничения:

  1. $AB$ и $BC$ перпендикулярны
  2. $AB$ и $BC$ имеют равную длину

Допустим, в момент времени $t_0$ система находилась в решенном состоянии, все геометрические ограничения были удовлетворены, этому состоянию соответствуют отрезки $A_0B_0$ и $B_0C_0$ на картинке ниже. Затем, с помощью мыши мы мгновенно изменили положение точки $C_0$, принадлежащей отрезку $B_0C_0$, после чего эта точка перешла в положение $C’$:

Далее должен отработать солвер, который должен найти новые (соответствующие моменту времени $t_1$) положения всех отрезков, что фактически означает поиск точек $A_1$, $B_1$ и $C_1$. У солвера есть две задачи:

  1. Необходимо, чтобы точка $C_1$ была как можно ближе к точке $C’$. Иногда (но не всегда) возможно обеспечить совпадение этих точек
  2. Необходимо, чтобы все существующие геометрические ограничения были удовлетворены

Вот так все будет выглядеть после того, как солвер завершит свою работу:

Запишем определенные выше задачи солвера более формальным языком:

\[dist(C', C_1) \rightarrow min\] \[\begin{cases} A_1B_1 \perp B_1C_1\\ \norm{\overline{A_1B_1}} = \norm{\overline{B_1C_1}} \end{cases} \Rightarrow \begin{pmatrix} dot(\overline{A_1B_1}, \overline{B_1C_1}) \\ \norm{\overline{A_1B_1}} - \norm{\overline{B_1C_1}} \end{pmatrix} =0\]

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

\[\begin{array}{rl} \min _\vec{x} & f(\vec{x}) \\ \text { subject to } & \vec{h}(\vec{x}) \geq 0 \\ & \vec{g}(\vec{x})=0 . \end{array}\]

Функция $\vec{h}(\vec{x})$ описывает ограничения-неравенства, $\vec{g}(\vec{x})$ – ограничения-равенства. У нас есть только ограничения-равенства, поэтому выражения принимают следующий вид:

\[\begin{array}{rl} \min _\vec{x} & f(\vec{x}) \\ \text { subject to } & \vec{g}(\vec{x})=0. \end{array}\]

Явно выразим вектор $\vec{x}$:

\[\vec{x}= \begin{pmatrix} A_x & A_y & B_x & B_y & C_x & C_y \end{pmatrix} ^\intercal\]

Явно выразим $f(\vec{x})$ через компоненты вектора $\vec{x}$:

\[f(\vec{x})=dist(C', C_1)=\sqrt{(C'_x - C_{1x})^2 + (C'_y - C_{1y})^2}\]

Явно выразим $\vec{g}(\vec{x})$ через компоненты вектора $\vec{x}$:

\[\vec{g}(\vec{x})= \begin{pmatrix} dot(\overline{A_1B_1}, \overline{B_1C_1}) \\ \norm{\overline{A_1B_1}} - \norm{\overline{B_1C_1}} \end{pmatrix} = \begin{pmatrix} (B_{1x} - A_{1x}) \cdot (C_{1x} - B_{1x}) + (B_{1y} - A_{1y}) \cdot (C_{1y} - B_{1y}) \\ \sqrt{(B_{1x} - A_{1x})^2 + (B_{1y} - A_{1y})^2} - \sqrt{(C_{1x} - B_{1x})^2 + (C_{1y} - B_{1y})^2} \end{pmatrix}\]

Теперь, когда мы сформулировали задачу нелинейного программирования, можно поговорить о том, как ее решать.

Оптимизация

Я решил не имплементровать алгоритмы оптимизации самостоятельно, а воспользоваться готовым Python пакетом scipy.optimize. В частности, из этого пакета нам потребуется функция minimize. Эта функция интересна тем, что принимает на вход параметр, определяющий тип солвера. Я попробовал разные типы, наиболее хорошо показал себя вариант “SLSQP”.

Теперь, думаю, основной принцип работы программы понятен. Но нужно также обсудить некоторые тонкости:

Параметризация геометрических примитивов

Из вышеприведенного примера понятно, что нам необходимо уметь преобразовывать нашу систему из набора отрезков и дуг в вектор чисел ($\vec{x}$) и обратно. Такое преобразование называется параметризацией. Чтобы параметризовать систему, нужно параметризовать ее составляющие и объединить результаты. Таким образом, нам нужно уметь параметризовывать отрезки и дуги.

Для отрезка можно использовать самый очевидный и простой способ, такой же, как в примере. Пусть у нас есть отрезок $AB$, тогда вектор $\vec{x}$, соответствующий ему, будет иметь вид:

\[\vec{x}_{AB}= \begin{pmatrix} A_x & A_y & B_x & B_y \end{pmatrix} ^\intercal\]

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

В данном случае, чтобы определить дугу, используются точки ее начала и конца, а также скаляр $d$. Важно отметить, что дуга всегда ориентирована по часовой стрелке. Вектор $\vec{x}$, соответствующий дуге, будет иметь вид:

\[\vec{x}_{arc}= \begin{pmatrix} P_{1x} & P_{1y} & P_{2x} & P_{2y} & d \end{pmatrix} ^\intercal\]

Теперь разберемся, что такое скаляр $d$. Точка $O$ является центром дуги, точка $M$ является серединой отрезка $P_1P_2$. Вектор $\overline{a}$ равен вектору $\overline{MO}$. Вектор $\overline{n}$ является нормалью к отрезку $P_1P_2$, он получен путем нормализации вектора $\overline{P_1P_2}$ и поворота его на 90° против часовой стрелки. Скаляр $d$ вычисляется так:

\[d=dot(\overline{a}, \overline{n})\]

Восстановление дуги из вектора $\vec{x}_{arc}$ также не представляет проблем. Так как точки $P_1$ и $P_2$ и так известны, необходимо найти только точку $O$. Для того, чтобы это сделать, нужно сначала найти точку $M$ как центр отрезка $P_1P_2$ а потом просто сместить эту точку на вектор $(\overline{n} \cdot d)$.

Замена переменных

Геометрические ограничения, которые встречается очень часто, это совпадение точек, то есть COINCIDENCE, а также горизонтальность и вертикальность точек, HORIZONTALITY и VERTICALITY. Снова рассмотрим небольшой пример:

У нас есть два отрезка, $AB$ и $CD$. Мы хотим наложить геометрическое ограничение совпадения точек на точки $B$ и $C$.

Самым простым и наивным способом это сделать является просто добавление двух элементов в $\vec{g}(\vec{x})$:

\[\vec{g}(\vec{x})= \begin{pmatrix} B_x - C_x \\ B_y - C_y \end{pmatrix}\]

Проблема этого подхода заключается в тот, что геометрических ограничений типа COINCIDENCE может быть очень много. Если учитывать их приведенным выше способом, это сильно усложнит жизнь солвера, так как в $\vec{g}(\vec{x})$ будет очень много элементов. Однако, есть способ этого избежать:

Изначально, вектор $\vec{x}$ для системы, изображенной выше, будет иметь вид:

\[\vec{x}= \begin{pmatrix} A_x & A_y & B_x & B_y & C_x & C_y & D_x & D_y \end{pmatrix} ^\intercal\]

Учитывая, что точки $B$ и $C$ совпадают, можно сделать такую замену:

\[\begin{cases} B_x = C_x = \alpha_0\\ B_y = C_y = \alpha_1\\ \end{cases}\]

После этого вектор $\vec{x}$ приобретет следующий вид:

\[\vec{x}= \begin{pmatrix} A_x & A_y & \alpha_0 & \alpha_1 & D_x & D_y \end{pmatrix} ^\intercal\]

Тот же самый подход можно очевидным образом использовать и для геометрических ограничений типа HORIZONTALITY и VERTICALITY.

Теперь поговорим о том, как реализовать этот подход с точки зрения программирования. Конечно, можно придумать много подходов, я опишу мой. Массив Vars описывает набор переменных, входящих в вектор $\vec{x}$. Также заведем массив Links и заполним его значениями -1. Это число, -1, будет означать, что переменная является “базовой”, то есть она напрямую входит в вектор $\vec{x}$:

Теперь, после добавления геометрического ограничения совпадения точек $B$ и $C$, добавим $\alpha_0$ и $\alpha_1$ в конец массива Vars. Также, добавим два числа -1 в конец массива Links. И теперь основной момент: заменим значения массива Links, соответствующие переменным $B_x$ и $C_x$, на индекс переменной $\alpha_0$, то есть 8. Также, заменим значения массива Links, соответствующие переменным $B_y$ и $C_y$, на индекс переменной $\alpha_1$, то есть 9:

Дальше все просто. Только те переменные, которым соответствует значение -1 в массиве Links, являются “базовыми”, то есть напрямую входят в вектор $\vec{x}$. Если же какой-то переменной соответствует какое-то другое число в массиве Links, то с помощью описанной структуры можно легко найти соответствующую ей “базовую” переменную.

Имплементация геометрических ограничений

Тут все довольно очевидно, достаточно просто каждое геометрическое ограничение представить в виде функции $\vec{g}(\vec{x})$, которая равна нулю в случае, когда геометрическое ограничение удовлетворено. Некоторые геометрические ограничения разобраны в разделе Пример, имплементацию остальных можно посмотреть в коде.