---

VecUtil.cc

Go to the documentation of this file.

/*
    File:           VecUtil.cc
    
    Function:       Contains a grab-bag of useful graphics
                    vector routines.

    Author(s):      Andrew Willmott

*/

#include "gcl/VecUtil.h"

Bool Refract(
        GCLReal         fromIndex, 
        GCLReal         toIndex,
        const Vector    &v,            
        const Vector    &n,            
        GCLReal         cosNV,         
        Vector          &refractDir
    )
{
    GCLReal kn, cos2, k;
    
    kn = fromIndex / toIndex;
    cos2 = 1.0 - sqr(kn) * (1.0 - sqr(cosNV));

    if (cos2 < 0.0)
        return(1);          // Total internal reflection!
            
    k = kn * cosNV - sqrt(cos2);
    
    refractDir = -kn * v + k * n;

    return(0);
}

Void CalcTriAreaNormal(
        const Point &a,
        const Point &b, 
        const Point &c,
        Vector      &n
    )
{
    n[0] = (a[1] - b[1]) * (a[2] + b[2]) + (b[1] - c[1]) * (b[2] + c[2])
        + (c[1] - a[1]) * (c[2] + a[2]);
    n[1] = (a[2] - b[2]) * (a[0] + b[0]) + (b[2] - c[2]) * (b[0] + c[0])
        + (c[2] - a[2]) * (c[0] + a[0]);
    n[2] = (a[0] - b[0]) * (a[1] + b[1]) + (b[0] - c[0]) * (b[1] + c[1]) 
        + (c[0] - a[0]) * (c[1] + a[1]);
}

Void UpdateBounds(const Point& pt, Point& min, Point& max)
{
    if (min[0] > pt[0])
        min[0] = pt[0];
    else if (max[0] < pt[0])
        max[0] = pt[0];

    if (min[1] > pt[1])
        min[1] = pt[1];
    else if (max[1] < pt[1])
        max[1] = pt[1];

    if (min[2] > pt[2])
        min[2] = pt[2];
    else if (max[2] < pt[2])
        max[2] = pt[2];
}

inline Double xCross(const Vector &a, const Vector &b, const Vector &c)
//  returns (a - b) x (c - a) . (1, 0, 0)
{ return((a[1] - b[1]) * (c[2] - b[2]) - (a[2] - b[2]) * (c[1] - b[1])); }
inline Double yCross(const Vector &a, const Vector &b, const Vector &c)
//  returns (a - b) x (c - a) . (0, 1, 0)
{ return((a[2] - b[2]) * (c[0] - b[0]) - (a[0] - b[0]) * (c[2] - b[2])); }
inline Double zCross(const Vector &a, const Vector &b, const Vector &c)
//  returns (a - b) x (c - a) . (0, 0, 1)
{ return((a[0] - b[0]) * (c[1] - b[1]) - (a[1] - b[1]) * (c[0] - b[0])); }

const GCLReal kEdgeFuzz = 0.0;

Bool PointIsInsideTriangle(
        const Point     &p0,      
        const Point     &p1,      
        const Point     &p2,      
        const Point     &point,
        Vector          *coords     
    )
{
    GCLReal temp, eps, a, b, c;
    Vector  normal;
    Int     normMajorAxis;
    
    // Calculate face normal & other info
    
    CalcTriAreaNormal(p0, p1, p2, normal);
    normal[0] = abs(normal[0]);
    normal[1] = abs(normal[1]);
    normal[2] = abs(normal[2]);
    normMajorAxis = MaxEltIndex(normal);            
    temp = normal[normMajorAxis];
    if (temp == 0.0)
        // degenerate triangle... can't hit this.
        return(false);
        
    eps = kEdgeFuzz * temp;
    switch(normMajorAxis)
    {
    case 0:
        a = xCross(p2, p1, point);
        if (a < eps)  return(false);

        b = xCross(p0, p2, point);
        if (b < eps)  return(false);

        c = temp - a - b;
        if (c < eps)  return(false);
        
        break;

    case 1:
        a = yCross(p2, p1, point);
        if (a < eps)  return(false);

        b = yCross(p0, p2, point);      
        if (b < eps)  return(false);

        c = temp - a - b;
        if (c < eps)  return(false);

        break;

    case 2:
        a = zCross(p2, p1, point);
        if (a < eps)  return(false);

        b = zCross(p0, p2, point);
        if (b < eps)  return(false);

        c = temp - a - b;
        if (c < eps)  return(false);
        
        break;
    }

    if (coords)
    {
        (*coords)[0] = a / temp;
        (*coords)[1] = b / temp;
        (*coords)[2] = c / temp;
    }
    
    return(true);
}

VecTrans ZViewMatrix3(Vector &dir)
{
    VecTrans M;

    M[0] = norm(MVector(dir[0], dir[1], dir[2]));
    M[1] = norm(MVector(dir[1], -dir[0], 0.0));
    M[2] = norm(MVector(-dir[0] * dir[2], -dir[1] * dir[2], 
                       sqr(dir[0]) + sqr(dir[1])));

#ifndef VL_ROW_ORIENT
    M = trans(M);
#endif

    return(M);
}

Void UnitSquareToUnitDisc(Coord in, Coord &out)
{
    GCLReal     phi, r;
    Coord       a = 2.0 * in - vl_1;    /* a is now on [-1,1]^2 */

    if (a[0] > -a[1])
    /* region 1 or 2 */
    {
        if (a[0] > a[1])
        /* region 1, also |a| > |b| */
        {
            r = a[0];
            phi = (vl_pi / 4.0) * (a[1] / a[0]);
        }
        else
        /* region 2, also |b| > |a| */
        {
            r = a[1];
            phi = (vl_pi / 4.0) * (2.0 - (a[0] / a[1]));
        }
    }
    else         
    /* region 3 or 4 */
    {
        if (a[0] < a[1]) 
        /* region 3, also |a| >= |b|, a != 0 */
        { 
            r = -a[0];
            phi = (vl_pi / 4.0) * (4.0 + (a[1] / a[0]));
        }
        else       
        /* region 4, |b| >= |a|, but a==0 and b==0 could occur. */
        {
            r = -a[1];
            if (a[1] != 0)
                phi = (vl_pi / 4.0) * (6.0 - (a[0] / a[1]));
            else
                phi = 0.0;
        }
    }

    out[0] = r * cos(phi);
    out[1] = r * sin(phi);
}

Void UnitDiskToUnitSquare(Coord in, Coord &out)
{
    GCLReal     r = len(in);
    GCLReal     phi = atan2(in[0], in[1]);
    Coord       a;

    if (phi < -vl_pi / 4.0) 
    /* in range [-pi/4,7pi/4] */
        phi += 2.0 * vl_pi;  
    if (phi < vl_pi / 4.0)
    /* region 1 */
    {  
        a[0] = r;
        a[1] = phi * a[0] / (vl_pi / 4.0);
    }
    else if (phi < 3.0 * vl_pi / 4.0)
    /* region 2 */
    {
        a[1] = r;
        a[0] = -(phi - vl_pi / 2.0) * a[1] / (vl_pi / 4.0);
    }
    else if (phi < 5.0 * vl_pi / 4.0)
    /* region 3 */
    {
        a[0] = -r;
        a[1] = (phi - vl_pi) * a[0] / (vl_pi / 4.0);
    }
    else
    /* region 4 */
    {
        a[1] = -r;
        a[0] = -(phi - 3.0 * vl_pi / 2.0) * a[1] / (vl_pi / 4.0);
    }

    out[0] = (a[0] + 1) / 2;
    out[1] = (a[1] + 1) / 2;
}

GCLReal ProjectToBox(Point &min, Point &max, Point &p, Vector &r)
{
    Vector  d;
    Int     i;

    for (i = 0; i < 3; i++)
        if (r[i] > 0.0)
            d[i] = (max[i] - p[i]) / r[i];
        else if (r[i] < 0.0)
            d[i] = (min[i] - p[i]) / r[i];
        else
            d[i] = HUGE_VAL;

    return(MinElt(d));
}

Vector RandomVector()
{
    Vector v;
    
    v[0] = vl_rand();
    v[1] = vl_rand();
    v[2] = vl_rand();

    return(v);
}

---