nerdcruft/libcruft-util
1
0
Fork 0
Code Issues Releases Activity

vector: add hughes-moeller make_basis

Browse Source
This commit is contained in:
Danny Robson 2018-04-18 21:46:26 +10:00
parent 999ec1c35c
commit b8253536ea
1 changed files with 40 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

49
vector.hpp
View File

@ -17,11 +17,13 @@
#ifndef CRUFT_UTIL_VECTOR_HPP
#define CRUFT_UTIL_VECTOR_HPP
#include "./coord/fwd.hpp"
#include "./coord.hpp"
#include "coord/fwd.hpp"
#include "coord/ops.hpp"
#include "coord.hpp"
#include "maths.hpp"
#include "debug.hpp"
#include "json/fwd.hpp"
#include "maths.hpp"
#include <cstddef>
#include <cmath>
@ -80,8 +82,6 @@ namespace util {
//-------------------------------------------------------------------------
// given a vector find two vectors which produce an orthonormal basis.
//
// we use frisvad's method, avoids explicit normalisation. a good
// alternative is hughes-moeller, but the paper is hard to find.
template <typename T>
std::pair<
util::vector<3,T>,
@ -89,6 +89,10 @@ namespace util {
>
make_basis (const util::vector<3,T> n)
{
#if 0
// frisvad's method avoids explicit normalisation. a good alternative
// is hughes-moeller, but the paper is hard to find.
CHECK (is_normalised (n));
// avoid a singularity
if (n.z < -T(0.9999999)) {
@ -101,10 +105,37 @@ namespace util {
const T a = 1 / (1 + n.z);
const T b = -n.x * n.y * a;
return {
{ 1 - n.x * n.x * a, b, -n.x },
{ b, 1 - n.y * n.y * a, -n.y }
};
const util::vector<3,T> v0 { 1 - n.x * n.x * a, b, -n.x };
const util::vector<3,T> v1 { b, 1 - n.y * n.y * a, -n.y };
CHECK (is_normalised (v0));
CHECK (is_normalised (v1));
return { v0, v1 };
#else
// huges-moeller isn't as fast, but is more accurate
if(util::abs (n.x) > util::abs (n.z))
{
// Normalization factor for b2
auto const a = rsqrt (n.x * n.x + n.y * n.y);
util::vector<3,T> b1 { -n.y * a, n.x * a, 0 };
// Cross product using that b2 has a zero component
util::vector<3,T> b0 { b1.y * n.z, -b1.x * n.z, b1.x * n.y - b1.y * n.x };
return { b0, b1 };
}
else
{
// Normalization factor for b2
auto const a = rsqrt (n.y * n.y + n.z * n.z);
util::vector<3,T> b1 { 0.0f, -n.z * a, n.y * a };
// Cross product using that b2 has a zero component
util::vector<3,T> b0 { b1.y * n.z - b1.z * n.y, b1.z * n.x, -b1.y * n.x };
return { b0, b1 };
}
#endif
}