Adding radixmath subpackage.

cmake/Dependencies.cmake

@@ -2,6 +2,7 @@ TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
    SUBPACKAGES_DIRS_CLASSIFICATIONS_OPTREQS
      bug  radixbug      SS OPTIONAL
      io   radixio       SS OPTIONAL
+     math radixmath     SS OPTIONAL
      plot radixplot     SS OPTIONAL
      widgets radixwidgets SS OPTIONAL
 )

radixmath/CMakeLists.txt

+TRIBITS_SUBPACKAGE(math)
+
+
+SET(SOURCE
+matrix.cc
+normal.cc
+point3d.cc
+ray.cc
+util.cc
+vector3d.cc
+)
+
+SET(HEADERS
+constants.hh
+matrix.hh
+normal.hh
+point3d.hh
+ray.hh
+util.hh
+vector3d.hh
+)
+
+#
+# Add library
+TRIBITS_ADD_LIBRARY(radixmath
+  SOURCES ${SOURCE}
+)
+
+#
+# Add test directory
+TRIBITS_ADD_TEST_DIRECTORIES(tests)
+
+INSTALL(FILES ${HEADERS} DESTINATION "include/radixmath/")
+TRIBITS_SUBPACKAGE_POSTPROCESS()

radixmath/cmake/Dependencies.cmake

+TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
+LIB_REQUIRED_PACKAGES radixbug
+LIB_OPTIONAL_PACKAGES
+TEST_REQUIRED_PACKAGES testframework
+TEST_OPTIONAL_PACKAGES
+LIB_REQUIRED_TPLS
+LIB_OPTIONAL_TPLS
+TEST_REQUIRED_TPLS
+TEST_OPTIONAL_TPLS
+)

radixmath/constants.hh

+/*
+ * File:   constants.hh
+ * Author: Jordan P. Lefebvre
+ *
+ * Created on August 13, 2012, 8:26 PM
+ */
+
+#ifndef RADIX_RADIXMATH_CONSTANTS_H
+#define	RADIX_RADIXMATH_CONSTANTS_H
+#include <cmath>
+#include <limits>
+namespace radix
+{
+    /**
+     * Type def double precision
+     */
+    typedef unsigned short index_t;
+    typedef double Real;
+    const Real kEpsilon = 1.0E-12;
+    const Real lookAheadEpsilon = 1.0E-9;
+    const Real halfKEpsilon = 0.5 * kEpsilon;
+    const Real kHugeValue = 1.0E10;
+    const Real maxReal = std::numeric_limits<Real>::max();
+    const Real minReal = std::numeric_limits<Real>::min();
+    const Real PI = 2 * std::acos( Real( 0 ) );//3.1415926535897932384;
+    const Real TWO_PI = PI * 2;//6.2831853071795864769;
+    const Real PI_ON_180 = PI / 180;//0.0174532925199432957;
+    const Real invPI = 1 / PI;//0.3183098861837906715;
+    const Real invTWO_PI = 1 / TWO_PI;//0.1591549430918953358;
+    const Real infinity = std::numeric_limits<Real>::infinity();
+    const Real sqrtHALF = std::sqrt(Real(0.5));
+    const Real sqrtTWO = std::sqrt(Real(2.0));
+    const Real sqrtTHREE = std::sqrt(Real(3.0));
+    /**
+     * Type def unsigned long to Identifier
+     */
+    typedef unsigned long Identifier;
+
+    /**
+     * @function clamp
+     * @param x 1d point in space
+     * @param min
+     * @param max
+     * @return
+     */
+    inline Real clamp(const Real x, const Real min, const Real max) {
+        return (x < min ? min : (x > max ? max : x));
+    }
+    inline bool isWithin(Real x, Real y, Real within=kEpsilon){
+        return std::abs(x-y) <= within;
+    }
+    /**
+     * Determine if the two values are within kepsilon of one-another
+     * @param x
+     * @param y
+     * @return true, iff abs(x-y) <= kEpsilon
+     */
+    inline bool isWithinKEpsilon(Real x, Real y){
+        return isWithin(x,y,kEpsilon);
+    }
+
+} // namespace radix
+#endif	/* RADIX_RADIXMATH_CONSTANTS_H */
+

radixmath/matrix.cc

+#include "radixbug/bug.hh"
+#include "radixmath/matrix.hh"
+#include <sstream>
+#include <cstring>
+#include <stdexcept>
+
+namespace radix
+{
+    const int MATH_MATRIX_COLUMN = 4;
+    const int MATH_MATRIX_ROW = 4;
+
+    Matrix::Matrix() {
+        for (int x = 0; x < MATH_MATRIX_COLUMN; x++) {
+            for (int y = 0; y < MATH_MATRIX_ROW; y++) {
+                if (x == y) {
+                    matrix[x][y] = 1.0;
+                } else {
+                    matrix[x][y] = 0.0;
+                }
+            }
+        }
+    }
+
+    Matrix::Matrix(const Matrix& matrix) {
+        for (int x = 0; x < MATH_MATRIX_COLUMN; x++) {
+            for (int y = 0; y < MATH_MATRIX_ROW; y++) {
+                this->matrix[x][y] = matrix.matrix[x][y];
+            }
+        }
+    }
+
+    Matrix::~Matrix() {
+    }
+
+    bool Matrix::equals(const Matrix &m, Real eps) const
+    {
+        for (int x = 0; x < MATH_MATRIX_COLUMN; x++) {
+            for (int y = 0; y < MATH_MATRIX_ROW; y++) {
+                if( !isWithin(this->matrix[x][y], m.matrix[x][y],eps) ) return false;
+            }
+        }
+        return true;
+    }
+
+    /**
+     * @brief Copy the given matrix to this matrix
+     * @param const Matrix& rhs - the right hand side of this assignment (matrix to copy from)
+     * @return Matrix& - this matrix after the copy occurs
+     */
+    Matrix& Matrix::operator=(const Matrix& rhs) {
+        if (this == &rhs) return (*this);
+
+        for (int x = 0; x < MATH_MATRIX_COLUMN; x++) {
+            for (int y = 0; y < MATH_MATRIX_ROW; y++) {
+                matrix[x][y] = rhs.matrix[x][y];
+            }
+        }
+        return (*this);
+    }
+
+    /**
+     * @brief Multiply this matrix by the given matrix and return the result
+     * @param const Matrix& matrix - the matrix to be multiplied against this matrix
+     * @return Matrix - the result of the multiplication
+     */
+    Matrix Matrix::operator*(const Matrix& matrix) const {
+        Matrix product;
+        Real sum = 0.0;
+
+        for (int x = 0; x < MATH_MATRIX_COLUMN; x++) {
+            for (int y = 0; y < MATH_MATRIX_ROW; y++) {
+                for (int j = 0; j < MATH_MATRIX_ROW; j++) {
+                    sum += this->matrix[x][j] * matrix.matrix[j][y];
+                }
+                product.matrix[x][y] = sum;
+                sum = 0.0;
+            }
+        }
+        return product;
+    }
+
+    /**
+     * @brief Divide this matrix by the given value and return the result
+     * @param const Real value - the value to divide this matrix by
+     * @return Matrix - the matrix result from the division
+     */
+    Matrix Matrix::operator/(const Real value) {
+        radix_check(value != 0.0);
+        for (int x = 0; x < MATH_MATRIX_COLUMN; x++) {
+            for (int y = 0; y < MATH_MATRIX_ROW; y++) {
+                matrix[x][y] = matrix[x][y] / value;
+            }
+        }
+        return (*this);
+    }
+
+    /**
+     * @brief Set this matrix as an identity matrix
+     */
+    void Matrix::identity() {
+        for (int x = 0; x < MATH_MATRIX_COLUMN; x++) {
+            for (int y = 0; y < MATH_MATRIX_ROW; y++) {
+                if (x == y) {
+                    matrix[x][y] = 1.0;
+                } else {
+                    matrix[x][y] = 0.0;
+                }
+            }
+        }
+    }
+
+    Matrix& Matrix::translate(const Real dx, const Real dy, const Real dz) {
+
+        Matrix translation_matrix; // temporary translation matrix
+
+        translation_matrix.matrix[0][3] = dx;
+        translation_matrix.matrix[1][3] = dy;
+        translation_matrix.matrix[2][3] = dz;
+
+        *this = (*this) * translation_matrix;
+        return *this;
+    }//translate
+
+    Matrix& Matrix::rotateX(const Real theta) {
+
+        Real sin_theta = std::sin(theta * PI_ON_180);
+        Real cos_theta = std::cos(theta * PI_ON_180);
+
+        Matrix xRotationMatrix; // temporary rotation matrix about x axis
+
+        xRotationMatrix.matrix[1][1] = cos_theta;
+        xRotationMatrix.matrix[1][2] = -sin_theta;
+        xRotationMatrix.matrix[2][1] = sin_theta;
+        xRotationMatrix.matrix[2][2] = cos_theta;
+
+        *this = xRotationMatrix * (*this);
+        return *this;
+    }//rotateX
+
+    Matrix& Matrix::rotateY(const Real theta) {
+
+        Real sin_theta = sin(theta * PI / 180.0);
+        Real cos_theta = cos(theta * PI / 180.0);
+
+        Matrix yRotationMatrix; // temporary rotation matrix about x axis
+
+        yRotationMatrix.matrix[0][0] = cos_theta;
+        yRotationMatrix.matrix[0][2] = sin_theta;
+        yRotationMatrix.matrix[2][0] = -sin_theta;
+        yRotationMatrix.matrix[2][2] = cos_theta;
+
+        *this = yRotationMatrix * (*this);
+        return *this;
+    }//rotateY
+
+    Matrix& Matrix::rotateZ(const Real theta) {
+        Real sin_theta = sin(theta * PI / 180.0);
+        Real cos_theta = cos(theta * PI / 180.0);
+
+        Matrix zRotationMatrix; // temporary rotation matrix about y axis
+
+        zRotationMatrix.matrix[0][0] = cos_theta;
+        zRotationMatrix.matrix[0][1] = -sin_theta;
+        zRotationMatrix.matrix[1][0] = sin_theta;
+        zRotationMatrix.matrix[1][1] = cos_theta;
+
+        *this = zRotationMatrix * (*this);
+        return *this;
+    }//rotateZ
+    std::string Matrix::toString()const{
+        std::stringstream stream;
+        for( int i = 0; i < MATH_MATRIX_COLUMN; i++){
+            for(int j = 0; j < MATH_MATRIX_ROW; j++){
+                stream<<matrix[i][j]<<" ";
+            }
+            stream<<std::endl;
+        }
+        return stream.str();
+    }
+}//namespace radix

radixmath/matrix.hh

+#ifndef RADIX_RADIXMATH_MATRIX_H
+#define RADIX_RADIXMATH_MATRIX_H
+#include <string>
+#include "radixmath/constants.hh"
+namespace radix
+{
+
+    /**
+     * @brief Matrix class for affine transformations
+     */
+    class Matrix {
+        // members
+    public:
+        /**
+         * @brief actual matrix
+         */
+    Real  matrix[4][4];
+
+        // constructors / deconstructor
+    public:
+        Matrix();
+        Matrix(const Matrix& matrix);
+        ~Matrix();
+
+        // operators
+    public:
+        bool equals(const Matrix& m, Real eps=kEpsilon)const;
+
+        /**
+         * @brief Copy the given matrix to this matrix
+         * @param const Matrix& rhs - the right hand side of this assignment (matrix to copy from)
+         * @return Matrix& - this matrix after the copy occurs
+         */
+        Matrix& operator=(const Matrix& rhs);
+
+        /**
+         * @brief Multiply this matrix by the given matrix and return the result
+         * @param const Matrix& matrix - the matrix to be multiplied against this matrix
+         * @return Matrix - the result of the multiplication
+         */
+        Matrix operator*(const Matrix& matrix) const;
+
+        /**
+         * @brief Divide this matrix by the given value and return the result
+         * @param const Real value - the value to divide this matrix by
+         * @return Matrix - the matrix result from the division
+         */
+        Matrix operator/(const Real value);
+
+        /**
+         * Apply a translation to this matrix
+         * @param dx - x component of translation
+         * @param dy - y component of translation
+         * @param dz - z component of translation
+         * @return
+         */
+        Matrix& translate(const Real dx, const Real dy, const Real dz);
+        /**
+         * Apply an x rotation to this matrix
+         * @param theta - angle (degrees) to rotate
+         * @return
+         */
+        Matrix& rotateX(const Real theta);
+        /**
+         * Apply an y rotation to this matrix
+         * @param theta - angle (degrees) to rotate
+         * @return
+         */
+        Matrix& rotateY(const Real theta);
+        /**
+         * Apply an z rotation to this matrix
+         * @param theta - angle (degrees) to rotate
+         * @return
+         */
+        Matrix& rotateZ(const Real theta);
+        // general methods
+
+        /**
+         * @brief Set this matrix as an identity matrix
+         */
+        void identity();
+        std::string toString()const;
+
+    };
+} // namespace radix
+#endif

radixmath/normal.cc

+/*
+ * File:   normal.cc
+ * Author: Jordan
+ *
+ * Created on August 9, 2012, 7:01 PM
+ */
+
+#include <cmath>
+#include <sstream>
+#include <iomanip>
+
+#include "radixmath/normal.hh"
+namespace radix
+{
+
+    /**
+     * Default constructor
+     */
+    Normal::Normal()
+    : x(0.0), y(0.0), z(0.0) {
+    }
+
+    /**
+     * Constructor
+     * @param a
+     */
+    Normal::Normal(Real a)
+    : x(a), y(a), z(a) {
+    }
+
+    /**
+     * Constructor
+     * @param _x
+     * @param _y
+     * @param _z
+     */
+    Normal::Normal(Real _x, Real _y, Real _z)
+    : x(_x), y(_y), z(_z) {
+    }
+
+    /**
+     * Copy constructor
+     * @param n
+     */
+    Normal::Normal(const Normal& n)
+    : x(n.x), y(n.y), z(n.z) {
+    }
+
+    /**
+     * Copy Vector3D constructor
+     * @param v
+     */
+    Normal::Normal(const Vector3D& v)
+    : x(v.x), y(v.y), z(v.z) {
+    }
+
+    /**
+     * Destructor
+     */
+    Normal::~Normal() {
+    }
+
+    Normal& Normal::operator=(const Normal& rhs) {
+        if (this == &rhs)
+            return (*this);
+
+        x = rhs.x;
+        y = rhs.y;
+        z = rhs.z;
+
+        return (*this);
+    }
+
+    Normal& Normal::operator=(const Point3D& rhs) {
+        x = rhs.x;
+        y = rhs.y;
+        z = rhs.z;
+        return (*this);
+    }
+
+    Normal& Normal::operator=(const Vector3D& rhs) {
+        x = rhs.x;
+        y = rhs.y;
+        z = rhs.z;
+        return (*this);
+    }
+
+    void Normal::normalize(void) {
+        Real length = std::sqrt(x * x + y * y + z * z);
+        x /= length;
+        y /= length;
+        z /= length;
+    }
+
+    /**
+     * @function operator*
+     * @brief multiplication by a matrix on the left
+     *
+     * a normal is transformed by multiplying it on the left by the transpose of the upper left 3 x 3
+     * partition of the inverse transformation matrix
+     */
+    Normal operator*(const Matrix& mat, const Normal& n) {
+        return (Normal(mat.matrix[0][0] * n.x + mat.matrix[1][0] * n.y + mat.matrix[2][0] * n.z,
+                mat.matrix[0][1] * n.x + mat.matrix[1][1] * n.y + mat.matrix[2][1] * n.z,
+                mat.matrix[0][2] * n.x + mat.matrix[1][2] * n.y + mat.matrix[2][2] * n.z));
+    }
+
+    std::string Normal::toString()const
+    {
+        std::stringstream stream;
+        stream<<std::setprecision(15);
+        stream<<x<<","<<y<<","<<z;
+        return stream.str();
+    }
+}//namespace radix

radixmath/normal.hh

+/*
+ * File:   normal.hh
+ * Author: Jordan P. Lefebvre, lefebvre.jordan@gmail.com
+ *
+ * Created on August 9, 2012, 7:01 PM
+ */
+
+#ifndef RADIX_RADIXMATH_NORMAL_HH_
+#define	RADIX_RADIXMATH_NORMAL_HH_
+
+#include "radixmath/matrix.hh"
+#include "radixmath/vector3d.hh"
+#include "radixmath/point3d.hh"
+#include "radixmath/constants.hh"
+
+#include <vector>
+
+namespace radix
+{
+
+    /**
+     * @class Normal
+     */
+    class Normal {
+    public:
+        typedef std::vector<Normal> List;
+
+        Real x, y, z;
+
+    public:
+
+        /**
+         * Default constructor
+         */
+        Normal();
+        /**
+         * Constructor
+         * @param a
+         */
+        Normal(Real a);
+        /**
+         * constructor
+         * @param _x
+         * @param _y
+         * @param _z
+         */
+        Normal(Real _x, Real _y, Real _z);
+        /**
+         * Copy constructor
+         * @param n
+         */
+        Normal(const Normal& n);
+        /**
+         * Copy {@code Vector3D} constructor
+         * @param v
+         */
+        Normal(const Vector3D& v);
+
+        /**
+         * Destructor
+         */
+        ~Normal();
+
+        /**
+         * Assignment of a vector to a normal
+         * @param rhs
+         * @return
+         */
+        Normal& operator=(const Vector3D& rhs);
+
+        /**
+         * Assignment operator
+         * @param rhs
+         * @return
+         */
+        Normal& operator=(const Normal& rhs);
+
+        bool operator==(const Normal& rhs)const{
+            return isWithinKEpsilon(x,rhs.x)
+                    && isWithinKEpsilon(y,rhs.y)
+                    && isWithinKEpsilon(z,rhs.z);
+        }
+
+        /**
+         * assignment of a point to a normal
+         * @param rhs
+         * @return
+         */
+        Normal& operator=(const Point3D& rhs);
+
+        /**
+         * @function unary minus
+         * @return
+         */
+        Normal operator-(void) const;
+
+        Normal operator+(const Normal& n) const;
+
+        /**
+         * @brief compound addition
+         * @param n
+         * @return
+         */
+        Normal& operator+=(const Normal& n);
+
+        /**
+         * @brief  dot product with a vector on the right
+         * @param a
+         * @return
+         */
+        Real operator*(const Vector3D& v) const;
+
+        /**
+         * @brief multiplication by a Real on the right
+         * @param a
+         * @return
+         */
+        Normal operator*(const Real a) const;
+        /**