source: GTP/trunk/App/Demos/Vis/CHC_revisited/Vector3.cpp @ 2751

#include "Matrix4x4.h"
#include "Vector3.h"
#include <float.h>



namespace CHCDemo
{

// Given min a vector to minimize and a candidate vector, replace
// elements of min whose corresponding elements in Candidate are
// smaller.  This function is used for finding objects' bounds,
// among other things.
void Minimize(Vector3 &min, const Vector3 &Candidate)
{
        if (Candidate.x < min.x)
                min.x = Candidate.x;
        if (Candidate.y < min.y)
                min.y = Candidate.y;
        if (Candidate.z < min.z)
                min.z = Candidate.z;
}

// Given min a vector to minimize and a candidate vector, replace
// elements of min whose corresponding elements in Candidate are
// larger.  This function is used for finding objects' bounds,
// among other things.
void Maximize(Vector3 &max, const Vector3 &Candidate)
{
        if (Candidate.x > max.x)
                max.x = Candidate.x;
        if (Candidate.y > max.y)
                max.y = Candidate.y;
        if (Candidate.z > max.z)
                max.z = Candidate.z;
}

// Project the vector onto the YZ, XZ, or XY plane depending on which.
// which      Coordinate plane to project onto
// 0          YZ
// 1          XZ
// 2          XY
// This function is used by the polygon intersection code.

void Vector3::ExtractVerts(float *px, float *py, int which) const
{
        switch (which) {
        case 0:
                *px = y;
                *py = z;
                break;
        case 1:
                *px = x;
                *py = z;
                break;
        case 2:
                *px = x;
                *py = y;
                break;
        }
}

// returns the axis, where the vector has the largest value
int
Vector3::DrivingAxis(void) const
{
        int axis = 0;
        float val = fabs(x);

        if (fabs(y) > val) {
                val = fabs(y);
                axis = 1;
        }

        if (fabs(z) > val)
                axis = 2;

        return axis;
}

// returns the axis, where the vector has the smallest value
int
Vector3::TinyAxis(void) const
{
        int axis = 0;
        float val = fabs(x);

        if (fabs(y) < val) {
                val = fabs(y);
                axis = 1;
        }

        if (fabs(z) < val)
                axis = 2;

        return axis;
}


// Construct a view vector ViewN, and the vector ViewU perpendicular
// to ViewN and lying in the plane given by ViewNxUpl
// the last vector of ortogonal system is ViewV, that is
// perpendicular to both ViewN and ViewU. 
// |ViewN| = |ViewU| = |ViewV| = 1
// The ViewN vector pierces the center of the synthesized image
//     ViewU vector goes from the center image rightwards
//     ViewV vector goes from the center image upwards
void
ViewVectors(const Vector3 &DirAt, const Vector3 &Viewer,
                        const Vector3 &UpL, Vector3 &ViewV, Vector3 &ViewU,
                        Vector3 &ViewN)
{
        Vector3 U, V, N;
        Vector3 Up = Normalize(UpL);

        N = -Normalize(DirAt);

        V = Normalize(Up - DirAt);
        V -= N * DotProd(V, N);
        V = Normalize(V);
        U = CrossProd(V, N);

        ViewU = U; // rightwards
        ViewV = V; // upwards
        ViewN = -N; // forwards
#ifdef GTP_DEBUG
        const float eps = 1e-3f;
        if (fabs(Magnitude(ViewU) - 1.0) > eps) {
                Debug << "ViewU magnitude error= " << Magnitude(ViewU) << "\n";
        }
        if (fabs(Magnitude(ViewV) - 1.0) > eps) {
                Debug << "ViewU magnitude error= " << Magnitude(ViewV) << "\n";
        }
        if (fabs(Magnitude(ViewN) - 1.0) > eps) {
                Debug << "ViewU magnitude error= " << Magnitude(ViewN) << "\n";
        }
#endif // GTP_DEBUG

        return;
}

// Given the intersection point `P', you have available normal `N'
// of unit length. Let us suppose the incoming ray has direction `D'.
// Then we can construct such two vectors `U' and `V' that
// `U',`N', and `D' are coplanar, and `V' is perpendicular
// to the vectors `N','D', and `V'. Then 'N', 'U', and 'V' create
// the orthonormal base in space R3.
void
TangentVectors(Vector3 &U,
                           Vector3 &V, // output
                           const Vector3 &normal, // input
                           const Vector3 &dirIncoming)
{
#ifdef GTP_DEBUG
        float d = Magnitude(normal); 
        if ( (d < 0.99) ||
                (d > 1.01) ) {
                        Debug << " The normal has not unit length = " << d << endl;
        }
        d = Magnitude(dirIncoming); 
        if ( (d < 0.99) ||
                (d > 1.01) ) {
                        Debug << " The incoming dir has not unit length = " << d << endl;
        }
#endif

        V = CrossProd(normal, dirIncoming);

        if (SqrMagnitude(V) < 1e-3) {
                // the normal and dirIncoming are colinear
                // we can/have to generate arbitrary perpendicular vector to normal.
                if (fabs(normal.x) < 0.6)
                        V.SetValue(0.0, -normal.z, normal.y);
                else {
                        if (fabs(normal.y) < 0.6)
                                V.SetValue(-normal.z, 0.0, normal.x);
                        else
                                V.SetValue(-normal.y, normal.x, 0.0);
                }
        }
        V = Normalize(V);

        U = CrossProd(normal, V);
#ifdef GTP_DEBUG
        d = SqrMagnitude(U);
        if ( (d < 0.99) ||
                (d > 1.01) ) {
                        Debug << "The size of U vector incorrect\n";
        }
#endif
        return;  
}


Vector3 UniformRandomVector()
{
        float r1 = RandomValue(0.0f, 1.0f);
        float r2 = RandomValue(0.0f, 1.0f);

        return UniformRandomVector(r1, r2);
}


Vector3 UniformRandomVector(const float r1, const float r2)
{
        float cosTheta = 1.0f - 2*r1;
        float sinTheta = sqrt(1 - sqrt(cosTheta));
        float fi = 2.0f * M_PI * r2;

        Vector3 dir(sinTheta*sin(fi), cosTheta, sinTheta * cos(fi));

        return dir;
}


Vector3 CosineRandomVector(const Vector3 &normal)
{
        float r1, r2;

        r1 = RandomValue(0.0f, 1.0f);
        r2 = RandomValue(0.0f, 1.0f);

        return CosineRandomVector(r1, r2, normal);
}


Vector3 CosineRandomVector(const float r1,
                                                   const float r2,
                                                   const Vector3 &normal)
{
        float theta = 2.0f * M_PI * r2;
        float radius = sqrt(r1);

        float x = radius * sin(theta);
        float z = radius * cos(theta);
        float y = sqrt(1.0f - x*x - z*z);

        Vector3 dir(x, y, z);

        //  return Normalize(dir);

        Matrix4x4 m = RotationVectorsMatrix(normal, Vector3(0,1,0));

        Matrix4x4 mi = Invert(m);

        m = m * RotationVectorsMatrix(Vector3(0,1,0), Normalize(dir)) * mi;

        return TransformNormal(m, normal);
}


Vector3 UniformRandomVector(const Vector3 &normal)
{
        float r1 = RandomValue(0.0f, 1.0f);
        float r2 = RandomValue(0.0f, 1.0f);
        
        float cosTheta = 1.0f - r1;
        float sinTheta = sqrt(1 - sqr(cosTheta));
        float fi = 2.0f*M_PI*r2;

        Vector3 dir(sinTheta*sin(fi),
                cosTheta,
                sinTheta*cos(fi));

        Matrix4x4 m = RotationVectorsMatrix(normal, Vector3(0,1,0));
        Matrix4x4 mi = Invert(m);

        m = m*RotationVectorsMatrix(Vector3(0, 1, 0), Normalize(dir)) * mi;

        return TransformNormal(m, normal);
}


bool Vector3::CheckValidity() const
{
#ifdef _WIN32
        return !(_isnan(x) || _isnan(y) || _isnan(z));
#else
        return !(isnan(x) || isnan(y) || isnan(z));
#endif  
}

}