RTVolumeRendering/Source/MathCommon.h

Go to the documentation of this file.

#ifndef MATHCOMMON_H
#define MATHCOMMON_H

#include <cmath>
#include <assert.h>

#ifndef VECTOR3_H
        #include "Vector3.h"
#endif

#ifndef POINT3_H
        #include "Point3.h"
#endif

#ifndef MATRIX_H
        #include "Matrix.h"
#endif

struct LineSegment
{
        Point3f p1; 
        Point3f p2; 
};

struct Ray
{
        Point3f p; 
        Vector3f d; 
};

struct Plane
{
        Point3f p; 
        Vector3f n; 
};

struct ScreenRect 
{
        int x1, y1; 
        int x2, y2; 
};

struct Sphere
{
        Point3f c; 
        float r; 
};

struct CubeDef
{
        Point3f c; 
        float dim; 
};

struct CubePlanes
{
        Plane plane[6];
};

struct CubeFrame
{
        LineSegment line[12];
};

template<class dataType, class timeType>
dataType Lerp(const dataType &d0, const dataType &d1, const timeType &t)
{
        return d0 + (d1 - d0) * t;
}

template<class T>
void Swap(T &v1, T &v2)
{
        T tmp = v1;
        v1 = v2;
        v2 = tmp;
}

template<class T>
T Clamp(const T &val, const T &low, const T &high)
{
        if(val < low)
                return low;
        else if(val > high)
                return high;
        return val;
}

inline CubePlanes CreatePlanes(const CubeDef &c)
{
        CubePlanes p;
        
        float d = c.dim / 2.0f;
        
        p.plane[0].n = Vector3f( 1.0f, 0.0f,  0.0f);  p.plane[0].p = Point3f( d, 0.0f, 0.0f);
        p.plane[1].n = Vector3f(-1.0f, 0.0f,  0.0f);  p.plane[1].p = Point3f(-d, 0.0f, 0.0f);
        
        p.plane[2].n = Vector3f( 0.0f,  1.0f,  0.0f); p.plane[2].p = Point3f(0.0f,  d, 0.0f);
        p.plane[3].n = Vector3f( 0.0f, -1.0f,  0.0f); p.plane[3].p = Point3f(0.0f, -d, 0.0f);
        
        p.plane[4].n = Vector3f( 0.0f,  0.0f,  1.0f); p.plane[4].p = Point3f(0.0f, 0.0f,  d);
        p.plane[5].n = Vector3f( 0.0f,  0.0f, -1.0f); p.plane[5].p = Point3f(0.0f, 0.0f, -d);
        
        return p;
}

inline CubeFrame CreateFrame(const CubeDef &c)
{
        CubeFrame f;
        
        float d = c.dim / 2.0f;
        
        Point3f dx(d, 0.0f, 0.0f);
        Point3f dy(0.0f, d, 0.0f);
        Point3f dz(0.0f, 0.0f, d);
        
        Point3f center = c.c;
        
        f.line[0].p1  = center + Point3f(-d, -d,  d); f.line[0].p2  = center + Point3f( d, -d,  d);
        f.line[1].p1  = center + Point3f(-d, -d,  d); f.line[1].p2  = center + Point3f(-d, -d, -d);
        f.line[2].p1  = center + Point3f(-d, -d,  d); f.line[2].p2  = center + Point3f(-d,  d,  d);
        
        f.line[3].p1  = center + Point3f( d, -d, -d); f.line[3].p2  = center + Point3f( d, -d,  d);
        f.line[4].p1  = center + Point3f( d, -d, -d); f.line[4].p2  = center + Point3f(-d, -d, -d);
        f.line[5].p1  = center + Point3f( d, -d, -d); f.line[5].p2  = center + Point3f( d,  d, -d);
        
        f.line[6].p1  = center + Point3f( d,  d,  d); f.line[6].p2  = center + Point3f(-d,  d,  d);
        f.line[7].p1  = center + Point3f( d,  d,  d); f.line[7].p2  = center + Point3f( d,  d, -d);
        f.line[8].p1  = center + Point3f( d,  d,  d); f.line[8].p2  = center + Point3f( d, -d,  d);
        
        f.line[9].p1  = center + Point3f(-d,  d, -d); f.line[9].p2  = center + Point3f( d,  d, -d);
        f.line[10].p1 = center + Point3f(-d,  d, -d); f.line[10].p2 = center + Point3f(-d,  d,  d);
        f.line[11].p1 = center + Point3f(-d,  d, -d); f.line[11].p2 = center + Point3f(-d, -d, -d);
        
        return f;
}

inline bool Inside(const Point3f &p, const Sphere &s)
{
        Vector3f d = p - s.c;
        return d.SquaredLength() < (s.r*s.r);
}

inline bool Inside(int x, int y, const ScreenRect &rect)
{
        return (x >= rect.x1 && x <= rect.x2 && y <= rect.y1 && y >= rect.y2);
}

inline int Intersect(const LineSegment &line, const Plane &plane, float &t)
{
        float d = Dot(plane.p, plane.n);
        float numerator = d - Dot(plane.n, line.p1);
        float denom = Dot(plane.n, line.p2 - line.p1);
        
        t = numerator / denom;
        
        return denom != 0.0f;
}

inline int Intersect(const Plane &p, const CubeFrame &c, Point3f points[6])
{
        int hits = 0;
        
        int i = 0;
        while(i < 12 && hits < 6)
        {
                float t;
                if(Intersect(c.line[i], p, t) && t >= 0.0f && t <= 1.0f)
                {
                        points[hits] = c.line[i].p1 + (c.line[i].p2 - c.line[i].p1)*t;
                        hits++;
                }
                i++;
        }
        
        return hits;
}

inline int Intersect(const LineSegment &line, const Sphere &s, Point3f points[2])
{
        const Point3f &p1 = line.p1;
        const Point3f &p2 = line.p2;
        const Point3f &p3 = s.c;
        
        Point3f pd = p2 - p1;
        
        // Tests if there is an intersection
        float un = Dot(p3 - p1, pd);
        float ud = Dot(pd, pd);
        float u = un / ud;
        
        // u must be in the range [0,1] for there to be an intersection.
        if(u <= 0 || u >= 1.0f)
                return 0;
        
        float a = ud;
        float b = 2.0f * Dot(pd, p1 - p3);
        float c = Dot(p3, p3) + Dot(p1, p1) - 2.0f * Dot(p3, p1) - s.r * s.r;
        
        float d = b * b - 4 * a * c;
        
        if(d < 0.0f)
                return 0;
        
        // the line is tangent to the sphere
        if(d == 0.0f)
        {
                float t = -b / (2.0f * a);
                
                points[0] = p1 + pd * t;
                
                return 1;
        }
        
        // two intersections
        float root = std::sqrt(d);
        float t1 = (-b - root) / (2.0f * a);
        float t2 = (-b + root) / (2.0f * a);
        
        points[0] = p1 + pd * t1;
        points[1] = p1 + pd * t2;
        
        return 2;
}

inline int Intersect(const Ray &r, const Sphere &s, Point3f points[2])
{       
        Vector3f diff = r.p - s.c;
        Point3f dir = r.d;
        
        // a == 1
        float b = 2.0f * Dot(r.d, diff);
        float c = Dot(diff, diff) - s.r * s.r;
        
        float d = b*b - 4*c;
        if(d < 0.0f)
                return 0;
        
        float root = std::sqrt(d);
        
        float t0 = (-b - root) / 2.0f;
        Point3f p0 = r.p + r.d * t0;
        
        float t1 = (-b + root) / 2.0f;
        Point3f p1 = r.p + r.d * t1;

        if(t0 < 0.0f && t1 < 0.0f)
                return 0;

        if(t0 < 0.0f || t1 < 0.0f)
        {
                if(t0 < 0.0f)
                        points[0] = Lerp(r.p, dir, t1);
                else
                        points[0] = Lerp(r.p, dir, t0);

                return 1;
        }

        points[0] = Lerp(r.p, dir, t0);
        points[1] = Lerp(r.p, dir, t1);
        return 2;
}

inline Point3f Average(const Point3f *points, int n)
{
        Point3f p(0.0f, 0.0f, 0.0f);
        for(int i = 0; i < n; i++)
        {
                p = p + points[i];
        }
        
        return p / (float)n;
}

inline float ReferanceAngle(const Vector3f &ref, const Vector3f &v1, const Vector3f &v2)
{
        Vector3f c = Cross(v1, v2);
        float angle = std::atan2(c.Length(), Dot(v1, v2));
        
        if(Dot(c, ref) < 0.0f)
                return -angle;
        return angle;
}

inline void SortPoints(Point3f *points, int n, const Vector3f &normal)
{
        assert(n <= 6);
        
        if(n < 3)
                return;
        
        Point3f ave = Average(points, n);
        Vector3f v1 = points[0] - ave;
        
        float angle[6];
        
        angle[0] = 0;
        for(int i = 1; i < n; i++)
        {
                Vector3f v2 = points[i] - ave;
                angle[i] = ReferanceAngle(normal, v1, v2);
        }
        
        // insertion sort
        for(int i = 0; i < n; i++)
        {
                float a = angle[i];
                Point3f p = points[i];
                
                int j = i - 1;
                while(j >= 0 && angle[j] > a)
                {
                        angle[j + 1] = angle[j];
                        points[j + 1] = points[j];
                        
                        j--;
                }
                
                angle[j + 1] = a;
                points[j + 1] = p;
        }
        
}

#endif

---