Грам-Шмидт неправильная реализация ортогонализации

Я нахожусь в процессе создания бесплатного движка 3D-игр с открытым исходным кодом на основе OpenGL3 (это не школьное задание, а скорее для развития личных навыков и для того, чтобы что-то вернуть сообществу открытого исходного кода). Я достиг стадии, когда мне нужно выучить много смежной математики, поэтому я читаю отличный учебник под названием «Математика для программирования 3D-игр и компьютерной графики, 3-е издание».

Я рано наткнулся на загадку, пытаясь выполнить упражнения из книги, поскольку моя попытка реализовать «алгоритм ортогонализации Грамма-Шмидта» в C ++ выдает неправильный ответ. Я не эксперт по математике (хотя я пытаюсь стать лучше), и у меня очень ограниченный опыт изучения математического алгоритма и его перевода в код (ограниченный некоторыми вещами, которые я узнал от Udacity.com). В любом случае, было бы действительно полезно, если бы кто-то мог взглянуть на мой неправильный код и дать мне подсказку или решение.

Вот:

/*
The Gram-Schmidt Orthogonalization algorithm is as follows:

Given a set of n linearly independent vectors Beta = {e_1, e_2, ..., e_n},
the algorithm produces a set Beta' = {e_1', e_2', ..., e_n'} such that
dot(e_i', e_j') = 0 whenever i != j.

A. Set e_1' = e_1
B. Begin with the index i = 2 and k = 1
C. Subtract the projection of e, onto the vectors e_1', e_2', ..., e_(i-1)'
from e_i, and store the result in e_i', That is,

dot(e_i, e_k')
e_i' = e_i - sum_over(-------------- e_k')
e_k'^2

D. If i < n, increment i and loop back to step C.
*/

#include <iostream>
#include <glm/glm.hpp>

glm::vec3 sum_over_e(glm::vec3* e, glm::vec3* e_prime, int& i)
{
int k = 0;
glm::vec3 result;

while (k < i-2)
{
glm::vec3 e_prime_k_squared(pow(e_prime[k].x, 2), pow(e_prime[k].y, 2), pow(e_prime[k].z, 2));
result += (glm::dot(e[i], e_prime[k]) / e_prime_k_squared) * e_prime[k];
k++;
}

return result;
}

int main(int argc, char** argv)
{
int n = 2;  // number of vectors we're working with
glm::vec3 e[] = {
glm::vec3(sqrt(2)/2, sqrt(2)/2, 0),
glm::vec3(-1, 1, -1),
glm::vec3(0, -2, -2)
};

glm::vec3 e_prime[n];
e_prime[0] = e[0];  // step A

int i = 0;  // step B

do  // step C
{
e_prime[i] = e[i] - sum_over_e(e, e_prime, i);

i++;    // step D
} while (i-1 < n);

for (int loop_count = 0; loop_count <= n; loop_count++)
{
std::cout << "Vector e_prime_" << loop_count+1 << ": < "<< e_prime[loop_count].x << ", "<< e_prime[loop_count].y << ", "<< e_prime[loop_count].z << " >" << std::endl;
}

return 0;
}

Этот код выводит:

Vector e_prime_1: < 0.707107, 0.707107, 0 >
Vector e_prime_2: < -1, 1, -1 >
Vector e_prime_3: < 0, -2, -2 >

но правильный ответ должен быть:

Vector e_prime_1: < 0.707107, 0.707107, 0 >
Vector e_prime_2: < -1, 1, -1 >
Vector e_prime_3: < 1, -1, -2 >

Редактировать: Вот код, который дает правильный ответ:

#include <iostream>
#include <glm/glm.hpp>

glm::vec3 sum_over_e(glm::vec3* e, glm::vec3* e_prime, int& i)
{
int k = 0;
glm::vec3 result;

while (k < i-1)
{
float e_prime_k_squared = glm::dot(e_prime[k], e_prime[k]);
result += ((glm::dot(e[i], e_prime[k]) / e_prime_k_squared) * e_prime[k]);
k++;
}

return result;
}

int main(int argc, char** argv)
{
int n = 3;  // number of vectors we're working with
glm::vec3 e[] = {
glm::vec3(sqrt(2)/2, sqrt(2)/2, 0),
glm::vec3(-1, 1, -1),
glm::vec3(0, -2, -2)
};

glm::vec3 e_prime[n];
e_prime[0] = e[0];  // step A

int i = 0;  // step B

do  // step C
{
e_prime[i] = e[i] - sum_over_e(e, e_prime, i);

i++;    // step D
} while (i < n);

for (int loop_count = 0; loop_count < n; loop_count++)
{
std::cout << "Vector e_prime_" << loop_count+1 << ": < "<< e_prime[loop_count].x << ", "<< e_prime[loop_count].y << ", "<< e_prime[loop_count].z << " >" << std::endl;
}

return 0;
}

1

Решение

Проблема, вероятно, в том, как вы определяете e_k'^2, Что касается векторной математики, то квадрат вектора обычно считается квадратом его нормы. Следовательно,

double e_prime_k_squared = glm::dot(e_prime_k, e_prime_k);

Более того, деление на вектор не определено (интересно, почему GLM это позволяет?), Так что если e_k'^2 это вектор, все это не определено.

2

Другие решения

Других решений пока нет …

По вопросам рекламы [email protected]