[WIP] 2D Geometric Constraint Solver

Содержание

  1. Вступление
  2. Постановка задачи
  3. Пример
  4. Оптимизация
  5. Параметризация геометрических примитивов
  6. Обозначения

Вступление

Я провожу кучу времени в 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: концентричность дуг

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

Пример

Допустим, у нас есть два отрезка, $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}\]

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

Оптимизация

Я решил не имплементровать алгоритмы оптимизации самостоятельно, а воспользоваться готовым пакетом 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\]

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

Обозначения

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