Матрицы, векторы, или истинный путь познания OpenGL

Изменено: . Разместил: FadeDemon | Меню
Количество просмотров3`820 / FD1856-6-0

Итак, вновь выдалась пара свободных дней, что-бы что-нибудь здесь написать. А всё потому, что вчера наконец-то сдался курсач по объектно-ориентированному программированию. И так получилось, что этому курсачу выпала честь стать темой сегодняшней статьи. Казалось-бы, при чем тут OpenGL? Однако, обо-всём по порядку…


Применение OpenGL к лабораторным по ООП
Рис. №973. Применение OpenGL к лабораторным по ООП

Выбирая тему для проекта, случайно я вспомнил одну винрарную игрушку — Colony Wars: Red Sun, и захотелось мне внезапно… Да, вы абсолютно правы, написать свой космический симулятор. Который бы по графике и возможностям, как минимум не уступал той, что была на PSX, а местами, даже предоставлял несколько бо́льшие возможности. Забегая вперед, надо сказать, что полностью задумку реализовать не удалось по причине нехватки времени патологической лени, однако в результате трехнедельных каторжных трудов получилась довольно внушительная платформа, на которой, в принципе, в будущем, можно на самом деле соорудить какую-нибудь незамысловатую космическую стрелялку. Платформа получилась на столько внушительной, что о некоторых вещах мне и захотелось поделиться с общественностью. Итак, обо всём вкратце и по порядку.


Шесть степеней свободы и матрицы


Даже после первых раздумий быстро обозначилась первая проблема: реализация честных 6-ти степеней свободы, без всякой халтуры. Всяких-там углов Эйлера и подобных способов хранения вращения объектов здесь уже было явно недостаточно, ну, или по крайней мере, такой способ был бы очень неудобным, особенно когда речь заходила бы о “конвертировании” этих самых углов в векторы направления. Поэтому, взгляд мой быстро направился в сторону матриц, поскольку именно квадратные матрицы размерностью 4 являются основной боевой единицей во всей идеологии OpenGL.


\[ \bold M_v = \left (\begin{matrix} x_0 & x_1 & x_2 & 0 \\ y_0 & y_1 & y_2 & 0 \\ z_0 & z_1 & z_2 & 0 \\ x_p & y_p & z_p & 1 \\ \end{matrix}\right ) \]
Матрица видовой модели

Давайте разберёмся, что-же такого интересного в этих матрицах есть, и зачем они вообще нужны. Если у вас есть опыт работы с OpenGL, то вам, наверно известны такие замечательные функции, как glRotate(), glTranslate(), glScale() и т.п. Их довольно удобно использовать, однако, задумывались ли вы, а что эти функции на самом деле делают? Функции эти, на самом деле, всего-лишь изменяют видовую матрицу. В дальнейшем, рассуждения, вообще говоря, пойдут только о видовой матрице, то-есть предполагается, что была вызвана glMatrixMode(GL_MODELVIEW). А изменяют они её при помощи такой интересной штуки, как аффинных преобразований, или если выражаться проще, путем умножения видовой матрицы на матрицу, соответствующую нужному преобразованию. Давайте разберемся, глядя на картинку выше, что же за магия в этой матрице храниться. На самом деле, всё достаточно просто: в видовой матрице храниться текущая система координат, представленная векторами {x0, y0, z0} (ось X), {x1, y1, z1} (ось Y), {x2, y2, z2} (ось Z), и смещение этой системы координат в пространстве относительно нулевой точки, представленное вектором {xp, yp, zp}. Нулевая точка, кстати, в данном случае, это точка, в которой находиться виртуальная “камера” сцены, на виду которой мы все собственно рисуем. Тут легко заметить, что единичная матрица, как этого и можно было ожидать, задаёт ортогональную систему координат, с единичными осями, то-есть такую, в которой мы все привыкли работать. Кстати надо напомнить, что по умолчанию в режиме перспективы (после вызова gluPerspective()) оси в OpenGL направлены так:

• Ось X — вправо от центра экрана — →
• Ось Y — вверх от центра экрана — ↑
• Ось Z — на нас, из глубины монитора

Этот простой факт несложно проверить, если, например, сделать видовую матрицу единичной, и просто попытаться рисовать, например, треугольники. В ортогональном режиме, направление осей зависит от порядка переданных в функцию glOrtho() аргументов, и, например, при таком порядке — glOrtho(0, 640, 0, 480, -100, 100) — таки совпадает с перспективным.


Как можно заметить, всё действительно просто и понятно. А как насчёт аффинных преобразований? На самом деле, там всё тоже довольно просто. Оказывается, имея систему координат в таком 4-х мерном виде, мы можем очень легко творить с ней всё что нам вздумается с помощью обыкновенного матричного умножения, просто умножая нашу матрицу на матрицы специального вида, задающие, собственно, наши преобразования. Преобразований, кстати, можно сразу назвать 3 штуки: смещение системы координат (по векторам её собственных осей), вращение её вокруг произвольных осей, и масштабирование. Чуть позже, мы рассмотрим их все. А пока, стоит заметить — для реализации шести степеней свободы, аффинных преобразований нам будет более чем достаточно.


Умножение матриц, и SSE4


Вообще, во всей этой системе с матрицами, есть одна не очень хорошая особенность — это сложность операции матричного умножения. Как мы знаем, умножение матриц подразумевает нахождение всех возможных скалярных произведений строк и столбцов этих матриц, которых в 4-х мерном случае, как не трудно заметить, будет целых 16 штук. Учитывая что на одно скалярное произведение приходиться 4 операции умножения, и 3 сложения, мы получаем угрожающие 64 умножения, и 48 сложений. На две только матрицы. На до ли говорить, как быстро это будет работать при большом количестве таких операций… Впрочем интересное решение есть.


Если делать умножение матриц “в тупую” просто написав, так сказать, длинные формулы, то компилятор это скорее всего скомпилирует в FPU-команды, которые не отличаются расторопностью, и которых действительно будет 64 + 48 штук. Однако, мы можем вспомнить, что в арсенале наших процессоров есть такая штука, как SSE, а в частности, такая команда, как DPPS, выполняющая ни что иное, как скалярное произведение векторов. Причем, если рассматривать её вариант с 4-х байтовыми float'ами, то в одно скалярное произведение как раз и влезает 4 операнда. Прямо то, что нам нужно. Учитывая, что эта команда делает 4 пакетных умножения, и 3 сложения, то чтобы перемножить матрицы, достаточно будет 16 раз вызвать эту волшебную команду, и просто собрать из кусочков результат. Я не замерял, но, думаю, вы сами догадаетесь, во сколько раз это будет быстрее наивного варианта. Реализацию такого умножения, мы здесь и приведём.


Для начала, мы вспомним, что SSE очень требовательно относится к выравниванию передаваемых ему адресов, поэтому сделаем вот такую штуку:


Код: C++
class T4Matrix {
	private:
		float d[20];
		float *m;
	...
}

void __fastcall T4Matrix::_init_mem(){
	this->m = (float *) (((unsigned long) & this->d[0] + 0x10) & ~0xf);
}

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


Код: C++
T4Matrix & __fastcall T4Matrix::operator *= (const T4Matrix & r){
	T4Matrix t(r);
	T4Matrix s(*this);
	s._transpose();

	float *p = this->m;
	float *sm = s.m;
	float *tm = t.m;

	// SSE4 matrix multiplication
	__asm {
		mov eax, sm;
		mov ecx, tm;
		mov edx, p;


		// load row 1
		movaps xmm0, [eax + 0x00]
		movaps xmm1, xmm0
		movaps xmm2, xmm0
		movaps xmm3, xmm0
		// multiply
		dpps xmm0, [ecx + 0x00], 0xf1
		dpps xmm1, [ecx + 0x10], 0xf1
		dpps xmm2, [ecx + 0x20], 0xf1
		dpps xmm3, [ecx + 0x30], 0xf1
		// integry result
		shufps xmm0, xmm0, 57
		movss xmm0, xmm1
		shufps xmm0, xmm0, 57
		movss xmm0, xmm2
		shufps xmm0, xmm0, 57
		movss xmm0, xmm3
		shufps xmm0, xmm0, 57
		// store row 1
		movaps [edx + 0x00], xmm0


		// load row 2
		movaps xmm0, [eax + 0x10]
		movaps xmm1, xmm0
		movaps xmm2, xmm0
		movaps xmm3, xmm0
		// multiply
		dpps xmm0, [ecx + 0x00], 0xf1
		dpps xmm1, [ecx + 0x10], 0xf1
		dpps xmm2, [ecx + 0x20], 0xf1
		dpps xmm3, [ecx + 0x30], 0xf1
		// integry result
		shufps xmm0, xmm0, 57
		movss xmm0, xmm1
		shufps xmm0, xmm0, 57
		movss xmm0, xmm2
		shufps xmm0, xmm0, 57
		movss xmm0, xmm3
		shufps xmm0, xmm0, 57
		// store row 2
		movaps [edx + 0x10], xmm0


		// load row 3
		movaps xmm0, [eax + 0x20]
		movaps xmm1, xmm0
		movaps xmm2, xmm0
		movaps xmm3, xmm0
		// multiply
		dpps xmm0, [ecx + 0x00], 0xf1
		dpps xmm1, [ecx + 0x10], 0xf1
		dpps xmm2, [ecx + 0x20], 0xf1
		dpps xmm3, [ecx + 0x30], 0xf1
		// integry result
		shufps xmm0, xmm0, 57
		movss xmm0, xmm1
		shufps xmm0, xmm0, 57
		movss xmm0, xmm2
		shufps xmm0, xmm0, 57
		movss xmm0, xmm3
		shufps xmm0, xmm0, 57
		// store row 3
		movaps [edx + 0x20], xmm0


		// load row 4
		movaps xmm0, [eax + 0x30]
		movaps xmm1, xmm0
		movaps xmm2, xmm0
		movaps xmm3, xmm0
		// multiply
		dpps xmm0, [ecx + 0x00], 0xf1
		dpps xmm1, [ecx + 0x10], 0xf1
		dpps xmm2, [ecx + 0x20], 0xf1
		dpps xmm3, [ecx + 0x30], 0xf1
		// integry result
		shufps xmm0, xmm0, 57
		movss xmm0, xmm1
		shufps xmm0, xmm0, 57
		movss xmm0, xmm2
		shufps xmm0, xmm0, 57
		movss xmm0, xmm3
		shufps xmm0, xmm0, 57
		// store row 3
		movaps [edx + 0x30], xmm0
	}

	this->transpose();
	return *this;
}

Здесь, как мы видим, присутствует транспонирование, и прежде чем понять, зачем оно здесь нужно, подумаем, как его можно тоже сделать быстрым. На самом деле, мне не пришло в голову ничего быстрее, чем вот такая реализация:


Код: C++
void __fastcall T4Matrix::_transpose(){
	float *p = this->m;
	// ultra-fast transposing
	__asm {
		mov eax, p

		mov edx, [eax + 1 * 4 + 0 * (4 * 4)]
		xchg edx, [eax + 0 * 4 + 1 * (4 * 4)]
		mov [eax + 1 * 4 + 0 * (4 * 4)], edx

		mov edx, [eax + 2 * 4 + 0 * (4 * 4)]
		xchg edx, [eax + 0 * 4 + 2 * (4 * 4)]
		mov [eax + 2 * 4 + 0 * (4 * 4)], edx

		mov edx, [eax + 3 * 4 + 0 * (4 * 4)]
		xchg edx, [eax + 0 * 4 + 3 * (4 * 4)]
		mov [eax + 3 * 4 + 0 * (4 * 4)], edx

		mov edx, [eax + 2 * 4 + 1 * (4 * 4)]
		xchg edx, [eax + 1 * 4 + 2 * (4 * 4)]
		mov [eax + 2 * 4 + 1 * (4 * 4)], edx

		mov edx, [eax + 3 * 4 + 1 * (4 * 4)]
		xchg edx, [eax + 1 * 4 + 3 * (4 * 4)]
		mov [eax + 3 * 4 + 1 * (4 * 4)], edx

		mov edx, [eax + 3 * 4 + 2 * (4 * 4)]
		xchg edx, [eax + 2 * 4 + 3 * (4 * 4)]
		mov [eax + 3 * 4 + 2 * (4 * 4)], edx
	}
}

Здесь, просто, нужно сразу понять, что менять местами надо на самом деле не так и много элементов (6 штук всего), и сделать это, конечно-же, можно без всяких циклов. А теперь, зачем оно нам нужно. Всё дело в том, что произведение матриц не коммутативно, то-есть, от перестановки их местами в произведении, меняется результат. А первичный вариант этого кода умножал не совсем в том порядке, в котором делает это, например, glMultMatrix(), но поскольку код мне переделывать было лень, вспомнилась формула:


\[ (\bold {A \cdot B})^\tau = \bold {B}^\tau \cdot \bold {A}^\tau \]
Свойства произведения

Как видим, достаточно транспонировать, и всё будет окей. Правда всё усугублялось тем, что OpenGL, по всей видимости, видовую матрицу хранит транспонированной, но в итоге код выше заработал отлично, и позволил полностью отказаться от glTranslate(), glRotate() и прочих функций. И да, такой подход позволял хранить вращения объектов в виде матриц, отказавшись от углов Эйлера, и подобных вещей. Давайте теперь разберёмся с реализацией аффинных преобразований.


Коротко об аффинных преобразованиях


На самом деле, здесь всё просто. Рассмотрим матрицы, которые задают операции сдвига (translate) и масштабирования (scale):


\[ \bold M_t = \left (\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ x & y & z & 1 \\ \end{matrix}\right )\, \, \, \, \, \, \, \, \bold M_s = \left (\begin{matrix} x & 0 & 0 & 0 \\ 0 & y & 0 & 0 \\ 0 & 0 & z & 0 \\ 0 & 0 & 0 & 1 \\ \end{matrix}\right ) \]
Матрицы сдвига и масштабирования

От нас требуется, всего-лишь, вместо x, y, z подставить нужные нам значения, например длины смещений по соответствующим осям, или нужные масштабы, а матричное умножение сделает за нас всё само. Но на самом деле, это не самое главное: эти 2 преобразования вполне легко реализуются и без матриц. Давайте мы рассмотрим ещё один вид преобразования, это поворот нашей системы координат вокруг произвольной оси в пространстве (отдельное спасибо Википедии за эту формулу):


\[ M(\vec{\mathbf{v}},\theta) = \begin{pmatrix} \cos \theta + (1 - \cos \theta) x^2 & (1 - \cos \theta) x y - (\sin \theta) z & (1 - \cos \theta) x z + (\sin \theta) y \\ (1 - \cos \theta) y x + (\sin \theta) z & \cos \theta + (1 - \cos \theta) y^2 & (1 - \cos \theta) y z - (\sin \theta) x \\ (1 - \cos \theta) z x - (\sin \theta) y & (1 - \cos \theta) z y + (\sin \theta) x & \cos \theta + (1 - \cos \theta) z^2 \end{pmatrix} \]
Матрица поворота вокруг произвольной оси

Вид здесь, конечно, немного посложнее, но от нас по прежнему требуется просто подставить x, y, z и θ (угол поворота, в радианах в данном случае), и перемножить матрицы. И система координат почти что магическим образом, покорно поворачивается. Даже можно привести простенькую, не претендующую на идеал, реализацию:


Код: C++
void __fastcall T4Matrix::rotate(const T4Vector & xiss, float ang, bool is_multiply){
	// sse = off
	float co = cosf(ang * M_PI / 180);
	float co1 = 1.0f - co;
	float si = sinf(ang * M_PI / 180);

	T4Matrix ro;
	float *du = ro.data();
	float *x = xiss.data();

	// compose a rotation matrix (by wikipedia)
	du[0 + 0 * 4] = co + co1 * x[0] * x[0];
	du[0 + 1 * 4] = co1 * x[0] * x[1] - si * x[2];
	du[0 + 2 * 4] = co1 * x[0] * x[2] + si * x[1];

	du[1 + 0 * 4] = co1 * x[1] * x[0] + si * x[2];
	du[1 + 1 * 4] = co + co1 * x[1] * x[1];
	du[1 + 2 * 4] = co1 * x[1] * x[2] - si * x[0];

	du[2 + 0 * 4] = co1 * x[2] * x[0] - si * x[1];
	du[2 + 1 * 4] = co1 * x[2] * x[1] + si * x[0];
	du[2 + 2 * 4] = co + co1 * x[2] * x[2];

	// multiply matrix
	if(is_multiply){
		*this *= ro;
	}else{
		*this = ro;
	}
}

Вот здесь, наступает момент, который явился ключевым в программировании шести степеней свободы. А именно: ведь не сложно догадаться, что можно вытаскивать из текущей матрицы векторы-направляющие координатных осей, и подсовывать их
той матрице поворота. Догадываетесь, что мы получим? А получим мы очень простую и лёгкую реализацию трех степеней свободы вращения: крен, тангаж и рысканье. Соответствуют эти вещи, поворотам вокруг оси, которая параллельна нашему взгляду (крен), вокруг оси, которая идет от нас вверх (рысканье), и которая предыдущим двум ортогональна (тангаж). Чтобы было понятнее, рысканье (yaw) соответствует поворотам головы влево-вправо, крен (roll) соответствует наклонам головы влево-вправо, а тангаж (pitch), соответствует опусканию взгляда вниз (например, когда мы смотрим под ноги), и поднятию его вверх. Что касается реализации движения в пространстве, то здесь всё по аналогии. Чтобы сделать движение вперед, по направлению взгляда, вытаскиваем из матрицы соответствующий вектор, умножаем его на скорость, и например, на fDeltaTime, и плюсуем к текущим координатам, которые, кстати, лучше всё-таки хранить отдельно от матрицы. Для космического корабля, который мы хотели запрограммировать, будет вполне достаточно, ибо не думаю, что последний должен уметь делать стрейфы, и смещаться вертикально вверх-вниз. Впрочем, реализовать это будет столь-же просто.


Типичный кубик Рубика на OpenGL
Рис. №977. Один из примеров удачного применения матриц повотора

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


Первая проблема, которая очень быстро настигнет нашу матрицу, это потери точности на арифметических операциях (вернее их накопление со временем), которые усугубляются тем, что тип чисел у нас пока-что float. Можно, конечно, переделать и на double, и на long double, но SSE к таким числам уже не прикрутишь. А тем временем, есть отличный способ, как от накопления ошибок избавиться.


Нам потребуется: Векторное произведение, 1 штука, Нормализация вектора, 1 штука. В чём суть? Для начала, разберемся в чём, собственно, проблема. Ошибки в точности при вращении нашей системы вокруг осей приводят к тому, что оси координат со временем теряют свою ортогональность, и начинают разъезжаться в длине, то-есть, их длина перестает быть равна единице, как у нормальных векторов. Однако, мы можем этому легко противостоять! Для начала вспомним, что если мы уже имеем 2 некоторых вектора, то ортогональный к ним, мы легко можем получить с помощью векторного произведения. А получив ортогональный двум другим вектор, мы можем найти ещё один ортогональный вектор, который, заменив собой один из исходных векторов, и будучи ортогональным двум другим ортогональным, приведёт систему координат в нужное нам ортогональное состояние. Остаётся только перед этим не забыть все векторы нормализовать. Или после этого, порядок этих действий совершенно не важен :).


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


\[ \overline{a} \times \overline{b} = \begin{vmatrix} \overline{\mathbf i} & \overline{\mathbf j} & \overline{\mathbf k} \\ a_x & a_y & a_z \\ b_x & b_y & b_z \end{vmatrix} = \left [\overline{\mathbf i}\left (a_y b_z - a_z b_y\right ) + \overline{\mathbf j}\left (a_z b_x - a_x b_z\right ) + \overline{\mathbf k}\left (a_x b_y - a_y b_x\right )\right ] \]
Векторное произведение, и как его можно расписать

Не сразу очевидно, но посчитать его можно очень быстро с помощью всё того-же SSE, использую поразительную мощность команды SHUFPS. Главное — правильно расставить маски перемешивания. Пример:


Код: C++
// vector product
T4Vector & __fastcall T4Vector::operator %= (const T4Vector & a){
	float *d = this->m;
	float *s = a.m;

	__asm {
		mov edx, d
		mov eax, s

		movaps xmm0, [eax]
		movaps xmm1, [edx]
		movaps xmm2, xmm0
		movaps xmm3, xmm1

		shufps xmm0, xmm0, 11001001b
		shufps xmm2, xmm2, 11010010b
		shufps xmm1, xmm1, 11010010b
		shufps xmm3, xmm3, 11001001b

		mulps xmm0, xmm1
		mulps xmm2, xmm3

		subps xmm0, xmm2

		movaps [edx], xmm0
	}

	return *this;
}

Что касается нормализации, то это уж совсем просто. Напомним, чтобы нормализовать вектор, то-есть, привести его к единичной длине, нужно найти его длину, и все компоненты на неё разделить. Кстати, тут очень кстати приходятся возможности интересной такой команды RSQRTSS, которая делает ни что иное, как быстро высчитывает x−½, что позволяет снова экономить наносекунды. Делается всё это безобразие, например, так:


Код: C++
void __fastcall T4Vector::normalize(){
	float *c = this->m;
	__asm {
		mov eax, c
		movaps xmm0, [eax]
		movaps xmm1, xmm0
		dpps xmm0, xmm0, 0xf1
		// check if zero
		xorps xmm2, xmm2
		ucomiss xmm0, xmm2
		jle getfuck
			rsqrtss xmm0, xmm0
			shufps xmm0, xmm0, 0
			mulps xmm1, xmm0
			movaps [eax], xmm1
		getfuck:
	}
}

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


Код: C++
void __fastcall T4Matrix::orthonormalize(){
	this->_transpose();
	T4Vector einz(this->m);
	// % is a vector product!
	T4Vector drei((this->m + 4) % einz);
	T4Vector zwei(einz % drei);

	einz.normalize();
	zwei.normalize();
	drei.normalize();

	this->set(0, einz);
	this->set(1, zwei);
	this->set(2, drei);
	this->_transpose();
}

Единственное, что не очень удобно, чтобы из матрицы быстро вытащить какую-либо конкретную ось, её надо предварительно транспонировать. Но это мелочи. Мы, зато, теперь умеем искусно фехтовать матрицами, и заставлять их делать то, что нам хочется :). А это главное…


Арбузы, и адский хаос
Рис. №974. Арбузы, и адский хаос

Ray casting. Расстреливаем треугольники из лучемёта


Вторая тема, которую хотелось-бы осветить, мало общего имеет с OpenGL, но является, на самом деле, не менее, важной, поскольку используется чуть менее, чем везде: от концепций рендеринга (Raycasting, Raytracing), до физики в игрушках вроде Counter-Strike. Суть задачи, которую мы сейчас рассмотрим, очень проста. Пусть в пространстве дан некоторый треугольник ABC в пространстве, и отрезок DE, там-же. Требуется, найти точку пересечения этого отрезка, и треугольника, или установить, что они не пересекаются. Как можно догадаться, вариации этой задачи как раз-таки и встречаются чуть менее, чем везде, начиная с того-же Counter-Strike. Ведь мы-же никогда не задумывались, а как на самом деле рассчитываются точки попаданий пуль в стены, или в головы врагов, когда мы даём очередь из M249? Сейчас, мы попробуем эту задачу решить.


Итак, начать нужно с того, что задачу можно разделить на 2 подзадачи: нахождение точки пересечения прямой и плоскости, которые задают отрезок и треугольник, и определения факта, принадлежит полученная точка треугольнику и отрезку, или нет. Почему так? А потому-что алгоритм решения первой задачи нам уже известен, и заключается он в нахождении решения системы 3-х линейных уравнений, уравнения плоскости, которую мы пересекаем, и 2-х уравнений плоскостей, которые в пересечении задают пересекающую прямую. Сложно? На самом деле всё просто. Итак, приступим.


\[ Ax + By + Cz + D = 0 \]
Нормальное уравнение плоскости

Итак, мы, кажется, знаем, что плоскость задаётся вот таким-вот симпатичным на вид уравнением, а котором A, B, C и D — некоторые коэффициенты, однозначно определяющие плоскость. Это конечно, замечательно, но сам собой возникает вопрос, а откуда, чёрт возьми, эти долбанные коэффициенты брать? Ведь мир так устроен, что у нас есть только 3 точки, которые задают наш треугольник, и должны как-то задавать и плоскость. Здесь на помощь приходит чёрная магия вуду 4 лаконичных формулы, которые мне удалось нагуглить, правда линк на оригинал навсегда потерян в bad-секторах моей памяти. В общем, отбросив лирику, мы просто скажем, что коэффициенты уравнения плоскости по трём точкам ищутся с помощью раскрытия вот таких симпатичных определителей:


\[ A = \left |\begin{matrix} 1 & y_1 & z_1 \\ 1 & y_2 & z_3 \\ 1 & y_3 & z_3 \\ \end{matrix}\right | \;\; B = \left |\begin{matrix} x_1 & 1 & z_1 \\ x_2 & 1 & z_3 \\ x_3 & 1 & z_3 \\ \end{matrix}\right | \;\; C = \left |\begin{matrix} x_1 & y_1 & 1 \\ x_2 & y_2 & 1 \\ x_3 & y_3 & 1 \\ \end{matrix}\right | \;\; D = -\left |\begin{matrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_3 & y_3 & z_3 \\ \end{matrix}\right | \;\; \]
Коэффициенты уравнения плоскости по трём точкам

Этот ужас, на самом деле, всё таки сто́ит попытаться раскрыть, чтобы сократить количество арифметических операций, и в итоге мы получим приблизительно это:

\[ A = y_1(z_2 - z_3) + y_2(z_3 - z_1) + y_3(z_1 - z_2) \]

\[ B = z_1(x_2 - x_3) + z_2(x_3 - x_1) + z_3(x_1 - x_2) \]

\[ C = x_1(y_2 - y_3) + x_2(y_3 - y_1) + x_3(y_1 - y_2) \]


Коэффициент D мы раскрывать не будем, а просто заметим, что он является определителем матрицы из трех векторов координат наших точек, взятый со знаком минус. Лучше, напомним, как, собственно, раскрывается определитель:

\[ \left |\begin{matrix} a_1 & b_1 & c_1 \\ a_2 & b_2 & c_2 \\ a_3 & b_3 & c_3 \end{matrix}\right | = a_1(b_2 c_3 - b_3 c_2) + a_2 (b_3 c_1 - b_1 c_3) + a_3 (b_1 c_2 - b_2 c_1) \]
Определитель матрицы 3-го порядка

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


Код: C++
float __fastcall T4Matrix::det3() const {
	float zu;
	float *ma = this->data();
	__asm {
		mov eax, ma

		movaps xmm1, [eax + 0x10]
		movaps xmm2, [eax + 0x20]

		movaps xmm0, [eax + 0x00]

		movaps xmm3, xmm1		// B
		movaps xmm4, xmm2		// C

		shufps xmm1, xmm1, 11001001b
		shufps xmm3, xmm3, 11010010b

		shufps xmm2, xmm2, 11010010b
		shufps xmm4, xmm4, 11001001b

		mulps xmm1, xmm2
		mulps xmm3, xmm4

		subps xmm1, xmm3
		dpps xmm0, xmm1, 0x71

		movss zu, xmm0
	}
	return zu;
}

Конечно-же, это очень напоминает векторное произведение, поскольку оно, как раз-таки, построено на определителе. Те-же SHUFPS, с почти такими-же масками, и никакой чёрной магии. Впрочем, для раскрытия первых трех определителей коэффициентов, нам лучше будет написать ещё более шуструю функцию, которая будет работать по раскрытым формулам. Будет она примерно такой:


Код: C++
T4Vector __fastcall T4Matrix::plane_eq_coeffs3() const {
	T4Vector o;
	T4Matrix tp = *this;
	tp._transpose();

	float *tpp = tp.data();
	float *opp = o.data();

	// find A, B, C coefficients in plane equation
	// three is reduced determinants fast calculation
	__asm {
		mov eax, tpp
		mov edx, opp

		movaps xmm1, [eax + 0x00]		// x
		movaps xmm2, [eax + 0x10]		// y
		movaps xmm0, [eax + 0x20]		// z

		movaps xmm3, xmm0				// z'
		movaps xmm4, xmm1				// z'
		movaps xmm5, xmm2				// y'

		shufps xmm0, xmm0, 11001001b
		shufps xmm1, xmm1, 11001001b
		shufps xmm2, xmm2, 11001001b

		shufps xmm3, xmm3, 11010010b
		shufps xmm4, xmm4, 11010010b
		shufps xmm5, xmm5, 11010010b

		subps xmm0, xmm3
		subps xmm1, xmm4
		subps xmm2, xmm5

		dpps xmm0, [eax + 0x10], 0x71	// xmm0[0] = A
		dpps xmm1, [eax + 0x20], 0x71	// xmm1[0] = B
		dpps xmm2, [eax + 0x00], 0x71	// xmm2[0] = C

		// integrity
		shufps xmm0, xmm0, 57
		movss xmm0, xmm1
		shufps xmm0, xmm0, 57
		movss xmm0, xmm2
		// shift by 2
		shufps xmm0, xmm0, 01001110b

		movaps [edx], xmm0
	}

	// D coeff
	o[3] = -this->det3();
	return o;
}

Ураа! Мы знаем коэффициенты в уравнении плоскости. Но постойте. Плоскость — это, конечно, хорошо. Но ведь у нас есть ещё и прямая, правда? Ээмм… Опять надо искать какие-то там коэффициенты. Но здесь, вроде-бы, всё более понятно. Метод поиска коэффициентов двух уравнений плоскостей, задающих прямую, по её направляющему вектору, или двум точкам, был мной выведен самостоятельно.


\[ {{x - x_0} \over {x_1 - x_0}} = {{y - y_0} \over {y_1 - y_0}} = {{z - z_0} \over {z_1 - z_0}} \]
Каноническое уравнение прямой в пространстве

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


\[ {{x - x_0} \over {x_1 - x_0}} = {{y - y_0} \over {y_1 - y_0}},\;\;{{y - y_0} \over {y_1 - y_0}} = {{z - z_0} \over {z_1 - z_0}} \]


Далее, для удобства записи, обозначим: N = x1 − x0, M = y1 − y0, K = z1 − z0, и попробуем с полученными уравнениями что-нибудь сделать:


\[ {{x - x_0} \over {N}} = {{y - y_0} \over {M}},\;\;\;\;{{y - y_0} \over {M}} = {{z - z_0} \over {K}} \] \[ {M(x - x_0) } = {N(y - y_0)},\;\;\;\;{K(y - y_0)} = {M(z - z_0)} \] \[ Mx - Mx_0 = Ny - Ny_0,\;\;\;\;Ky - Ky_0 = Mz - Mz_0 \] \[ \underline{Mx -Ny = Mx_0 - Ny_0},\;\;\;\;\underline{Ky - Mz = Ky_0 - Mz_0} \]

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


\[ \left\lbrace \begin{aligned} & Ax & + & \,By & + & \,Cz & = & \,-D \\ & (y_1 - y_0)x & + & \,(x_1 - x_0)y & + & \,0z & = & \,(y_1 - y_0)x_0 - (x_1 - x_0)y_0 \\ & 0x & + & \,(z_1 - z_0)y & + & \,(y_1 - y_0)z & = & \,(z_1 - z_0)y_0 - (y_1 - y_0)z_0 \end{aligned} \right. \]
Система уравнений, задающая пересечение плоскости, и прямой

Здесь надо пояснить, что коэффициенты A, B, C, D подставляются те самые, что мы нашил выше из определителей, а векторы {x0, y0, z0} и {x1, y1, z1} есть ни что иное, как координаты точек, которыми задаётся наша прямая. Всё, что остаётся нам сделать, это просто взять, и решить нашу систему уравнений. Ну, или установить, что она не решается. Я, вообще, не имею ничего против метода Гаусса, но в данном случае, нам вполне хватит старых добрых формул Крамера, учитывая что определители реализованы десятком SSE инструкций.


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


\[ \left\lbrace \begin{aligned} A_0x + B_0y + C_0z = D_0 \\ A_1x + B_1y + C_1z = D_1 \\ A_2x + B_2y + C_2z = D_2 \end{aligned} \right. \]
Типичная система линейных уравнений

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


\[ \Delta = \left |\begin{matrix} A_0 & B_0 & C_0 \\ A_1 & B_1 & C_1 \\ A_2 & B_2 & C_2 \end{matrix}\right |\;\; \Delta_x = \left |\begin{matrix} D_0 & B_0 & C_0 \\ D_1 & B_1 & C_1 \\ D_2 & B_2 & C_2 \end{matrix}\right |\;\; \Delta_y = \left |\begin{matrix} A_0 & D_0 & C_0 \\ A_1 & D_1 & C_1 \\ A_2 & D_2 & C_2 \end{matrix}\right |\;\; \Delta_z = \left |\begin{matrix} A_0 & B_0 & D_0 \\ A_1 & B_1 & D_1 \\ A_2 & B_2 & D_2 \end{matrix}\right |\;\; \]
Определители, необходимые в методе Крамера

А как только определители будут посчитаны, неизвестные x, y и z мы сразу-же можем отыскать по очень лаконичным и понятным формулам:


\[ x = {\Delta_x \over \Delta}\;\;\;y = {\Delta_y \over \Delta}\;\;\;z = {\Delta_z \over \Delta} \]
Решение методом Крамера

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


Код: C++
bool __fastcall T4Matrix::cramer3(const T4Vector & freem, T4Vector & solution) const {
	float det = this->det3();
	if(fabs(det) <= +0.0f){
		// no solution
		return false;
	}

	T4Matrix t1(*this);
	t1.transpose();
	T4Matrix t2(t1);

	t2.set(0, freem);
	solution[0] = t2.det3() / det;

	t2 = t1;
	t2.set(1, freem);
	solution[1] = t2.det3() / det;

	t2 = t1;
	t2.set(2, freem);
	solution[2] = t2.det3() / det;

	return true;
}

Переменная freem здесь есть ни что иное, как передаваемый вектор свободных членов (D0, D1, D2), поскольку матрица у нас, хоть и 4 × 4, функцией она рассматривается как матрица 3 × 3, а свободные члены откуда-то брать надо.


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


Код: C++
// intersection with plane
bool __fastcall T4Vector::intersect(const T4Vector & a, const T4Vector & b, 
	const T4Vector & c, const T4Vector & p, T4Vector & w) const {

	// compose matrix
	T4Matrix op(a, b, c);

	// find A, B, C coefficients in plane equation
	T4Vector coe(op.plane_eq_coeffs3());
	T4Vector p2(*this + p);

	// free members vector
	T4Vector freem(
		-coe[3],
		(p2[1] - p[1]) * p[0] - (p2[0] - p[0]) * p[1],
		(p2[2] - p[2]) * p[1] - (p2[1] - p[1]) * p[2]
	);

	// compose a linear equations system matrix, and find solution via Cramer's formulas
	op.set(0, coe);
	op.set(1, T4Vector(p2[1] - p[1], p[0] - p2[0], 0));
	op.set(2, T4Vector(0, p2[2] - p[2], p[1] - p2[1]));

	// TRY TO FIND ANY SOLUTION!!11
	return op.cramer3(freem, w);
}

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


Код: C++
float __fastcall T4Vector::length() const {
	float *c = this->m;
	float r;
	__asm {
		mov eax, c
		movaps xmm0, [eax]
		dpps xmm0, xmm0, 0xf1
		sqrtss xmm0, xmm0
		movss r, xmm0
	}
	return r;
}

bool T4Vector::intheline(const T4Vector & p1, const T4Vector & p2) const {
	float lila = (p1 - p2).length();
	return (p1 - *this).length() < lila && (p2 - *this).length() < lila;
}

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


\[ \theta\left (\vec a, \vec b\right ) = \arccos \left ({{\vec a \cdot \vec b} \over {\left |a\right | \cdot \left |b\right |}}\right ) \]
Угол между векторами, в радианах

С точки зрения кода, я реализовал это, перегрузив очень подходящий по смыслу оператор ^, и получилось примерно вот так:


Код: C++
float __fastcall operator ^ (const T4Vector & l, const T4Vector & r){
	float *d = l.m;
	float *s = r.m;

	float rt;

	__asm {
		mov edx, d
		mov eax, s

		movaps xmm0, [edx]
		movaps xmm1, [eax]
		movaps xmm2, xmm0

		// dot product
		dpps xmm2, xmm1, 0xf1

		// lengths
		dpps xmm0, xmm0, 0xf1
		dpps xmm1, xmm1, 0xf1

		// reverse square root
		rsqrtss xmm0, xmm0
		rsqrtss xmm1, xmm1

		mulss xmm2, xmm0
		mulss xmm2, xmm1

		movss rt, xmm2
	}
	return rt;
}

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


Определяем, лежит-ли точка в треугольнике…
Рис. №972. Определяем, лежит-ли точка в треугольнике…

И что-же мы видим здесь? Во первых, нам не важно, что точки у нас расположены в пространстве, условия поставленной задачи гарантируют, что проверяемая точка находится в той-же плоскости, что и три точки, задающие треугольник. Поэтому, нам ничего не мешает перейти к двумерному варианту задачи. А здесь мы имеем следующие 2 факта:

1. Точка лежит внутри угла, образуемого двумя векторами только тогда, когда оба угла между вектором, идущим от центра угла к точке, меньше, чем угол между двумя исходными векторами.
Предполагается, что угол всегда меньше 180°, иначе понятие “внутри угла” стало бы слишком расплывчатым.

2. Если для двух углов в треугольнике выполняется условие (1), то точка гарантированно лежит внутри этого треугольника.

Не знаю даже, как назвать эти два утверждения, но они, в целом, отлично работают, и пример на картинке даже демонстрирует это. Действительно, видно, что точка P1 не принадлежит треугольнику, а точка P2 принадлежит. И при этом, углы BCP2 и P2CA меньше, чем BCA, а углы BAP2 и P2AC меньше, чем BAC. А если мы посмотрим на точку P1, то сразу-же и заметим, что угол P1AC очевидно больше угла BAC, что нашей гипотезе, в целом, и соответствует. Что-же касается второго утверждения, то тут надо просто увидеть, что пересечение областей видимости любых двух углов, исходный треугольник собой и образует, и проверять условие для третьего угла совершенно не нужно. Думаю, это легко заметить, даже просто посмотря на картинку. Осталось последнее: код, который всё это делает. Вот, собственно, он:


Код: C++
// затычка, возникшая из-за проблем с точностью
float __fastcall acosex(float arg){
	if(arg < -1.0f) arg = -1.0f;
	if(arg > 1.0f) arg = 1.0f;
	return acos(arg);
}


bool T4Vector::intriangle(const T4Vector & a, const T4Vector & b, const T4Vector & c) const {

	T4Vector ab = b - a;
	T4Vector ac = c - a;

	T4Vector ca = a - c;
	T4Vector cb = b - c;

	T4Vector ap = *this - a;
	T4Vector cp = *this - c;

	// base angles
	float aco1 = acosex(ab ^ ac);
	float aco2 = acosex(ca ^ cb);

	return acosex(ac ^ ap) < aco1 && acosex(ap ^ ab) < aco1 && 
		acosex(ca ^ cp) < aco2 && acosex(cp ^ cb) < aco2;

}

Немного полезного кода


Напоследок, выложим сюда полную версию кода модуля, который отвечал у меня за математику, и реализован, как я уже писал, с использованием SSE4. Представлена вся математика там двумя классами: T4Vector и T4Matrix, в которые инкапсулируется много всяких полезных функций и операторов. Думаю, что всё это может оказаться полезным для тех, кто ищет действительно быструю реализацию матриц и векторов для OpenGL.

Модуль матрично-векторной математики на SSE4

На этом, дорогие товарищи, и закончим. До скорых встреч :) Напоследок, картинки, демонстрирующие тот самый мой почти готовый курсач:


Самодельный космический симулятор
Рис. №975. Самодельный космический симулятор

Космический симулятор, и взрыв
Рис. №976. Космический симулятор, и взрыв
Ключевые слова Программизм, математика, OpenGL, космический симулятор, матрица, DPPS, крен, тангаж, рысканье, Векторное произведение, Нормализация, RSQRTSS, метод Крамера, SSE4, аффинные преобразования
QR-код
Интересное
А теперь, внимание, вопрос: а можно ли обменять значение двух переменных так, чтобы третью явно или косвенно не использовать? Оказывается, можно.

Читать »»
Случайные фото