Implementation of analytic projective geometry for computer graphics

Full Text

1. Введение

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

В проективной геометрии непропорционально большую роль играет визуализация геометрических структур. Авторам не удалось найти полной системы компьютерной алгебры, специализирующейся на методах проективной геометрии, что несомненно связано с ее спецификой. Ранее мы использовали методы проективной геометрии для обоснования применения гиперкомплексных чисел для описания движений в разных вариантах физических пространств [1–4]. В качестве следующего шага мы использовали подход геометрической алгебры (с использованием систем компьютерной алгебры) для исследования тех же объектов под другим углом зрения [5–7]. На следующем шаге мы решили вернуться к прямому применению теоретических основ к исследованиям. В качестве первого этапа данной научной программы предлагается частичная компьютеризация исследований. Алгебраические формулы реализованы в виде структур данных и наборов функций, связанных с геометрическими структурами. Затем с помощью языка Asymptote сделана программная реализация для визуализации полученных алгебраических структур. В данном случае к структурам относятся точки, прямые и плоскости.

Язык Asymptote имеет структуру системы компьютерной алгебры, но вместо символьного представления алгебраических объектов в результате работы Asymptote получается их визуальное представление.

Система Asymptote [8–12] — специализированный высокоуровневый язык программирования для создания векторных изображений научной направленности (геометрические чертежи, схемы, диаграммы и т.д.). Язык использует C-подобный синтаксис с некоторым числом дополнений и улучшений¹, что по замыслу разработчиков делает его более удобным для освоения и использованием по сравнению с MetaPost [14] и PostScript.

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

В данной статье мы не будем излагать основ языка Asymptote. Заинтересованный читатель может найти описание всех его возможностей в официальном руководстве на сайте [15], базовые примеры приведены в обучающих заметках [11], также множество примеров с комментариями на русском языке можно посмотреть в источниках [16–18]. Следует отметить, что большая часть примеров относится к двумерным изображениям. В статье акцент сделан на описании реализации посредством Asymptote структур данных, таких как проективные точки, прямые и плоскости для трехмерных графиков, приведены примеры исходного кода с пояснениями.

1.1. Структура статьи

В начале статьи приведена реализация структур данных, описывающих проективные точки (раздел 2), прямые (раздел 3) и плоскости (раздел 4). Изложение сопровождается примерами исходного кода с пояснением всех специфических конструкций, ключевых слов и функций языка Asymptote. В разделе 5 приводятся несколько примеров использования созданных структур, а также описываются некоторые технические тонкости создания именно трехмерных изображений, так как в официальном руководстве данный аспект освещен довольно скупо.

1.2. Обозначения и соглашения

Напомним основные используемые обозначения проективной геометрии. Формулы приведены в однородных координатах, обозначения взяты из монографии [19]:

• $\left\{v\right|p\times v\}$  — прямая, проходящая через точку P по направлению v;
• $\{p_2-p_1|p_1\times p_2\}$  — прямая, проходящая через две точки P₁ и P₂;
• $\left\{p\right|0\}$ — прямая, проходящая через начало координат и точку P;
• $\{w_1{p}_2-w_2{p}_1|p_1\times p_2\}$ — прямая, проходящая через две точки ${\overrightarrow{p}}_1=\left(p_1\right|w_1)$ и ${\overrightarrow{p}}_2=\left(p_2\right|w_2)$;
• $\{n_1\times n_2|d_1{n}_2-d_2{n}_1\}$ — прямая линия пересечения двух плоскостей $\left[n_1\right|d_1]$ и $\left[n_2\right|d_2]$;
• $(m\times n+dv|-(n,v))$ — точка пересечения плоскости $\left[n\right|d]$ и прямой $\left\{v\right|m\}$;
• $[m_1\times m_2|(v_2,m_1)]$ — точка пересечения двух прямых $\left\{v_1\right|m_1\}$ и $\left\{v_2\right|m_2\}$;
• $(d_1{n}_3\times n_2+d_2{n}_1\times n_3+d_3{n}_2\times n_1|(n_1,n_2,n_3))$ — точка пересечения трех плоскостей $\left[n_1\right|d_1]$, $\left[n_2\right|d_2]$ и $\left[n_3\right|d_3]$;
• $(v\times m|(v,v))$ — точка, ближайшая к началу координат на прямой $\left\{v\right|m\}$;
• $(-dn|(n,n))$ — точка, ближайшая к началу координат на плоскости $\left[n\right|d]$;
• $[v\times u|-(u,m)]$ — плоскость, содержащая прямую $\left\{v\right|m\}$ и направление $u$;
• $[v\times p+m|-(p,m)]$ — плоскость, содержащая прямую $\left\{v\right|m\}$ и точку $\left(p\right|1)$;
• $\left[m\right|0]$ — плоскость, содержащая прямую $\left\{v\right|m\}$ и начало координат;
• $[v\times p+wm|-(p,m)]$ — плоскость, содержащая прямую $\left\{v\right|m\}$ и точку $\left(p\right|w)$;
• $[v\times u|(u,v,p)]$ — плоскость, содержащая точку $\left(p\right|1)$ и направления v и u;
• $[m\times v|(m,m)]$ — плоскость с прямой $\left\{v\right|m\}$, наиболее отдаленная от O;
• $[-wp|(p,p)]$ — плоскость с точкой $\left(p\right|w)$, наиболее отдаленная от O;
• $\frac{(v_1,m_2)+(v_2,m_1)}{v_1\times v_2}$ — расстояние между прямыми $\left\{v_1\right|m_1\}$ и $\left\{v_2\right|m_2\}$;
• $\frac{v\times p+m}{v}$ — расстояние между прямой $\left\{v\right|m\}$ и точкой $\left(p\right|1)$;
• $m/v$ — расстояние от прямой $\left\{v\right|m\}$ до начала координат;
• $\frac{(n,p)+d}{n}$ — расстояние от плоскости $\left[n\right|d]$ до точки $\left(p\right|1)$;
• $d/n$ — расстояние от плоскости $\left[n\right|d]$ до начала координат.

2. Создание структур на примере однородных координат точки

Для создания пользовательских структур в языке Asymptote используется синтаксис, во многом похожий на соответствующий синтаксис языка Си, однако имеющий некоторую специфику, которую мы рассмотрим на примере структуры для однородных координат.

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

Язык Asymptote имеет встроенный тип triple, который используется для задания координат точек в трехмерном пространстве. Структура данных для однородных координат отсутствует, по крайней мере недоступна пользователям, поэтому следует ее создать.

Назовем структуру Vec4, по аналогии с соответствующим типом данных из GLSL (OpenGL Shading Language, Graphics Library Shader Language):

struct Vec4 {
real x; real y; real z; real w;
pair xy; triple xyz;
real[] xyzw;
}

Структура данных должна содержать всего 4 числа типа real, однако для удобства использования кроме полей x, y, z, w, соответствующих однородным координатам $x:y:z:w$, мы задаем поля, позволяющие получить только координаты xy в виде типа данных pair, трехмерные декартовы координаты xyz в виде типа данных triple и все координаты в виде массива xyzw. Кроме удобства при манипуляциях с указанными объектами такой подход также повторяет программный интерфейс языка GLSL.

После определения структуры Vec4 Asymptote автоматически создает одноименную функцию, обязательными аргументами которой являются все поля, перечисленные при определении структуры. Данная функция эквивалентна конструктору по умолчанию в терминах языка C++ или Java.

В нашем случае придется перечислить все семь полей, хотя они частично дублируют друг друга. Данная проблема решается заданием специализированных операторов инициализации (конструкторов), для чего следует перегрузить специальный оператор init. Внутри этого оператора можно использовать ключевое слово this, которое позволяет сослаться на экземпляр структуры.

Для нашего примера создадим три специализированных конструктора:

• Конструктор Vec4 из набора чисел:

void operator init(real x, real y,
$\hookrightarrow$ real z, real w) {
this.x = x;
this.y = y;
this.z = z;this.w = w;
this.xy = (x, y);
this.xyz = (x, y, z);
this.xyzw = new real[] {this.x,
$\hookrightarrow$ this.y, this.z, this.w};
}

• Конструктор Vec4 из декартовых координат точки и отдельного скаляра w:

void operator init(triple p, real w)
$\hookrightarrow$ {
this.x = p.x;
this.y = p.y;
this.z = p.z;
this.w = w;
this.xy = (p.x, p.y);
this.xyz = (p.x, p.y, p.z);
this.xyzw = new real[] {this.x,
$\hookrightarrow$ this.y, this.z, this.w};
}

• Конструктор Vec4 из массива:

void operator init(real[] v) {
this.x = v[0];
this.y = v[1];
this.z = v[2];
this.w = v[3];
this.xy = (v[0], v[1]);
this.xyz = (v[0], v[1], v[2]);
this.xyzw = v[0:5];
}

Все перечисленные операторы следует размещать внутри пространства имен структуры (т.е. между открывающей и закрывающей фигурными скобками структуры Vec4).

Для удобства манипулирования данными можно также перегрузить оператор cast, который ответственен за неявное преобразование типов данных:

• Преобразование типа Vec4 в массив:

real[] operator cast(Vec4 v) {
return v.xyzw;
}

• Преобразование типа Vec4 в трехмерные декартовы координаты:

triple operator cast(Vec4 v) {
return v.xyz / v.w;
}

• Преобразование массива в тип Vec4:

Vec4 operator cast(real[] v) {
return Vec4(v);
}

При преобразовании данных в тип triple осуществляется дополнительное деление на весовую координату w. В случае несобственной (идеальной) точки данное преобразование сработает некорректно, так как у идеальной точки координата w = 0 и произойдет деление на ноль. Однако проверку на равенство нулю координаты w всегда можно произвести отдельно.

3. Проективные прямые в виде структуры Line3D

3.1. Структура и конструкторы

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

struct Line3D {
// Направляющий вектор и момент:
triple V, m;
// Единичный направляющий вектор:
triple v;
// Произвольная точка прямой:
triple P;
}

Для полноценного использования созданной структуры в ее области имен необходимо задать несколько дополнительных операторов инициализации. Так, прямая, проходящая через две собственные точки, заданные в декартовых координатах, вычисляется по формуле:

$m=p_1\times p_2.$ (1)

Соответствующий конструктор будет выглядеть следующим образом:

void operator init (triple P2, triple
$\hookrightarrow$ P1) {
this.V = P2 - P1;
this.v = unit(this.V);
this.m = cross(P1, P2) /
$\hookrightarrow$ length(this.V);
this.P = P1;
}

Следует обратить внимание, что для вычисления единичного вектора v мы использовали встроенную в Asymptote функцию unit, для вычисления векторного произведения — функцию cross, а для вычисления нормы вектора — функцию length. Обратите внимание, что при вычислении m используется функция (1) в комбинации с соотношением 

$$\begin{gathered} d_{O{O}_{\perp }}=O{O}_{\perp }= \\ =\frac{v\times m}{v}=\frac{vm\sin \frac{\pi }{2}}{v^2}=\frac{m}{v}, \end{gathered}$$

т.е. она делится на исходную длину вектора v, вычисленного как p₂ – p₁. В этом случае длина m имеет геометрический смысл расстояния от начала координат до прямой.

Следующий конструктор инициализирует структуру по заданным векторам v и m:

void operator init (triple keyword V,
$\hookrightarrow$ triple keyword m) {
this.V = V;
this.v = unit(this.V);
this.m = m / length(this.V);
this.P = cross(this.v, this.m);
}

Вектор v может быть произвольной длины, его исходное значение записывается в поле V структуры, а его нормированное значение — в поле v. Здесь вектор m нормализован путем деления на собственную длину ненормированного касательного вектора ||v||.

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

triple V = (1,1,1); triple m = (1, 1,
$\hookrightarrow$ 1);
Line3D L = Line3D(V=V, m=m);

Если осуществить вызов через Line3D L = Line(V, m), то будет вызван предыдущий конструктор по двум точкам, заданным в декартовых координатах.

В реализации конструктора, который инициализирует прямую, проходящую через две точки, заданные в проективных координатах, используется следующая формула:

$\{w_1{p}_2-w_2{p}_1|p_1\times p_2\}.$ (2)

Соответствующий конструктор будет выглядеть следующим образом:

void operator init (Vec4 P2, Vec4 P1) {
this.V = P1.w * P2.xyz - P2.w *
$\hookrightarrow$ P1.xyz;
this.v = this.V / length(this.V);
this.m = cross(P1.xyz, P2.xyz) /
$\hookrightarrow$ length(this.V);
this.P = P1.xyz / P1.w;
}

Мы неявно подразумеваем, что по крайней мере P₁ — собственная точка, т.е. ее координата w не равна нулю. Подчеркнем, что формула (3) — универсальна и работает для всех возможных вариантов точек. Например, если ${\overrightarrow{p}}_1=(x,y,z,w)$ — собственная точка, а ${\overrightarrow{p}}_2=\left(v\right|0)=(v_x,v_y,v_z\mathrm{,0})$ — несобственная точка, то координаты Плюккера прямой вычисляются как $\overrightarrow{l}=\left\{wv\right|p\times v\}.$

Произведя нормализацию путем деления на w|v|, получим формулу:

$\left\{v\right|m\}=\{v_x,v_y,v_z|m_x,m_y,m_z\}=\left\{v\right|p\times v\}.$ (3)

Формуле (3) соответствует следующий конструктор (прямая, проходящая через точку по направляющему вектору):

void operator init (Vec4 P, triple V) {
this.V = V;
this.v = this.V / length(this.V);
this.m = cross(P.xyz, this.V) /
$\hookrightarrow$ length(this.V);
this.P = P.xyz;
}

Если же обе точки ${\overrightarrow{p}}_1=\left(v_1\right|0)$ и ${\overrightarrow{p}}_2=\left(v_2\right|0)$ являются несобственными, то получаем несобственную прямую $\left\{0\right|v_1\times v_2\}$. Несобственные прямые могут встречаться при рассмотрении взаимного расположения плоскостей. Если плоскости параллельны, то они будут пересекаться как раз по несобственным прямым. Пример формулы (3) иллюстрирует преимущество проективного подхода, который позволяет одной формулой охватить сразу все возможные варианты.

Осталось создать еще два конструктора, которые реализуют формулу $\left\{v\right|0\}.$

Соответствующие конструкторы имеют вид:

• Прямая, проходящая через начало координат и точку P (в декартовых координатах):

void operator init (Vec4 P) {
this.V = P.xyz;
this.v = unit(this.V);
this.m = (0, 0, 0);
this.P = P.xyz;
}

• Прямая, проходящая через начало координат и точку P:

void operator init (triple P) {
this.V = P;
this.v = unit(this.V);
this.m = (0, 0, 0);
this.P = P;
}

Для визуализации прямой необходимо знать хотя бы две ее точки. Для реализации параметрического уравнения прямой определим в пространстве имен структуры Line3D следующую функцию:

triple get_point(real t) {
return this.P + this.v * t;
}

Удобно также иметь возможность вычислить сразу массив точек прямой, что крайне просто делается с помощью поэлементного цикла for:

triple[] get_points(real[] T) {
triple[] res;
for (real t : T) {
res.push(get_point(t));
}
return res;
}

Данный вариант цикла for похож на аналогичный цикл из языка Python, позволяет перебирать элементы из массива и проводить с ними манипуляции.

3.2. Взаимное расположение прямых

В предыдущем пункте мы задали поля структуры, набор конструкторов и несколько методов. Теперь зададим несколько функций, аргументами которых будут экземпляры структуры Line3D. Данные функции разместим в том же файле, что и структура Line3D, но вне ее пространства имен (за закрывающей фигурной скобкой).

В первую очередь реализуем функцию вычисления взаимного момента двух прямых l₁ и l₂ по простой формуле $M=(v_1,m_2)+(v_2,m_1).$

Реализация соответствующей функции имеет вид:

real reciprocal_moment(Line3D line01,
$\hookrightarrow$ Line3D line02) {
return dot(line01.m, line02.v) +
$\hookrightarrow$ dot(line01.v, line02.m);
}

Здесь мы использовали встроенную функцию dot, которая вычисляет скалярное произведение двух векторов.

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

• если M = 0, то прямые лежат в одной плоскости;
• если M > 0, то переход от l₁ к l₂ осуществляется правым поворотом;
• если M < 0, то переход от l₁ к l₂ получается левым поворотом.

Заметим, что делать проверку на равенство нулю следует с учетом погрешности представления вещественных чисел числами с плавающей точкой.

Далее определим функцию, реализующую формулу

$d=d_2-d_1=\frac{(v_1,m_2)+(v_2,m_1)}{v_1\times v_2}.$ (4)

Соотношение (4) позволяет вычислить длину d взаимного перпендикуляра между двумя скрещивающимися прямыми. Так как у нас уже есть функция вычисления взаимного момента, то вычисления длины перпендикуляра будет осуществляться также в одну строчку:

real reciprocal_perp_length(Line3D line01,
$\hookrightarrow$ Line3D line02) {
return reciprocal_moment(line01, line02)
$\hookrightarrow$ / length(cross(line01.v, line02.v));
}

Предыдущую функцию можно было бы сделать более общей, добавив проверку скрещиваемости прямых, и в случае, если прямые лежат в одной плоскости, применить формулу

$d_{Q{Q}_{\perp }}=\frac{||m+v\times q||}{||v||}=\frac{||(p-q)\times v||}{||v||}$

в виде

$d=\frac{v\times p_2+m_1}{v_1},$

где  ${\overrightarrow{l}}_1=\left\{v_1\right|m_1\}$ — координаты Плюккера первой прямой, ${\overrightarrow{p}}_2=\left(p_2\right|1)$ — однородные координаты произвольной точки второй прямой. Однако мы не стали усложнять функцию дополнительными внутренними проверками. Также реализуем формулу

$[m_1\times m_2|(v_2,m_1)]=[m_2\times m_1\left|\right(v_1,m_2\left)\right].$ (5)

При этом сделаем проверку на равенство нулю взаимного момента, так как в этом случае две прямые будут лежать в одной плоскости (будут компланарны):

// line-line-meet
Vec4 intersection(Line3D line01, Line3D
$\hookrightarrow$ line02) {
real M = reciprocal_moment(line01,
$\hookrightarrow$ line02);
assert( abs(M) < 1e-8, "Прямые не
$\hookrightarrow$ компланарны" );
return Vec4( cross(line01.m, line02.m),
$\hookrightarrow$ dot(line02.v, line01.m));
// return Vec4( cross(line02.m,
$\hookrightarrow$ line01.m), dot(line01.v, line02.m));
}

Как указывалось выше, проверку на равенство нулю следует сделать с некоторой погрешностью. В нашей реализации мы допускаем довольно большую погрешность, так как основной целью является визуализация полученных объектов, где малые погрешности не критичны.

Для вычисления координат Плюккера прямой, содержащей взаимный перпендикуляр, воспользуемся довольно-таки громоздкой формулой:

$$\begin{gathered} v_{\perp }=v_1\times v_2, \\ {m}_{\perp }=m_1\times l_2-m_2\times v_1+ \\ +\frac{(v_1,m_2)+(v_2,m_1)}{v_1\times v_2}(v_1,v_2)v_1\times v_2. \end{gathered}$$

Компьютерная реализация имеет вид:

Line3D reciprocal_perp_line(Line3D line01,
$\hookrightarrow$ Line3D line02) {
triple V = cross(line01.v, line02.v);
triple m = cross(line01.m, line02.v) -
$\hookrightarrow$ cross(line02.m, line01.v) +
$\hookrightarrow$reciprocal_perp_length(line01,
$\hookrightarrow$line02) * dot(line01.v, line02.v) *
$\hookrightarrow$cross(line01.v, line02.v) /
$\hookrightarrow$length(cross(line01.v, line02.v));
return Line3D(V=V, m=m);
}

4. Проективные плоскости в виде структуры Plane3D

4.1. Структура и конструкторы

Теперь рассмотрим структуру, реализующую задание плоскости в однородном виде. Для полного описания плоскости достаточно вектора $\overrightarrow{\pi }=\left[n\right|d],$ где n — единичный вектор нормали к плоскости, а d — ориентированное расстояние от начала координат до плоскости.

В модуле three языка Asymptote существует функция plane со следующей сигнатурой:

path3 plane(triple u, triple v,
$\hookrightarrow$ tripleO=O);

Здесь u и v — направляющие векторы плоскости, а 0 — точка их приложения. Данная функция возвращает замкнутый прямоугольный контур 0 -- 0+u -- 0+u+v -- 0+v -- cycle, который затем можно отобразить в виде плоскости с помощью функции surface. Поэтому в структуре Plane3D полезно хранить два направляющих вектора v и u, а также произвольную точку P. Полностью поля структуры будут выглядеть следующим образом:

struct Plane3D {
// Вектор нормали (ненормированный):
triple N;
// Единичный вектор нормали:
triple n;
// Направляющие векторы (касательные
$\hookrightarrow$ векторы):
triple U, V;
// Единичные направляющие векторы:

triple u, v;
// Ориентированное расстояние от
$\hookrightarrow$ начала координат до прямой:
real d;
// Произвольная точка плоскости:
triple P;
}

Зададим несколько сравнительно тривиальных операторов инициализации. Первый из них создает плоскость, содержащую прямую $\overrightarrow{l}=\left\{v\right|n\}$ и направляющий вектор v:

void operator init (Line3D line, triple
$\hookrightarrow$ U) {
this.U = U;
this.V = line.V;
this.u = unit(U);
this.v = line.v;
this.N = cross(this.v, this.u);
this.n = unit(this.N);
this.d = -dot(this.u, line.m);
this.P = -this.d * this.n;
}

В качестве точки плоскости берется ближайшая к началу координат. Вычисление проводится по формуле $[v\times u|-(u,m)].$

Второй конструктор инициализирует плоскость, проходящую через точку и прямую. Точка может быть задана однородными координатами, тогда работает формула $[v\times p+wm|-(p,m)].$

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

$[v\times p+m|-(p,m)].$

Так как случай однородных координат более общий, то мы ограничимся только им в реализации конструктора: 

void operator init (Line3D line, Vec4 P)
$\hookrightarrow$ {
this.U = line.P - P.xyz;
this.V = line.V;

this.u = unit(U);
this.v = line.v;
this.N = cross(line.v, P.xyz) + P.w *
$\hookrightarrow$ line.m;
this.n = unit(this.N);
this.d = -dot(P.xyz, line.m);
this.P = P.xyz;
}

Третий конструктор инициализирует плоскость, проходящую через прямую и начало координат:

void operator init (Line3D line) {
this.U = line.P;
this.V = line.V;
this.u = unit(U);
this.v = line.v;
this.N = line.m;
this.n = unit(this.N);
this.d = 0;
this.P = line.P;
}

Вычисления проводятся по следующей формуле: $\left[m\right|0].$

Четвертый конструктор создает плоскость, проходящую через заданную точку и два направляющих вектора. Это эквивалентно записи параметрического уравнения плоскости. Вычисления проводятся по формуле (8). В Asymptote нет смешанного произведения, но оно естественно заменяется на комбинацию скалярного и векторного произведений:

void operator init (triple P, triple U,
$\hookrightarrow$ triple V) {
this.U = U;
this.V = V;
this.u = unit(U);
this.v = unit(V);
this.N = cross(U, V);
this.n = unit(this.N);
this.d = dot(V, cross(U, P)) /
$\hookrightarrow$ length(this.N);
this.P = P;
}

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

void operator init (real A, real B, real
$\hookrightarrow$ C, real D) {
this.N = (A, B, C);
this.n = unit(this.N);
triple w;
if ( A <= B && A <= C) {
w = this.n + (1, 0, 0);
} else if (B <= A && B <= C) {
w = this.n + (0, 1, 0);
} else {
w = this.n + (0, 0, 1);
}
this.u = unit(cross(w, this.n));
this.v = cross(this.n, this.u);
this.U = u;
this.V = v;
this.d = D / length(this.N);
this.P = -this.d * this.n;
}

По известному единичному нормальному вектору плоскости n можно вычислить направляющие векторы плоскости u и v. Для этого можно использовать простую процедуру, указанную в [20, с. 29]. Возьмем произвольный вектор w, не коллинеарный n. Чтобы найти такой вектор, можно взять наименьшую компоненту вектора n и увеличить ее на 1, например $w=(n_x,n_y,n_z+1)$, если n_z < n_x и n_z < n_y. После этого можно вычислить вектор u: 

$u=\frac{w\times n}{||w\times n||}.$

Из свойств векторного произведения следует, что . При вычислениях с помощью компьютера полезно также удостовериться, что значение ||w × n|| не слишком велико или мало, чтобы избежать лишних погрешностей, связанных с машинной точностью.

На последнем шаге находим v как v = n × u. Так как ||n|| = ||u|| = 1, то и v — единичный вектор. Таким образом мы получили ортонормированную тройку векторов $⟨n,u,v⟩$.

Точку плоскости можно найти как проекцию начала координат на плоскость $O{O}_{\perp }=-dn$ при условии, что n — единичный вектор: $(-dn|{||n||}^2).$

В области имен структуры зададим еще две функции. Первая из них фактически является параметрическим уравнением плоскости: 

triple get_point(real s, real t) {
return this.P + this.u * s + this.v *
$\hookrightarrow$ t;
}

в виде прямоугольника:

guide3 get_contour() {
triple P = this.P - 0.5(this.u +
$\hookrightarrow$ this.v);
return P -- P + this.u -- P + this.u
$\hookrightarrow$ + this.v -- P +this.v -- cycle;
}

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

4.2. Взаимное расположение плоскостей

Имея в наличии структуру Plane3D, можно перейти к решению задачи определения взаимного положения двух и трех плоскостей.

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

$\left\{v\right|m\}=\{n_1\times n_2|d_1{n}_2-d_2{n}_1\}.$ (6)

Возможны следующие варианты:

• если v ≠ 0, то прямая собственная;
• если v = 0, то прямая несобственная;

- если d₁ = d₂, то плоскости совпадают;

- если d₁ ≠ d₂, то плоскости параллельны.

Отметим также, что если две плоскости совпадают, то необязательно n₁ = n₂, так как возможен вариант, когда n₁ = –n₂. В этом случае лицевая и внутренняя стороны плоскостей отличаются.

В случае трех плоскостей вариантов может быть больше. Все они изображены на рис. 1 и 2.

 

Рис. 1. Варианты непараллельных плоскостей.

 

Рис. 2. Варианты параллельного расположения плоскостей. Вариант (в) — предельный случай варианта (б).

 

Чтобы выписать формулы для вычислений, введем следующие обозначения. Плоскости π₁, π₂ и π₃ запишем в координатах Плюккера в следующем виде: $\overrightarrow{\pi }1 = [n1 | d1], \overrightarrow{\pi }2 = [n2 | d2] и \overrightarrow{\pi }3 = [n3 | d3].$ Каждая пара этих трех плоскостей пересекается по некоторой прямой l_ij, где i, j = 1, 2, 3 и i < j.

Тогда по формуле (14) запишем

${\overrightarrow{l}}_{ij}=\{n_i\times n_j|d_i{n}_j-d_j{n}_i\}.$

То есть получается следующее соотношение: $v_{ij}=n_i\times n_j$ и $m_{ij}=d_1{n}_2-d_2{n}_1$.

Все возможные варианты расположения трех плоскостей можно определить по следующей схеме:

• точка пересечения является собственной, т.е. $(n_1,n_2,n_3)\cancel{=}0;$
• точка пересечения является несобственной, т.е. $(n_1,n_2,n_3)=0$:

- все прямые являются собственными, параллельными и различными (не совпадают), т.е. $v_{12}\times v_{13}=v_{13}\times v_{23}=v_{23}\times v_{12}=0$ и m₁₂ ≠ m₂₃ ≠ m₁₃;

- все прямые являются собственными, параллельными и совпадают, т.е. $v_{12}\times v_{13}=v_{13}\times v_{23}=v_{23}\times v_{12}=0$ и m₁₂ = m₂₃ = m₁₃;

- две прямые являются собственными, параллельными и различными (не совпадают), т.е. v₁₃ × v₂₃ = 0 и m₁₃ ≠ m₂₃, а одна прямая является несобственной, т.е. v₁₂ = 0;

- все прямые являются несобственными, т.е. $v_{12}=v_{23}=v_{13}=\mathrm{0,}$ а плоскости не совпадают, т.е. d₁ ≠ d₂ ≠ d₃;

- все прямые являются несобственными, т.е. $v_{12}=v_{23}=v_{13}=\mathrm{0,}$ а плоскости совпадают, т.е. d₁ = d₂ = d₃.

Однородные координаты точки пересечения трех плоскостей вычисляются по формуле:

$\begin{array}{c}p=\left(\begin{array}{l}{d}_1{n}_3\times n_2+d_2{n}_1\times n_3+\\ +d_3{n}_2\times n_1|\left(n_1,n_2,n_3\right)\end{array}\right)=\\ =\left({\Delta }_x,{\Delta }_y,{\Delta }_z|\Delta \right).\end{array}$ (15)

Обозначения Δ_x, Δ_y, Δ_z представляют собой определители решения системы уравнений  методом Краммера:

$\left\{\begin{array}{c}{n}_{1x}x+n_{1y}y+n_{1z}z=-d_1,\\ n_{2x}x+n_{2y}y+n_{2z}z=-d_2,\\ n_{3x}x+n_{3y}y+n_{3z}z=-d_3.\end{array}\right.$

$\begin{array}{lll}\Delta & =\left|\begin{array}{ccc}{n}_{1x}& n_{1y}& n_{1z}\\ n_{2x}& n_{2y}& n_{2z}\\ n_{3x}& n_{3y}& n_{3z}\end{array}\right|,{\Delta }_x& =\left|\begin{array}{ccc}-d_1& n_{1y}& n_{1z}\\ -d_2& n_{2y}& n_{2z}\\ -d_3& n_{3y}& n_{3z}\end{array}\right|,\\ {\Delta }_y& =\left|\begin{array}{ccc}{n}_{1x}& -d_1& n_{1z}\\ n_{2x}& -d_2& n_{2z}\\ n_{3x}& -d_3& n_{3z}\end{array}\right|,{\Delta }_z& =\left|\begin{array}{ccc}{n}_{1x}& n_{1y}& -d_1\\ n_{2x}& n_{2y}& -d_2\\ n_{3x}& n_{3y}& -d_3\end{array}\right|.\end{array}$

Обозначение (Δ_x, Δ_y, Δ_z | Δ) задает точку в однородных координатах:

$\overrightarrow{p}=\left(\frac{{\Delta }_x}{\Delta },\frac{{\Delta }_y}{\Delta },\left.\frac{{\Delta }_z}{\Delta }\right|1\right)=:\left({\Delta }_x,{\Delta }_y,{\Delta }_z|\Delta \right).$

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

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

Line3D intersection(Plane3D plane01,
$\hookrightarrow$ Plane3D plane02) {
Line3D line;
line.m = plane01.d * plane02.n -
$\hookrightarrow$ plane02.d * plane01.n;
line.V = cross(plane01.n, plane02.n);
// Если плоскости пересекаются по
$\hookrightarrow$ несобственной прямой, то v = 0
if (abs(line.V) < 1e-8) {
line.v = (0, 0, 0);
} else {
line.v = unit(line.V);
// проведем нормировку, разделив на
$\hookrightarrow$ длину направляющего вектора
line.m = line.m / length(line.V);
line.P = cross(line.v, line.m) /
$\hookrightarrow$ dot(line.v, line.v);
}
return line;
}

Необходимые проверки были встроены непосредственно в код. Функция abs вычисляет норму вектора (можно заменить на функцию length). Еще раз следует заметить, что в случае использования чисел с плавающей запятой необходимо учитывать машинную погрешность и использовать не строгое равенство, а проверять, что модуль числа не превосходит некоторое малое число.

Код для вычисления точки пересечения трех плоскостей будет выглядеть следующим образом:

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

Еще одна функция реализует формулу для вычисления однородных координат точки пересечения плоскости $\left[n\right|d]$ и прямой $\left\{v\right|m\}$:

$\begin{array}{c}\left(-\left.\frac{\left(m\times n+dv\right)}{\left(n,v\right)}\right|1\right)=\left(-\left(m\times n+dv\right)|\left(n,v\right)\right)=\\ =\left(\left(m\times n+dv\right)|-\left(n,v\right)\right).\end{array}$ (16)

Здесь также можно поступить двояко: вернуть тип triple и встроить все проверки внутрь функции, или же вернуть тип Vec4 и возложить ответственность за проверки на пользователя функции. Здесь мы приводим первый вариант функции:

triple intersection(Plane3D plane, Line3D
$\hookrightarrow$ line) {
assert( dot(line.v, plane.n) > 1e-8,
$\hookrightarrow$"Плоскость и прямая не
$\hookrightarrow$пересекаются!" );
return -(cross(line.m, plane.n) +
$\hookrightarrow$plane.d * line.v) / dot(line.v,
$\hookrightarrow$plane.n);
}

5. Пример визуализации

В данном разделе рассмотрим пример использования созданных структур и встроенных в Asymptote средств для визуализации комплексного изображения пересечения плоскостей и прямых. Приведем полный код программы и прокомментируем все его части.

Программа рисует три плоскости, которые пересекаются в одной точке. Изображаются прямые пересечения каждой пары плоскостей. Точка пересечения плоскостей одновременно является точкой пересечения прямых, по которым пересекаются данные плоскости. Результат показан на рис. 3.

 

Рис. 3. Пересечение трех плоскостей.

 

Третья плоскость специально выбрана так, чтобы ее нормальный вектор совпадал с направляющим вектором прямой пересечения первой и второй плоскостей, т.е. $n_3=n_1\times n_2=v_{12}$. Это изображение может служить иллюстрацией к доказательству формулы нахождения прямой, по которой пересекаются две плоскости.

Файл с исходным кодом Asymptote имеет расширение .asy. Для его запуска будем использовать команду asy. Вначале следует импортировать все используемые модули как встроенные, так и реализованные нами.

import settings;
settings.outformat = "pdf" ;
settings.render = 15;

include "Vec4.asy" ;
include "Line3D.asy" ;
include "Plane3D.asy" ;
import three;
import graph3;
unitsize(3cm);

Модуль settings позволяет настроить некоторые параметры работы Asymptote. В частности мы указываем, что следует создавать изображение в формате pdf, а также выставляем качество созданного трехмерного изображения. Значения больше нуля приведут к тому, что изображение получится растровым. В случае трехмерных чертежей мы вынуждены использовать растровый формат, так как иначе некорректно рассчитываются наложения поверхностей друг на друга. Качество растра следует подбирать опытным путем.

Далее импортируем три созданные нами структуры модуль three для поддержки трехмерных векторов и модуль graph3 для визуализации трехмерных поверхностей и кривых. Также выставляем единицу измерения по умолчанию в 3 сантиметра — такое значение было подобрано опытным путем исходя из критерия читаемости изображения при его публикации в статье или презентации.

В модуле three определена функция plane, которая принимает в качестве аргументов два направляющих вектора плоскости и точку их приложения, а в качестве результата возвращает контур плоскости в виде прямоугольника. Точка приложения векторов будет служить правым верхним углом данного прямоугольника. Далее этот контур можно использовать для отображения плоскости с помощью функции surface.

В первую очередь задаем направляющие векторы:

triple u_01 = X;
triple v_01 = Y;
triple u_02 = rotate(30, Y) * u_01;
triple v_02 = v_01;
triple u_03 = u_01;
triple v_03 = rotate(-90, X) * v_01;

Для первой плоскости в качестве направляющих выберем орты осей Ox и Oy, используя заданные в модуле graph3 векторы X, Y.

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

// Произвольная точка, к которой будем
$\hookrightarrow$ привязывать направляющие векторы
triple P_0 = (0, -0.5, 1);
// Точка плоскости, от которой будут
$\hookrightarrow$ откладываться направляющие векторы при
$\hookrightarrow$рисовании
triple P_01 = P_0 - u_01;
triple P_02 = P_0 - u_02;
triple P_03 = P_0 - X + 0.5Y + 0.5Z;

Далее можно создать контуры и поверхности плоскостей:

// Создаем контуры плоскостей
path3 pl_01 = plane(u=2u_01, v=v_01,
$\hookrightarrow$ O=P_01);
path3 pl_02 = plane(u=2u_02, v=v_02,
$\hookrightarrow$ O=P_02);
path3 pl_03 = plane(u=2u_03, v=v_03,
$\hookrightarrow$ O=P_03);
surface pls_01 = surface(pl_01);
surface pls_02 = surface(pl_02);
surface pls_03 = surface(pl_03);

Чтобы созданные объекты были отрисованы, необходимо воспользоваться функцией draw. Покажем как отрисовать первую плоскость:

draw(
s=pls_01,
nu=1, nv=1,
surfacepen=lightgrey+opacity(.8),
meshpen=black+thin(),
light=nolight
);

Разберем передаваемые в данную функцию параметры:

• в параметр s передается объект surface;
• параметры nu и nv задают количество точек сетки сплайновой поверхности;
• параметр surfacepen позволяет настроить перо для рисования поверхности;
• параметр meshpen отдельно настраивает перо для сетки, в нашем случае от сетки остается только контур плоскости;
• с помощью параметра light настраивается свет.

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

draw(pl_01, p=0.25bp+black,
$\hookrightarrow$ L=Label("$\pi_1$" ,
$\hookrightarrow$position=Relative(0.80)));
draw(pl_02, p=0.25bp+black,

$\hookrightarrow$L=Label("$\pi_2$" ));
draw(pl_03, p=0.25bp+black,
$\hookrightarrow$L=Label("$\pi_3$" ,
$\hookrightarrow$position=Relative(0.3)));

Здесь параметр p задает перо, параметр L — подпись к изображаемому объекту. Для создания подписей можно использовать нотацию LaTeX, что делает Asymptote удобным для создания иллюстраций к научным текстам.

Далее перейдем к использованию созданных нами структур. В первую очередь зададим все три плоскости:

Plane3D plane01 = Plane3D(P_01, u_01,
$\hookrightarrow$ v_01);
Plane3D plane02 = Plane3D(P_02, u_02,
$\hookrightarrow$ v_02);
Plane3D plane03 = Plane3D(P_03, u_03,
$\hookrightarrow$ v_03);

Затем найдем прямые пересечения данных плоскостей:

Line3D line12 = intersection(plane01,
$\hookrightarrow$ plane02);
Line3D line13 = intersection(plane01,
$\hookrightarrow$ plane03);
Line3D line23 = intersection(plane02,
$\hookrightarrow$ plane03);

Прямые отображаются в виде отрезков, поэтому следует знать две точки каждой прямой. Подбираются точки эмпирически, а чтобы получить координаты точек используем метод get_point. Например, для первой прямой:

triple l12_sp = line12.get_point(-1.6);
triple l12_ep = line12.get_point(+0.9);

Далее отображаем прямые. Так, для первой прямой:

draw(l12_sp--l12_ep, p=black,
$\hookrightarrow$L=Label("$l_{12}$" ,
$\hookrightarrow$position=BeginPoint));

Теперь находим точку пересечения и отображаем ее:

triple intersection_point =
$\hookrightarrow$intersection(plane01, plane02,
$\hookrightarrow$plane03);
dot(intersection_point);

Наконец, откладываем нормальные векторы всех трех плоскостей от точки их пересечения:

draw(
intersection_point -- intersection_point
$\hookrightarrow$ + plane01.n,
arrow=Arrow3(size=4),
p=black,
L=Label(s="$\vb{n}_1$" , align=E,
$\hookrightarrow$ position=Relative(0.9))
);

6. Заключение

Мы рассмотрели реализацию аналитической проективной геометрии на языке Asymptote. Были реализованы однородные координаты, координаты Плюккера прямой и однородные координаты плоскости. Подробно описаны три созданные нами структуры, инициализирующие операторы и функции. Фундаментальным преимуществом по сравнению с классической аналитической геометрией является существенное упрощение вычислений, так как большинство функций просто повторяют проективные формулы и вычислительная часть зачастую умещается в несколько строк.

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

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

Второе направление — техническое. Оно заключается в том, чтобы реализовать набор функций для визуализации точек, прямых и плоскостей, заданных непосредственно в проективном виде.

Источники финансирования

Публикация выполнена в рамках проекта № 021934-0-000 Системы грантовой поддержки научных проектов РУДН (Геворкян М.Н., Королькова А.В., Севастьянов Л.А.) и при поддержке Программы стратегического академического лидерства РУДН (Кулябов Д.С.).

 

¹ Сам язык в своем синтаксисе скорее ориентируется на C++ [13]. Например, переменные в Asymptote можно объявлять в любом месте программы и инициализировать сразу при объявлении, а функции могут принимать параметры по умолчанию.

About the authors

M. N. Gevorkyan

RUDN University

Author for correspondence.
Email: gevorkyan-mn@rudn.ru
Russian Federation, 6 Miklukho-Maklaya St, Moscow, 117198

A. V. Korol’kova

Joint Institute for Nuclear Research

Email: korolkova-av@rudn.ru
Russian Federation, 6 ul. Zholio-Kyuri 6, Dubna, Moscow oblast, 141980

D. S. Kulyabov

RUDN University; Joint Institute for Nuclear Research

Email: kulyabov-ds@rudn.ru
Russian Federation, 6 Miklukho-Maklaya St, Moscow, 117198; 6 ul. Zholio-Kyuri 6, Dubna, Moscow oblast, 141980

L. A. Sevast’yanov

RUDN University; Joint Institute for Nuclear Research

Email: sevastianov-la@rudn.ru
Russian Federation, 6 Miklukho-Maklaya St, Moscow, 117198; 6 ul. Zholio-Kyuri 6, Dubna, Moscow oblast, 141980

References

Supplementary files