Elite Frontier

Back in 2002, I started writing C++ with a unique style of combining Modern C++ Design and Domain Driven Design. My endeavor was to write a 3D game engine. It was a fun exercise.

First, define the Math Trait. We will use the Real Number System

template <typename T>
struct RealNumberSystem
{
    //3 tuple Reals
    static int const N=3;    

    typedef T Radian;
    typedef T Theta;
    typedef T Rho;
    typedef T Radius;
    typedef T Scalar;
};

Next, define the Geometry Trait. We will use Standard Euclidean. The Tuple4 typedef is interesting.

struct StandardEuclidean
{
    static int const RightHanded=1;
    static int const LeftHanded=0;

    typedef Vector1x3<float> RowVector;
    typedef Vector3x1<float> ColumnVector;
    typedef Matrix3x3<float> Matrix3;

    typedef ColumnVector Origin;
    typedef ColumnVector Axis;
    typedef ColumnVector AnyVector;
    typedef ColumnVector ComponentVector;
    typedef Matrix3    Basis;

    typedef ColumnVector Point;
    typedef ColumnVector Normal;
    typedef ColumnVector Direction;

    //Four dimensions
    typedef Vector<float,4,1,Storage2Dim> Tuple4;

};

The following is the most critical class, which I think is the linch-pin that holds the math and the geometry together. Without this class, there is NO mathematical representation of a geometrical point, therefore the domain semantics within the code would eventually break down. I think this is the most conceptually beautiful code I have ever written:

//Mathematical representation of Geometrical Point

#include "Vector.h" 

template <typename T>
class Cartesian : public Vector3x1<T>
{
public:

    Cartesian(T x, T y, T z) : Vector3x1<T>(x,y,z)
    {
    }
};

Finally, define the Coordinate Systems and transformations between the systems

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//CoordinateSystem Template Class Declaration

template <typename T,
          template <typename T1> class MathTrait=RealNumberSystem,
          typename GeometryTrait=StandardEuclidean>
class Cartesian3DCoordinateSystem
{
public:
    //Direct Access
    typename GeometryTrait::Axis m_U[MathTrait<T>::N];
    typename GeometryTrait::Origin m_P;

    Cartesian3DCoordinateSystem()
    {
        if(GeometryTrait::RightHanded)
        {
            m_P.Set(0,0,0);
            m_U[0].Set(1,0,0);
            m_U[1].Set(0,1,0);
            m_U[2].Set(0,0,1);
        }
        else
        {
            m_P.Set(0,0,0);
            m_U[0].Set(0,0,1);
            m_U[1].Set(0,1,0);
            m_U[2].Set(1,0,0);
        }
    }

    //Create this Coordinate System's Basis to match X
    Cartesian3DCoordinateSystem(const typename GeometryTrait::AnyVector& X)
    {
        if(GeometryTrait::RightHanded)
        {
            m_P.Set(0,0,0);
            m_U[0].Set(1,0,0);
            m_U[1].Set(0,1,0);
            m_U[2].Set(0,0,1);
        }
        else
        {
            m_P.Set(0,0,0);
            m_U[0].Set(0,0,1);
            m_U[1].Set(0,1,0);
            m_U[2].Set(1,0,0);
        }

        typename MathTrait<T>::Scalar Y[MathTrait::N]={0};
        for (int i=0; i<MathTrait<T>::N; ++i);
        {
            Y[i]=VectorColumn3::Dot(m_U[i],(X-m_P));
            m_U[i]*=Y[i];
        }

    }

    typename GeometryTrait::Basis GetBasis()
    {
        typename GeometryTrait::Basis R(m_U[0],m_U[1],m_U[2]);
        return R;
    }
};

template <typename T,
          template <typename> class MathTrait,
          typename GeometryTrait>
Vector3x1<T> Cylindrical2Cartesian(const CylindricalCoordinateSystem<T,MathTrait,GeometryTrait>& c)
{
    //m_Angle must be in Radians
    return (Vector3x1<float>((c.m_R*sin(c.m_Angle)),(c.m_R*cos(c.m_Angle)),0)); 
}

template <typename T,
          template <typename T> class MathTrait=RealNumberSystem,
          typename GeometryTrait=StandardEuclidean>
class CylindricalCoordinateSystem
{
public:
    //Direct Access
    typename GeometryTrait::Origin m_Z;
    typename MathTrait<T>::Theta m_Angle;
    typename MathTrait<T>::Radius m_R;    //Radians

    CylindricalCoordinateSystem()
    {
        m_Angle=T();
        m_R=T();
    }

    void SetAngleDegrees(T angle)
    {
        const float pi=3.141592653f;
        const float radconv=pi/180;
        m_Angle=angle*radconv;
    }

    //The axis parameters are perpendicular
    CylindricalCoordinateSystem(T x, T y, typename GeometryTrait::Origin Z)
    {
        //x and y are relative to Z
        x=Z.x-x;
        y=Z.y=y;

        m_R=sqrt((pow(x,2)+pow(y,2)));
        m_Angle=(acos(x/m_R));
        if(y<=0) m_Angle=-m_Angle;
    }
};

template <typename T,
          template <typename T> class MathTrait,
          typename GeometryTrait>
Vector3x1<T> Spherical2Cartesian(const SphericalCoordinateSystem<T,MathTrait,GeometryTrait>& c)
{
    return (Vector3x1<float>((c.m_R*sin(c.m_PolarAngle)*cos(c.m_Azimuth)),(c.m_R*sin(c.m_PolarAngle)*cos(c.m_Azimuth)),(c.m_R*cos(c.m_PolarAngle)))); 
}

template <typename T,
          template <typename T1> class MathTrait=RealNumberSystem,
          typename GeometryTrait=StandardEuclidean>
class SphericalCoordinateSystem
{
public:
    //Direct Access
    typename MathTrait<T>::Radius m_R;
    typename MathTrait<T>::Theta m_Azimuth;
    typename MathTrait<T>::Rho m_PolarAngle;

    SphericalCoordinateSystem(typename GeometryTrait::Point P)
    {
        m_R=sqrt((pow(P.x,2)+pow(P.y,2)+pow(P.z,2)));
        m_PolarAngle=acos(P.z/m_R);
        m_Azimuth=acos(x/(sqrt(pow(P.x,2)+pow(P.y,2))));
    }
};

Define a linear component and a plane. I often wonder how things would be different, if we used: typename GeometryTrait = NonEuclidean

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//LinearComponent Template Class Definition

template <typename T,
          typename GeometryTrait=StandardEuclidean>
class LinearComponent
{
public:
    //Direct access
    typename GeometryTrait::Point m_B;
    typename GeometryTrait::ColumnVector m_M;

    LinearComponent(typename const GeometryTrait::Point& B, typename const GeometryTrait::ColumnVector& M) : 
    m_B(B), m_M(M)
    {
        m_M.Normalize();
    }
};

////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
//Plane Template Class Definition

template <typename T,
          typename GeometryTrait=StandardEuclidean>
class Plane
{

public:
    //Three dimensional representation
    typename GeometryTrait::Point m_P;
    typename GeometryTrait::Normal m_N;

    //Four dimensional representation
    typename GeometryTrait::Tuple4    m_ND;    // <N,D>

    Plane() : m_P(),m_N() 
    {

    }

    Plane(const T& px, const T& py, const T& pz,
          const T& nx, const T& ny, const T& nz, const T& d) : m_P(px,py,pz),m_N(nx,ny,nz),m_ND(nx,ny,nz,d)
    {
        m_N.Normalize();

        assert( (Vector<float,3,1>::Dot(m_N,m_P) + (m_ND.m_Storage.M[3][0])==0));
    }

    Plane(const T& A, const T& B, const T& C, const T& D) : m_N(A,B,C), m_P(), m_ND(A,B,C,D)
    {
        assert(m_N.Magnitude()>0);

        (m_ND.m_Storage.M[3][0])/=m_N.Magnitude();
        m_N.Normalize();

        //Create a dummy Point
        m_P=m_N * -m_ND.m_Storage.M[3][0];

        //D==-(N dot P)
        T ndotp=(Vector<T,3,1>::Dot(m_N,m_P));
        T d=m_ND.m_Storage.M[3][0];
        assert(fabs(d + ndotp) < 0.001f); 

        // N dot P + D == 0
        T dd=(Vector<T,3,1>::Dot(m_N,m_P) + (m_ND.m_Storage.M[3][0]));
        assert(fabs(dd) < 0.001f);    //dd==zero
    }

    void SetNormal(const T& x, const T& y, const T& z)
    {
        m_N.Set(x,y,z);

        m_N.Normalize();

        m_ND.Set(m_N.x,m_N.y,m_N.z,0);    //Direction
    }

    void SetPoint(const T& x, const T& y, const T& z)
    {
        m_P.Set(x,y,z);
    }

    void SetD(const T& d)
    {
        m_ND.m_Storage.M[3][0]=d;
    }

    T GetD() const
    {
        return m_ND.m_Storage.M[3][0];
    }

    Plane(typename const GeometryTrait::Point& P, typename const GeometryTrait::Normal& N) : m_P(P), m_N(N)
    {
        //Make calculations easier by normalizing
        m_N.Normalize();

        m_ND.Set(N.x,N.y,N.z,-(Vector<float,3,1>::Dot(m_N,P)));

        assert( (Vector<float,3,1>::Dot(m_N,m_P) + (m_ND.m_Storage.M[3][0])==0));
    }
};

I’ve left out the Matrix template and Vector template code, since it’s pretty boring.