15.01.2019 16:13
AkhramovichSA 26 10400 Источник
Если вы открыли данную статью, то наверняка уже слышали о кватернионах, и возможно даже используете их в своих разработках. Но пора подняться на уровень выше — к бикватернионам.

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

Определение бикватерниона


Бикватернион — гиперкомплексное число, имеющее размерность 8. В англоязычных статьях и литературе они называются — «dual quaternion», в русскоязычной литературе встречается еще названия «дуальный кватернион» или «комплексный кватернион».

Основное отличие от кватернионов заключается в том, что кватернион описывает ориентацию объекта в пространстве, а бикватернион еще и положение объекта в пространстве.

Бикватернион можно представить в виде двух кватернионов:

$\widetilde{\textbf{q}} = \begin{bmatrix} \textbf{q}_1 \\ \textbf{q}_2 \end{bmatrix},$


$\textbf{q}_1$ — действительная часть, определяет ориентацию объекта в пространстве;
$\textbf{q}_2$ — дуальная часть, определяет положение объекта в пространстве.

Бикватернион также называют еще комплексным кватернионом, в этом случае его представляют в виде кватерниона, каждый компонент которого представляет собой дуальное число (не путать с комплексным). Дуальное число $A = a_1 + \epsilon a_2$, где $a_1$ и $a_2$ — действительные числа, а $\epsilon$ — символ (комплексность) Клиффорда, обладающий свойством $\epsilon^2 = 0$. Не будем углубляться в математику, так как нас интересует больше прикладная часть, поэтому далее — бикватернион будем рассматривать как два кватерниона.

Геометрическая интерпретация бикватерниона


По аналогии с кватернионом, с помощью которого можно задать ориентацию объекта, бикватернионом можно задать еще и положение. Т.е. бикватернион задает сразу две величины — положение и ориентацию объекта в пространстве. Если их рассматривать в динамике, то бикватернион определяет две величины — линейную скорость перемещения и угловую скорость вращения объекта. На рисунке ниже показан геометрический смысл бикватерниона.



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

Скалярные характеристики


Рассмотрим основные скалярные характеристики. Здесь надо обратить внимание на то, что они возвращают не обычные вещественные числа, а дуальные.

1. Норма бикватерниона

$\|\widetilde{\textbf{q}}\| = \|\textbf{q}_1\| + \epsilon(q_{1_0} q_{2_0} + \textbf{q}_1^T \textbf{q}_2)$



2. Модуль бикватерниона

$|\widetilde{\textbf{q}}| = |\textbf{q}_1| + \epsilon\frac{q_{1_0} q_{2_0} + \textbf{q}_1^T \textbf{q}_2}{|\textbf{q}_1|}$



Основные операции


Рассмотрим основные операции работы с бикватернионами. Как вы можете заметить, они очень похожи на аналогичные операции работы с кватернионами.

1. Бикватернионное сопряжение

$\widetilde{\textbf{q}}^* = \begin{bmatrix} \textbf{q}_1^* \\ \textbf{q}_2^* \end{bmatrix}$



2. Бикватернионное сложение и вычитание

$\widetilde{\textbf{q}} \pm \widetilde{\textbf{p}} = \begin{bmatrix} \textbf{q}_1 \pm \textbf{p}_1 \\ \textbf{q}_2 \pm \textbf{p}_2 \end{bmatrix}$



Сложение и вычитание бикватернионов коммутативно (слагаемые можно менять местами).

3. Умножение действительного числа на бикватернион

$a\widetilde{\textbf{q}} = \widetilde{\textbf{q}}a = \begin{bmatrix} a\textbf{q}_1 \\ a\textbf{q}_2 \end{bmatrix}$



4. Бикватернионное умножение

$\widetilde{\textbf{q}} \otimes \widetilde{\textbf{p}} = \begin{bmatrix} \textbf{q}_1 \otimes \textbf{p}_1 \\ \textbf{q}_1 \otimes \textbf{p}_2 + \textbf{q}_2 \otimes \textbf{p}_1 \end{bmatrix}$


Бикватернионное умножение некоммутативно (при изменении порядка сомножителей результат бикватернионного умножения разный).

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

5. Обратный бикватернион

$\widetilde{\textbf{q}}^{-1} = \frac{\widetilde{\textbf{q}}^*}{\|\widetilde{\textbf{q}}\|}$



Определение бикватерниона через углы ориентации и вектор положения


Для начала определим системы координат, в которых будем рассматривать ориентацию и положение объекта в пространстве. Это необходимо сделать для задания действительной части бикватерниона (кватерниона ориентации), последовательность поворотов которого влияет на получаемый кватернион из углов ориентации. Здесь будем руководствоваться самолетными углами — рысканья $\psi$, тангажа $\vartheta$ и крена $\gamma$.

Определим базовую систему координат. Представьте, что вы стоите на поверхности Земли и смотрите в направлении Севера.

Точка Oo — начало системы координат, располагается в точке начала движения объекта.
Ось OoYg — направленна вертикально вверх, и противоположна к направлению вектора силы тяжести.
Ось OoXg — направлена в сторону Севера, по касательной местного меридиана.
Ось OoZg — дополняет систему до правой и направлена вправо, в сторону Востока.

Вторая система координат — связанная. Представьте, например, самолет или другой объект.
Точка O — начало системы координат, как правило, располагается в точке центра масс объекта.
Ось OY — направлена вертикально вверх, и перпендикулярна горизонтальной плоскости объекта.
Ось OX — направлена вперед, к передней точке объекта.
Ось OZ — дополняет систему до правой.

Положение объекта в пространстве задается радиус-вектором начала (точка O) связанной системы координат относительно неподвижной базовой системы координат. Ориентация связанной системы координат относительно базовой определяется тремя последовательными поворотами на:

угол рысканья $\psi$ — поворот вокруг оси OY,
угол тангажа $\vartheta$ — поворот вокруг оси OZ,
угол крена $\gamma$ — поворот вокруг оси OX.

Для первоначального определения бикватерниона необходимо задать действительную и дуальную части бикватерниона. Ориентация и положение объекта задается относительно некой базовой системы координат с помощью углов ориентации $\psi, \vartheta, \gamma$ и вектора положения центра масс $r = (r_x, r_y, r_z)^T$.

Действительную часть $\textbf{q}_1$ можно задать с помощью формулы:

$ \textbf{q}_1 = \begin{bmatrix} \cos\frac{\psi}{2} \cos\frac{\vartheta}{2} \cos\frac{\gamma}{2} & - & \sin\frac{\psi}{2} \sin\frac{\vartheta}{2} \sin\frac{\gamma}{2} \\ \cos\frac{\psi}{2} \cos\frac{\vartheta}{2} \sin\frac{\gamma}{2} & + & \sin\frac{\psi}{2} \sin\frac{\vartheta}{2} \cos\frac{\gamma}{2} \\ \cos\frac{\psi}{2} \sin\frac{\vartheta}{2} \sin\frac{\gamma}{2} & + & \sin\frac{\psi}{2} \cos\frac{\vartheta}{2} \cos\frac{\gamma}{2} \\ \cos\frac{\psi}{2} \sin\frac{\vartheta}{2} \cos\frac{\gamma}{2} & - & \sin\frac{\psi}{2} \cos\frac{\vartheta}{2} \sin\frac{\gamma}{2} \end{bmatrix} $


Обратите внимание, если у вас последовательность поворотов другая, то выражения будут тоже другими.

Дуальная часть $\textbf{q}_2$ определяется выражением:

$\textbf{q}_2 = \frac12 \textbf{r} \otimes \textbf{q}_1$



Вычисление углов ориентации и вектора положения из бикватерниона. Обратное преобразование


Вычислить углы ориентации можно из действительной части бикватерниона $\textbf{q}_1$:

$\psi = \arctan{\frac{2(q_0 q_2 - q_1 q_3)}{q_0^2 + q_1^2 - q_2^2 - q_3^2}}$


$\vartheta = \arcsin{(2(q_1 q_2 + q_0 q_3))}$


$\gamma = \arctan{\frac{2(q_0 q_1 - q_2 q_3)}{q_0^2 - q_1^2 + q_2^2 - q_3^2}}$



Положение объекта определяется выражением:

$\textbf{r} = 2 \textbf{q}_2 \otimes \textbf{q}_1^{-1}$


в результате получается вектор в кватернионной форме $\textbf{r} = (0, r_x, r_y, r_z)^T$

Поворот и перемещение вектора бикватернионом


Одно из замечательных свойств бикватернионов — это поворот и перемещение вектора из одной системы координат в другую. Пусть OoXgYgZg — неподвижная базовая система кординат, а OXYZ — связанная система координат объекта. Тогда ориентацию и положение объекта относительно базовой системы координат можно задать бикватернионом $\widetilde{\textbf{q}}$. Если задан вектор $\textbf{r}$ в связанной системе координат, тогда можно получить вектор $\textbf{r}_0$ в базовой системе координат с помощью формулы:

$\textbf{r}_0 = \widetilde{\textbf{q}} \otimes \textbf{r} \otimes \widetilde{\textbf{q}}^{-1}$


и обратно:

$\textbf{r} = \widetilde{\textbf{q}}^{-1} \otimes \textbf{r}_0 \otimes \widetilde{\textbf{q}}$


где $\textbf{r}$ — вектор в бикватернионной форме, $\textbf{r} = (1, 0, 0, 0, 0, r_x, r_y, r_z)$

JavaScript-библиотека работы с бикватернионами


Все вышеперечисленные операции работы с бикватернионами реализованы в javascript-библиотеке, в зависимости от ваших задач ее можно реализовать на других языках программирования. Основные функции работы с бикватернионами:
Функция Описание
DualQuaternion.dq Тело бикватерниона в виде массива 8-ми чисел
DualQuaternion(dq0, dq1, dq2, dq3, dq4, dq5, dq6, dq7) Конструктор, который определяет бикватернион путем заданием всех восьми чисел
DualQuaternion.fromEulerVector(psi, theta, gamma, v) Получить бикватернион путем задания ориентации объекта углами Эйлера и вектора положения объекта
DualQuaternion.getEulerVector() Получить углы Эйлера и вектор положения из бикватерниона
DualQuaternion.getVector() Получить вектор положения из бикватерниона
DualQuaternion.getReal() Получить действительную часть бикватерниона (определяет ориентацию объекта в пространстве)
DualQuaternion.getDual() Получить дуальную часть бикватерниона (определяет положение объекта в пространстве)
DualQuaternion.norm() Получить норму бикватерниона в виде дуального числа
DualQuaternion.mod() Получить модуль бикватерниона в виде дуального числа
DualQuaternion.conjugate() Получить сопряженный бикватернион
DualQuaternion.inverse() Получить обратный бикватернион
DualQuaternion.mul(DQ2) Бикватернионное умножение
DualQuaternion.toString() Преобразовать бикватернион в строчку, например, для вывода в отладочную консоль

Файл dual_quaternion.js
/**
 *
 * Author 2017, Akhramovich A. Sergey (akhramovichsa@gmail.com)
 * see https://github.com/infusion/Quaternion.js
 */

// 'use strict';

/**
 * Dual Quaternion constructor
 *
 * @constructor
 * @returns {DualQuaternion}
 */
function DualQuaternion(dq0, dq1, dq2, dq3,  dq4, dq5, dq6, dq7) {
	if (dq0 === undefined) {
		this.dq = [1, 0, 0, 0,  0, 0, 0, 0];
	} else {
		this.dq = [dq0, dq1, dq2, dq3,  dq4, dq5, dq6, dq7];
	}
	return this;
};

// Получение кватерниона по углам Эйлера и вектору положения
DualQuaternion['fromEulerVector'] = function(psi, theta, gamma, v) {
	var q_real = new Quaternion.fromEuler(psi, theta, gamma);
	var q_v    = new Quaternion(0, v[0], v[1], v[2]);
	var q_dual = q_v.mul(q_real);

	return new DualQuaternion(
		q_real.q[0],     q_real.q[1],     q_real.q[2],     q_real.q[3],
		q_dual.q[0]*0.5, q_dual.q[1]*0.5, q_dual.q[2]*0.5, q_dual.q[3]*0.5);
};

DualQuaternion.prototype = {
	'dq': [1, 0, 0, 0,  0, 0, 0, 0],
	
	/**
	 * Получение углов Эйлера (psi, theta, gamma) и вектора положения из бикватерниона
	 */
	'getEulerVector': function() {
		var euler_angles = this.getReal().getEuler();

		var q_dual = this.getDual();
		var q_dual_2 = new Quaternion(2.0*q_dual.q[0], 2.0*q_dual.q[1], 2.0*q_dual.q[2], 2.0*q_dual.q[3]);
		var q_vector = q_dual_2.mul(this.getReal().conjugate());

		return [euler_angles[0], euler_angles[1], euler_angles[2],
			q_vector.q[1], q_vector.q[2], q_vector.q[3]];
	},

	/**
	 * Получение только вектора положения из бикватерниона
	 */
	'getVector': function() {
		var euler_vector = this.getEulerVector();
		return [euler_vector[3], euler_vector[4], euler_vector[5]];
	},

	/**
	 * Получить действительную часть бикватерниона
	 * @returns {Quaternion}
	 */
	'getReal': function() {
		return new Quaternion(this.dq[0], this.dq[1], this.dq[2], this.dq[3]);
	},

	/**
	 * Получить дуальную часть бикватерниона
	 * @returns {Quaternion}
	 */
	'getDual': function() {
		return new Quaternion(this.dq[4], this.dq[5], this.dq[6], this.dq[7]);
	},

	/**
	 * Норма бикватерниона
	 * Внимание! Возвращает дуальное число!
	 */
	'norm': function() {
		return [Math.pow(this.dq[0], 2) + Math.pow(this.dq[1], 2) + Math.pow(this.dq[2], 2) + Math.pow(this.dq[3], 2),
			this.dq[0]*this.dq[4] + this.dq[1]*this.dq[5] + this.dq[2]*this.dq[6] + this.dq[3]*this.dq[7]];
	},

	/**
	 * Модуль бикватерниона
	 * Внимание! Возвращает дуальное число!
	 */
	'mod': function() {
		var q_real_mod = Math.sqrt(Math.pow(this.dq[0], 2) + Math.pow(this.dq[1], 2) + Math.pow(this.dq[2], 2) + Math.pow(this.dq[3], 2));
		return [q_real_mod,
			(this.dq[0]*this.dq[4] + this.dq[1]*this.dq[5] + this.dq[2]*this.dq[6] + this.dq[3]*this.dq[7])/q_real_mod];
	},

	/**
	 * Сопряженный бикватернион
	 * DQ' := (dq0, -dq1, -dq2, -dq3,  dq4, -dq5, -dq6, -dq7)
	 */
	'conjugate': function() {
		return new DualQuaternion(this.dq[0], -this.dq[1], -this.dq[2], -this.dq[3],  this.dq[4], -this.dq[5], -this.dq[6], -this.dq[7]);
	},

	// Вычислить обратный бикватернион
	'inverse': function() {
		var q_real_norm = new Quaternion(this.dq[0], this.dq[1], this.dq[2], this.dq[3]).norm();
		
		var dq_norm_inv = [q_real_norm, - (this.dq[0]*this.dq[4] + this.dq[1]*this.dq[5] + this.dq[2]*this.dq[6] + this.dq[3]*this.dq[7])/q_real_norm];

		var dq_conj = this.conjugate();
		
		// Умножение бикватерниона на дуальное число
		return new DualQuaternion(
			dq_norm_inv[0] * dq_conj.dq[0], 
			dq_norm_inv[0] * dq_conj.dq[1],
			dq_norm_inv[0] * dq_conj.dq[2],
			dq_norm_inv[0] * dq_conj.dq[3],
			dq_norm_inv[0] * dq_conj.dq[4] + dq_norm_inv[1] * dq_conj.dq[0],
			dq_norm_inv[0] * dq_conj.dq[5] + dq_norm_inv[1] * dq_conj.dq[1],
			dq_norm_inv[0] * dq_conj.dq[6] + dq_norm_inv[1] * dq_conj.dq[2],
			dq_norm_inv[0] * dq_conj.dq[7] + dq_norm_inv[1] * dq_conj.dq[3]);
	},

	/**
	 * Бикватернионное умножение
	 * q1_real*q2_real, q1_real*q2_dual + q1_dual*q2_real
	 */
	'mul': function(DQ2) {
		var q1_real = this.getReal();
		var q1_dual = this.getDual();
		var q2_real = DQ2.getReal();
		var q2_dual = DQ2.getDual();

		var q_res_real   = q1_real.mul(q2_real);
		var q_res_dual_1 = q1_real.mul(q2_dual);
		var q_res_dual_2 = q1_dual.mul(q2_real);

		return new DualQuaternion(
			q_res_real.q[0],
			q_res_real.q[1],
			q_res_real.q[2],
			q_res_real.q[3],
			q_res_dual_1.q[0] + q_res_dual_2.q[0],
			q_res_dual_1.q[1] + q_res_dual_2.q[1],
			q_res_dual_1.q[2] + q_res_dual_2.q[2],
			q_res_dual_1.q[3] + q_res_dual_2.q[3]);
	},

	/**
	 * Преобразование вектора бикватернионом
	 */
	'transformVector': function (v) {
		var dq_res = this.mul(new DualQuaternion(1, 0, 0, 0,  0, v[0], v[1], v[2])).mul(this.conjugate());

		return [dq_res.dq[5], dq_res.dq[6], dq_res.dq[7]];
	},

	/**
	 * Преобразовать в строку, для отладки
	 */
	'toString': function() {
		return '[' + 
			this.dq[0].toString() + ', ' + this.dq[1].toString() + ', ' + this.dq[2].toString() + ', ' + this.dq[3].toString() + ', ' +
			this.dq[4].toString() + ', ' + this.dq[5].toString() + ', ' + this.dq[6].toString() + ', ' + this.dq[7].toString() + ']';
	}

}

/*
// TEST:
var dq1 = new DualQuaternion.fromEulerVector(0 * Math.PI/180.0, 0 * Math.PI/180, 0 * Math.PI/180, [10, 20, 30]);

console.log(dq1.toString());
console.log('getEulerVector = ', dq1.getEulerVector());

console.log('norm = ', dq1.norm());
console.log('mod = ', dq1.mod());
console.log('conjugate = ', dq1.conjugate().dq);
console.log('inverse = ', dq1.inverse().dq);

var dq2 = new DualQuaternion.fromEulerVector(0 * Math.PI/180.0, 0 * Math.PI/180, 0 * Math.PI/180, [10, 0, 0]);
console.log('mul = ', dq1.mul(dq2).dq);

console.log('transformVector ??? = ', dq1.transformVector([0, 0, 0]));
*/


Пример работы с бикватернионами


Для лучшего понимания основ применения бикватернионов в качестве примера, рассмотрим небольшую игру. Задается прямоугольная область — карта. По карте плавает корабль, на котором установлено поворотное орудие. Здесь необходимо учесть, что для корабля базовая система координат является система координат карты, а для орудия базовая система координат — корабль. Все объекты отрисовываются в системе координат карты и здесь будет интересно увидеть, как можно перейти от системы координат орудия в систему координат карты с помощью свойства бикватернионного умножения. Управление движением корабля осуществляется клавишами W, A, S, D. Направление орудия задается курсором мышки.



Корабль и орудие описываются двумя классами: Ship и Gun. В конструкторе класса корабля задается его форма в виде бикватернионных точек, начальная ориентация и положение на карте в виде бикватерниона this.dq_pos.

Так же заданы бикватернионные приращения при управлении кораблем. При движении вперед-назад (клавиши W, S) будет изменяться только дуальная часть бикватерниона, а при управлении вправо-влево (клавиши A, D) будет меняться действительная и дуальная часть бикватерниона, которая задает угол поворота.

function Ship(ctx, v) {
    this.ctx    = ctx;
    this.dq_pos = new DualQuaternion.fromEulerVector(0*Math.PI/180, 0, 0, v);

    // Форма корабля
    this.dq_forward_left    = new DualQuaternion.fromEulerVector(0, 0, 0, [ 15, 0, -10]);
    this.dq_forward_right   = new DualQuaternion.fromEulerVector(0, 0, 0, [ 15, 0,  10]);
    this.dq_backward_left   = new DualQuaternion.fromEulerVector(0, 0, 0, [-15, 0, -10]);
    this.dq_backward_right  = new DualQuaternion.fromEulerVector(0, 0, 0, [-15, 0,  10]);
    this.dq_forward_forward = new DualQuaternion.fromEulerVector(0, 0, 0, [ 30, 0,  0]);

    // Приращения текущей позиции при управлении
    this.dq_dx_left     = new DualQuaternion.fromEulerVector( 1*Math.PI/180, 0, 0, [0, 0, 0]);
    this.dq_dx_right    = new DualQuaternion.fromEulerVector(-1*Math.PI/180, 0, 0, [0, 0, 0]);
    this.dq_dx_forward  = new DualQuaternion.fromEulerVector(0, 0, 0, [ 1, 0, 0]);
    this.dq_dx_backward = new DualQuaternion.fromEulerVector(0, 0, 0, [-1, 0, 0]);

    return this;
};

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

Ship.prototype = {
    'ctx': 0,
    'dq_pos':  new DualQuaternion.fromEulerVector(0, 0, 0, 0, 0, 0),

    /**
     * Нарисовать кораблик
     */
    'draw': function() {

        // Переместить все точки кораблика с помощью бикватернионного умножения
        v_pos             = this.dq_pos.getVector();
        v_forward_left    = this.dq_pos.mul(this.dq_forward_left).getVector();
        v_forward_right   = this.dq_pos.mul(this.dq_forward_right).getVector();
        v_backward_left   = this.dq_pos.mul(this.dq_backward_left).getVector();
        v_backward_right  = this.dq_pos.mul(this.dq_backward_right).getVector();
        v_forward_forward = this.dq_pos.mul(this.dq_forward_forward).getVector();

        // Непосредственно рисование
        ctx.beginPath();
        ctx.moveTo(v_backward_left[0],   v_backward_left[2]);
        ctx.lineTo(v_forward_left[0],    v_forward_left[2]);
        ctx.lineTo(v_forward_left[0],    v_forward_left[2]);
        ctx.lineTo(v_forward_forward[0], v_forward_forward[2]);
        ctx.lineTo(v_forward_right[0],   v_forward_right[2]);
        ctx.lineTo(v_backward_right[0],  v_backward_right[2]);
        ctx.lineTo(v_backward_left[0],   v_backward_left[2]);
        ctx.stroke();
        ctx.closePath();
    }
};

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

function Gun(ctx, ship, v) {
    this.ctx  = ctx;
    this.ship = ship;

    // Позиция орудия относительно корабля
    this.dq_pos = new DualQuaternion.fromEulerVector(0, 0, 0, v);

    // Форма орудия
    this.dq_forward  = new DualQuaternion.fromEulerVector(0, 0, 0, [20, 0, 0]);
    this.dq_backward = new DualQuaternion.fromEulerVector(0, 0, 0, [ 0, 0, 0]);

    // Вращение орудия при управлении
    this.dq_dx_left     = new DualQuaternion.fromEulerVector( 1*Math.PI/180, 0, 0, [0, 0, 0]);
    this.dq_dx_right    = new DualQuaternion.fromEulerVector(-1*Math.PI/180, 0, 0, [0, 0, 0]);

    return this;
};

В классе орудия также реализована только одна функция его отрисовки Ship.draw(). Орудие отображается в виде линии, которая задается двумя точками this.dq_backward и this.dq_forward. Для определения координат точек орудия применяется операция бикватернионного умножения.

Gun.prototype = {
    'ctx':  0,
    'ship': 0,
    'dq_pos': new DualQuaternion.fromEulerVector(0, 0, 0, [0, 0, 0]),

    /**
     * Нарисовать орудие
     */
    'draw': function() {

        // Переместить орудие относительно корабля
        v_pos        = this.ship.dq_pos.getVector();
        v_forward    = this.ship.dq_pos.mul(this.dq_backward).mul(this.dq_forward).getVector();
        v_backward   = this.ship.dq_pos.mul(this.dq_backward).getVector();
        
        // Непосредственно рисование
        ctx.beginPath();
        ctx.moveTo(v_backward[0], v_backward[2]);
        ctx.lineTo(v_forward[0],  v_forward[2]);
        ctx.stroke();
        ctx.closePath();
    }
};

Обработка управления кораблем и орудием реализована через события. За нажатие и отжатие клавиш управления кораблем отвечают четыре переменные leftPressed, upPressed, rightPressed, downPressed, которые обрабатываются в основном цикле программы.

leftPressed  = false;
rightPressed = false;
upPressed    = false;
downPressed  = false;
dq_mouse_pos = new DualQuaternion.fromEulerVector(0, 0, 0, [0, 0, 0]);

document.addEventListener("keydown",   keyDownHandler,   false);
document.addEventListener("keyup",     keyUpHandler,     false);
document.addEventListener("mousemove", mouseMoveHandler, false);

// Обработка нажатия клавиш управления
function keyDownHandler(e) {
    if      (e.keyCode == 37 || e.keyCode == 65 || e.keyCode == 97)  { leftPressed  = true; } // влево  A
    else if (e.keyCode == 38 || e.keyCode == 87 || e.keyCode == 119) { upPressed    = true; } // вверх  W
    else if (e.keyCode == 39 || e.keyCode == 68 || e.keyCode == 100) { rightPressed = true; } // вправо D
    else if (e.keyCode == 40 || e.keyCode == 83 || e.keyCode == 115) { downPressed  = true; } // вниз   S
}

// Обработка отжатия клавиш управления
function keyUpHandler(e) {
    if      (e.keyCode == 37 || e.keyCode == 65 || e.keyCode == 97)  { leftPressed  = false; } // влево  A
    else if (e.keyCode == 38 || e.keyCode == 87 || e.keyCode == 119) { upPressed    = false; } // вверх  W
    else if (e.keyCode == 39 || e.keyCode == 68 || e.keyCode == 100) { rightPressed = false; } // вправо D
    else if (e.keyCode == 40 || e.keyCode == 83 || e.keyCode == 115) { downPressed  = false; } // вниз   S
}

Одна из самых интересных функций, с точки зрения применения бикватернионных операций, это управление орудием корабля в направлении указателя мышки. Сначала координаты указателя мышки определяются в бикватернион dq_mouse_pos. Затем вычисляется бикватернион положения мышки относительно корабля с помощью бикватернионного умножения. От бикватерниона мышки отнимается бикватернион корабля dq_mouse_pos_about_ship = ship_1.dq_pos.inverse().mul(dq_mouse_pos);
(Прим.: операции последовательного бикватернионного умножения читайте справа налево). И наконец, определяется угол между векторами орудия и мышки. Начальной точке орудия gun_1.dq_backward присваивается полученное значение.

function mouseMoveHandler(e) {
    var relativeX = e.clientX - canvas.offsetLeft;
    var relativeY = e.clientY - canvas.offsetTop;

    // Обрабатывать события только когда курсор мышки находится в игровой области
    if (relativeX > 0 && relativeX < canvas.width && 
        relativeY > 0 && relativeY < canvas.height) {
        
        // Бикватернион положения мышки
        dq_mouse_pos = new DualQuaternion.fromEulerVector(0, 0, 0, [relativeX, 0, relativeY]);

        // Бикватернион положения мышки относительно корабля
        // Направление орудия. От координат мышки отнимается координаты корабля
        // Последовательность бикватернионного умножения важна
        // DQ_ship^(-1) * DQ_mouse
        dq_mouse_pos_about_ship = ship_1.dq_pos.inverse().mul(dq_mouse_pos);

        // Угол между векторами орудия и мышки
        q_gun_mouse = new Quaternion.fromBetweenVectors(gun_1.dq_forward.getVector(), dq_mouse_pos_about_ship.getVector());

        dq_gun_mouse = new DualQuaternion(q_gun_mouse.q[0], q_gun_mouse.q[1], q_gun_mouse.q[2], q_gun_mouse.q[3], 0, 0, 0, 0);

        gun_1.dq_backward = dq_gun_mouse;
        
        // console.log(dq_gun_mouse.getEulerVector());
        // console.log(relativeX + ' ' + relativeY + ' ' + gun_1.dq_forward.toString());
    }
}

В основном теле программы инициализируются объекты корабля и орудия ship_1 и gun_1, выводится отладочная информация и осуществляется обработка управления кораблем.

var canvas = document.getElementById("myCanvas");
var ctx    = canvas.getContext("2d");

ship_1 = new Ship(ctx, [100, 0, 100]);
gun_1  = new Gun(ctx, ship_1, [0, 0, 0]);

function draw() {
    ctx.clearRect(0, 0, canvas.width, canvas.height);
    ship_1.draw();
    gun_1.draw();

    // Debug info
    ship_euler_vector    = ship_1.dq_pos.getEulerVector();
    ship_euler_vector[0] = ship_euler_vector[0]*180/Math.PI;
    ship_euler_vector[1] = ship_euler_vector[1]*180/Math.PI;
    ship_euler_vector[2] = ship_euler_vector[2]*180/Math.PI;
    ship_euler_vector    = ship_euler_vector.map(function(each_element){ return each_element.toFixed(2); });
    ship_dq              = ship_1.dq_pos.dq.map(function(each_element){ return each_element.toFixed(2); });

    gun_dq = ship_1.dq_pos.mul(gun_1.dq_backward).dq.map(function(each_element){ return each_element.toFixed(2); });

    ctx.font = "8pt Courier";
    ctx.fillText("Ship: " + ship_dq + " | psi, theta, gamma, vector:" + ship_euler_vector, 10, 20);
    ctx.fillText("Gun:  " + gun_dq, 10, 40);

    // Управление корабликом
    if (leftPressed)  { ship_1.dq_pos = ship_1.dq_pos.mul(ship_1.dq_dx_left);     }
    if (rightPressed) { ship_1.dq_pos = ship_1.dq_pos.mul(ship_1.dq_dx_right);    }
    if (upPressed)    { ship_1.dq_pos = ship_1.dq_pos.mul(ship_1.dq_dx_forward);  }
    if (downPressed)  { ship_1.dq_pos = ship_1.dq_pos.mul(ship_1.dq_dx_backward); }

    requestAnimationFrame(draw);
}

draw();

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

Пример работы с бикватернионами

Заключение


У вас может возникнуть вопрос: зачем применять такой сложный математический аппарат, когда можно обойтись стандартными средствами для перемещения и вращения объектов? Одним из основных преимуществ заключается в том, что бикватернионная форма записи является более эффективной в вычислительном плане, так как все операции работы с бикватернионами после раскрытия выражений являются линейными. В данном видео Geometric Skinning with Approximate Dual Quaternion Blending показано насколько эффективнее вычисления с использованием бикватернионов в сравнении с другими методами.

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

  1. Челноков Юрий Николаевич. Кватернионные и бикватернионные модели и методы механики твердого тела и их приложения. Геометрия и кинематика движения. — монументальный теоретический труд.
  2. Гордеев Вадим Николаевич. Кватернионы и бикватернионы с приложениями в геометрии и механике. — написана более понятным языком и показаны применения в задачах формообразования криволинейных пространственных структур.

Комментарии (26)


  1. kovserg
    15.01.2019 21:56
    #19618564

    … более эффективной в вычислительном плане, так как все операции работы с бикватернионами после раскрытия выражений являются линейными.

    И где видно что количество операций сложения, умножения и обращений к памяти уменьшилось или хотя бы сравнение по производительности на тестовых примерах.
    Чем бикватернион лучше, просто вектора и кватерниона?
    Или в случае кораблей 2d вектора и одного угла?

    ps: Еще такой вопрос чем отличаются норма и модуль бикватерниона и почему они содержат символ Клиффорда?


    1. AkhramovichSA Автор
      16.01.2019 03:11
      #19619308
      +1

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

      В статье это не показано, так как это отдельная тема, но у меня есть публикация где такой подробный анализ проводится.
      Я привел ссылку на видео где показано сравнение моделирования 3D-графики между обычными методами и методами с использованием бикватернионов.
      Из своего опыта скажу, что при использовании бикватернионов вы убираете всю тригонометрию. Тригонометрические функции используются только для задания начальных условий, ну и для вывода отладочной информации в человекопонятном виде.
      Возможно стоит подробно раскрыть тему операции бикватернионного умножения в отдельной статье.


      Чем бикватернион лучше, просто вектора и кватерниона?
      Или в случае кораблей 2d вектора и одного угла?

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


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


      ps: Еще такой вопрос чем отличаются норма и модуль бикватерниона и почему они содержат символ Клиффорда?

      Бикватернион, по сути своей — дуальный кватернион q? = q1 + ?q2.
      Но часто его представляют в виде двух кватернионов, как это сделано в статье, поэтому здесь символ Клиффорда опущен, хотя его наличие надо всегда подразумевать у дуальной части q2.


      Если совсем утрировать, то норма — это сумма квадратов компонентов бикватерниона, а модуль — корень квадратный из суммы квадратов компонентов бикватерниона.
      Эти величины используются в других операциях, например, в операции вычисления обратного бикватерниона (в коде dual_quaternion.js в методе inverse как раз реализована операция деления бикватерниона на дуальное число), или при использовании численных методов интегрирования бикватерниона приходится его нормировать.


      1. Keyten
        16.01.2019 13:47
        #19621148

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


        1. AkhramovichSA Автор
          16.01.2019 14:05
          #19621254

          Я старался не делать громоздких выводов. Но, раскрытие этих формул есть в коде javascript-библиотеки, методы «norm» и «mod».


    1. Xop
      16.01.2019 10:54
      #19620122
      +1

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

      На самом деле — дуальные кватернионы немного "дороже" по операциям, чем просто пара "вектор переноса + обычный кватернион". Главное их преимущество — почти полное отсутствие артефактов при скиннинге в местах, где на вершины влияет более одной кости. В видео приведенном в статье про это рассказывается с 0:55. И там кстати видно, что интерполяция на дуальных кватернионов работает чуть медленнее (~200 FPS) чем "обычный" метод (~270 FPS)


  1. vlad9486
    15.01.2019 23:51
    #19618938
    +2

    Октонион — это другое. Эта штука в статье — дуальный кватернион.


    1. AkhramovichSA Автор
      16.01.2019 00:58
      #19619108
      +1

      Да, вы правы. Просто в статьях иногда эти понятия путают (из статьи убрал).
      Также неверно называть его «комплексный кватернион».


  1. netricks
    16.01.2019 11:19
    #19620286

    Большое спасибо за обзор. Он очень полезен.

    У меня такой вопрос:
    Есть ли преимущество в использовании бикватерниона перед преобразованием, описываемым кватернионом (для ориентации) и 3-вектором (для трансляции). Есть ли преимущество по скорости вычисления комбинированного перемещения?

    P.S. А, хотя вижу комментарий выше.


    1. AkhramovichSA Автор
      16.01.2019 13:52
      #19621174

      Я выложил статью на github в которой проведено такое исследование.
      Ben Kenwright — A Beginners Guide to Dual-Quaternions.pdf
      Примерный перевод главы «11. Результаты»:

      Бикватернион объединяет перемещение и вращение в одну переменную состояния. Это переменная состояния предлагает надежный, однозначный, вычислительно эффективный способ представления преобразования твердого тела.
      Вычислительные затраты различных матриц и бикватерниона:
      Матрица4х4: 64 умножений + 48 сложений.
      Матрица4х3: 48 умножений + 32 сложений.
      Бикватернион: 42 умножений + 38 сложений.

      В наших тестах мы обнаружили, что метод умножения бикватернионных преобразований в среднем на 10% быстрее по сравнению с умножением матриц. Мы не использовали преимущества архитектуры ЦП, используя параллельные методы, такие как SIMD, которые могут еще больше улучшить скорость вычислений, как продемонстрировал Паллави [MEHU10] (как для матриц, так и для умножения кватернионов).

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


      1. netricks
        16.01.2019 14:10
        #19621280

        Случай кватернион + 3вектор тут не рассматривается, но если подумать, комбинация преобразований в этой форме включает довольно затратный поворот 3вектора. Надо будет посчитать на досуге.


        1. lgorSL
          18.01.2019 01:59
          #19629488

          Ну пусть мы комбинируем две пары (поворот, вектор) (q1,r1) и (q2, r2).
          В результате получится что-то вроде (q2 * q1, r1 + q1 * r2 * q1^{-1}) (я мог местами перепутать индексы, но на рассуждения это не влияет.)


          1. перемножение кватернионов: вроде бы 16 умножений и 12 сложений.
          2. поворот вектора кватернионом можно реализовать как два умножения кватернионов. Там одна из координат будет 0, можно это использовать и не делать лишних умножнений. (https://habr.com/en/post/255005/) Получается, первое "умножение " можно сделать за 12 умножений и 8 сложений, инвертирование кватерниона 3 или даже 1 одно вычитание (сложение), смотря как инвертировать. При перемножении от результата нам вроде бы нужны только три компоненты (это же повернутый вектор будет), поэтому 12 умножений и 9 сложений. Сложение с вектором — ещё 3 сложения.

          Итого, если я не ошибаюсь, получается 40 умножений и 33 сложения. Хм. Капельку меньше, чем для бикватернионов.


          P.S. Вообще говоря, такие оценки не несут большой ценности, но мне самому было интересно. Важно насколько хорошо происходящее сочетается с инструкциями процессора, которые, например, позволяют вектор из четырёх чисел складывать с другим вектором или поточечно умножать — и тогда те же матрицы 4*4 могут оказаться не хуже чем 4*3 или бикватернионы.


      1. netricks
        16.01.2019 14:13
        #19621290

        А что за операция рассмотрена? Преобразование точки или комбинация?
        Надо будет почитать статью.

        P.S. Если преобразование точки, то 3вектор+кватернион однозначно медленнее матрицы 3х4.


      1. kovserg
        17.01.2019 21:13
        #19628702

        Смущают цифры

        Вычислительные затраты различных матриц и бикватерниона:
        Матрица4х4: 64 умножений + 48 сложений.
        Матрица4х3: 48 умножений + 32 сложений.
        Бикватернион: 42 умножений + 38 сложений.

        У меня получается примерно так:
        Умножение вектора на матрицу 3*4 = 9 умножений + 9 сложений
        Умножение матрицы на матрицу 3*4 = 36 умножений + 27 сложений
        Умножение кватерниона на кватернион = 10 умножений + 6 сложений
        Поворот вектора кватернионом = 15 умножений + 9 сложений
        Преобразование с поворотом и сдвигом через кватернион => 15 умножений + 15 сложений

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

        Может я где-то не прав?
        R'=A*R

x'=A00*x+A01*y+A02*z+A03
y'=A10*x+A11*y+A12*z+A13
z'=A20*x+A21*y+A22*z+A23

=> *9 +9

C=A*B

| C00 C01 C02 C03 |   | A00 A01 A02 A03 |   | B00 B01 B02 B03 |
| C10 C11 C12 C13 | = | A10 A11 A12 A13 | * | B10 B11 B12 B13 |
| C20 C21 C22 C23 |   | A20 A21 A22 A23 |   | B20 B21 B22 B23 |
| 0   0   0     1 |   | 0   0   0     1 |   | 0   0   0     1 |

C00 = A00*B00 + A01*B10 + A02*B20
C01 = A00*B01 + A01*B11 + A02*B21
C02 = A00*B02 + A01*B12 + A02*B22
C03 = A00*B03 + A01*B13 + A02*B23 + A03

C10 = A10*B00 + A11*B10 + A12*B20
C11 = A10*B01 + A11*B11 + A12*B21
C12 = A10*B02 + A11*B12 + A12*B22
C13 = A10*B03 + A11*B13 + A12*B23 + A03

C20 = A20*B00 + A21*B10 + A22*B20
C21 = A20*B01 + A21*B11 + A22*B21
C22 = A20*B02 + A21*B12 + A22*B22
C23 = A20*B03 + A21*B13 + A22*B23 + A03

=> *36 +27

R'=q*R/q

Rx'=(1-q1*q1)*Rx-q1*q2*Ry-q1*q3*Rz
Ry'=-q1*q2*Rx+(1-q2*q2)*Ry-q2*q3*Rz
Rz'=-q1*q3*Rx-q2*q3*Ry+(1-q3*q3)*Rz

=> *15 +9

R'=R1+q*(R-R0)/q
=> *15 +15


        1. 0serg
          17.01.2019 22:08
          #19628858

          Там речь судя по всему о комбинировании пребразований (умножение матриц / кватернионов / etc) а не о применении преобразования к вектору / векторам. Правда я как-то с трудом себе представляю где вопрос производительности комбинирования преобразований имел бы какое-то практическое значение.

          В Ваших расчетах вроде всё верно.


          1. MATov
            17.01.2019 23:01
            #19629030

            Правда я как-то с трудом себе представляю где вопрос производительности комбинирования преобразований имел бы какое-то практическое значение.

            В геймдеве, например.


            1. 0serg
              17.01.2019 23:05
              #19629042

              Сколько там матриц надо перемножить чтобы посчитать финальную? Десяток? Ну и на что ускорение подобного преобразования на 10% повлияет? Вот к миллионам векторов потом эту матрицу применять — там да, нужно время. Но как раз там-то как верно заметил kovserg быстрее матриц ничего нет. И по числу операций и по удобству приложения к современным вычислительным платформам.


  1. slopestyler
    16.01.2019 14:40
    #19621448

    Видел применение дуальных чисел для расчёта производных. Интересно, можно ли тут это использовать и быстро считать производные по компонентам кватернионов.


    1. netricks
      16.01.2019 15:10
      #19621612

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


      1. lgorSL
        18.01.2019 02:11
        #19629510

        нет, там ничего избыточного не производится.
        (f, f') * (g, g') = (f * g, f * g' + f' * g)
        Умножением для дуальных чисел будет обычное умножение и сразу же вычисление производной от произведения. Если в будущем планируется использовать и то и другое, то этот способ оптимальнее некуда.


        1. netricks
          18.01.2019 10:20
          #19630286

          Если посчитать количество операций на расчет функции + производной и на расчет функции в дуальных числах, окажется, что второе вычислительно тяжелее


          1. lgorSL
            18.01.2019 11:44
            #19630762

            Можно пример? Мне кажется, что будет то же самое.


            1. netricks
              18.01.2019 12:24
              #19631060

              Возьмём функцию (x^2)*(x+3). Вычислим ее производную в точке 'a' в дуальных числах.
              x->(a,1)
              x^2 -> (a,1)*(a,1) -> (a*a, a*1 + a*1) -> (b, bi)
              x+3 -> (a,1) + (3,0) -> (a+3, 1+0) -> (c, ci)
              x^2 * (x+3) -> (b,bi)*(c,ci) -> (b*c, b*ci + c*bi) -> (d, di)
              result: (d, di).imag -> di
              Итого: 6 умножений, 4 сложения.

              Вычисляем предберя производную: 2x*(x+3) + x^2 = 3*x*x + 6*x = 3*x*(x+2)
              Итого: 1 сложение, 2 умножения. (и столько же, если мы хотим и значение функции тоже (Итого: 2 сложения, 4 умножения))

              Проблема в том, что, если вы погоните дуальные числа через кватернионы, то каждая операция умножения (коих в кватернионном умножении 16) будет превращаться в три. Стоимость кватернионного умножения в дуальных числах — 48 аппаратных умножений, не считая сложений.


    1. mayorovp
      16.01.2019 15:41
      #19621774

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


      1. lgorSL
        18.01.2019 02:06
        #19629504

        А если сделать "дуальные-дуальные-кватернионы"? Ну просто смотреть на дуальные кватернионы как на один объект, задающий ориентацию в пространстве, и поверх прикрутить дуальность для производных. Мне кажется, что такое может работать.


        1. netricks
          18.01.2019 10:58
          #19630468

          Это будет работать, но, как я пишу выше, это огроменный оверхед по производительности. Это вдвойне излишне, учитывая, что обычно в контексте задач кинематики (для которых обычно применяются кватернионы и бикватернионы)
          g'(M*x) == M*g'(x) — производная трансформированного вектора равна трансформированной производной вектора (хотя, это, конечно… не всегда)…


  1. dMac
    17.01.2019 17:34
    #19627700

    Капец я тупой…