Adding axis aligned bounding box (aabb) and interval.

radixgeometry/CMakeLists.txt

@@ -2,11 +2,15 @@ TRIBITS_SUBPACKAGE(geometry)
 
 
 SET(SOURCE
+aabb.cc
 grid.cc
+interval.cc
 )
 
 SET(HEADERS
+aabb.hh
 grid.hh
+interval.hh
 )
 
 #

radixgeometry/aabb.cc

+/*
+ * File:   aabb.cc
+ * Author: Jordan P. Lefebvre, lefebvre.jordan@gmail.com
+ *
+ * Created on August 13, 2012, 8:23 PM
+ */
+
+#include "radixgeometry/aabb.hh"
+#include "radixgeometry/interval.hh"
+
+#include <sstream>
+#include <iomanip>
+#include <algorithm>
+namespace radix
+{
+AABB::AABB(void)
+    : x0(-1), x1(1), y0(-1), y1(1), z0(-1), z1(1) {
+}
+
+AABB::AABB(const Real _x0, const Real _x1,
+           const Real _y0, const Real _y1,
+           const Real _z0, const Real _z1)
+    : x0(_x0), x1(_x1), y0(_y0), y1(_y1), z0(_z0), z1(_z1) {
+}
+
+AABB::AABB(const Point3D p0, const Point3D p1)
+    : x0(p0.x), x1(p1.x), y0(p0.y), y1(p1.y), z0(p0.z), z1(p1.z) {
+}
+
+AABB::AABB(const AABB& aabb)
+    : x0(aabb.x0), x1(aabb.x1), y0(aabb.y0), y1(aabb.y1), z0(aabb.z0), z1(aabb.z1) {
+}
+
+AABB& AABB::operator=(const AABB& rhs) {
+    if (this == &rhs)
+        return (*this);
+
+    x0 = rhs.x0;
+    x1 = rhs.x1;
+    y0 = rhs.y0;
+    y1 = rhs.y1;
+    z0 = rhs.z0;
+    z1 = rhs.z1;
+
+    return (*this);
+}//operation=
+
+AABB::~AABB(void) {
+}
+
+bool AABB::hit(const Ray& ray, Real& t) const {
+
+    Real tx_min, ty_min, tz_min;
+    Real tx_max, ty_max, tz_max;
+    Real t_near, t_far;
+    if( ray.d.x == 1 || ray.d.x == -1){
+        if( ray.o.y < y0 || ray.o.y > y1 ) return false;
+        if( ray.o.z < z0 || ray.o.z > z1 ) return false;
+        if( ray.d.x > 0 ){
+            tx_min = (x0-ray.o.x) * ray.d.x;
+            tx_max = (x1-ray.o.x) * ray.d.x;
+        }else{
+            tx_min = (x1-ray.o.x) * ray.d.x;
+            tx_max = (x0-ray.o.x) * ray.d.x;
+        }
+        if( tx_max < kEpsilon ) return false;
+        if( tx_min <= kEpsilon ){
+            t = tx_max;
+        }else{
+            t = tx_min;
+        }
+        return true;
+    }
+    if(  ray.d.y == 1 || ray.d.y == -1){
+        if( ray.o.x < x0 || ray.o.x > x1 ) return false;
+        if( ray.o.z < z0 || ray.o.z > z1 ) return false;
+        if( ray.d.y > 0 ){
+            ty_min = (y0-ray.o.y) * ray.d.y;
+            ty_max = (y1-ray.o.y) * ray.d.y;
+        }else{
+            ty_min = (y1-ray.o.y) * ray.d.y;
+            ty_max = (y0-ray.o.y) * ray.d.y;
+        }
+        if( ty_max < kEpsilon ) return false;
+        if( ty_min <= kEpsilon ){
+            t = ty_max;
+        }else{
+            t = ty_min;
+        }
+        return true;
+    }
+    if(  ray.d.z == 1 || ray.d.z == -1 ){
+        if( ray.o.y < y0 || ray.o.y > y1 ) return false;
+        if( ray.o.x < x0 || ray.o.x > x1 ) return false;
+        if( ray.d.z > 0 ){
+            tz_min = (z0-ray.o.z) * ray.d.z;
+            tz_max = (z1-ray.o.z) * ray.d.z;
+        }else{
+            tz_min = (z1-ray.o.z) * ray.d.z;
+            tz_max = (z0-ray.o.z) * ray.d.z;
+        }
+        if( tz_max < kEpsilon ) return false;
+        if( tz_min <= kEpsilon ){
+            t = tz_max;
+        }else{
+            t = tz_min;
+        }
+
+        return true;
+    }
+    Real a = 1.0 / ray.d.x;
+    if (a >= 0) {
+        Real a1 = x0 - ray.o.x;
+        Real a2 = x1 - ray.o.x;
+        tx_min = (a1) * a;
+        tx_max = (a2) * a;
+    } else {
+        Real a1 = x0 - ray.o.x;
+        Real a2 = x1 - ray.o.x;
+        tx_min = (a2) * a;
+        tx_max = (a1) * a;
+    }
+
+    Real b = 1.0 / ray.d.y;
+    if (b >= 0) {
+        Real b1 = y0 - ray.o.y;
+        Real b2 = y1 - ray.o.y;
+        ty_min = (b1) * b;
+        ty_max = (b2) * b;
+    } else {
+        Real b1 = y0 - ray.o.y;
+        Real b2 = y1 - ray.o.y;
+        ty_min = (b2) * b;
+        ty_max = (b1) * b;
+    }
+
+    Real c = 1.0 / ray.d.z;
+    if (c >= 0) {
+        Real c1 = z0 - ray.o.z;
+        Real c2 = z1 - ray.o.z;
+        tz_min = (c1) * c;
+        tz_max = (c2) * c;
+    } else {
+        Real c1 = z0 - ray.o.z;
+        Real c2 = z1 - ray.o.z;
+        tz_min = (c2) * c;
+        tz_max = (c1) * c;
+    }
+
+
+    //int face_in, face_out;
+    // find the largest entering t value
+    if (tx_min > ty_min) {
+        t_near = tx_min;
+        //face_in = (a >= 0.0) ? 0 : 3;
+    } else {
+        t_near = ty_min;
+        //face_in = (b >= 0.0) ? 1 : 4;
+    }
+
+    if (tz_min > t_near) {
+        t_near = tz_min;
+        //face_in = (c >= 0.0) ? 2 : 5;
+    }
+
+    //find smallest exiting t value
+    if (tx_max < ty_max) {
+        t_far = tx_max;
+        //face_out = (a >= 0.0) ? 3 : 0;
+    } else {
+        t_far = ty_max;
+        //face_out = (b >= 0.0) ? 4 : 1;
+    }
+    if (tz_max < t_far) {
+        t_far = tz_max;
+        //face_out = (c >= 0.0) ? 5 : 2;
+    }
+    if (t_near < t_far && t_far > kEpsilon) {
+        if (t_near > kEpsilon) {
+            t = t_near;
+        }else{
+            t = t_far;
+        }
+        return true;
+    }
+    return false;
+}//hit
+
+bool AABB::inside(const Point3D& p) const {
+    return ((p.x > x0 && p.x < x1)
+            && (p.y > y0 && p.y < y1)
+            && (p.z > z0 && p.z < z1));
+}//inside
+bool AABB::on(const Point3D& p) const {
+    return ( ((p.x == x0 || p.x == x1)
+              && (p.y >= y0 && p.y <= y1)
+              && (p.z >= z0 && p.z <= z1))
+             || ((p.y == y0 || p.y == y1)
+                 && (p.x >= x0 && p.x <= x1)
+                 && (p.z >= z0 && p.z <= z1))
+             || ((p.z == z0 || p.z == z1)
+                 && (p.x >= x0 && p.x <= x1)
+                 && (p.y >= y0 && p.y <= y1)));
+}
+std::string AABB::toString()const{
+    std::stringstream stream;
+    stream << std::setprecision(15)
+           << "aabb x0="<<x0<<",x1="<<x1
+           << ",y0="<<y0<<",y1="<<y1
+           << ",z0="<<z0<<",z1="<<z1;
+    return stream.str();
+}
+AABB AABB::operator+(const AABB& rhs){
+    // union with a null AABB results in this
+    if( isNull(rhs) ) return *this;
+    return AABB(std::min(rhs.x0,this->x0), std::max(rhs.x1,this->x1)
+                ,std::min(rhs.y0,this->y0), std::max(rhs.y1,this->y1)
+                ,std::min(rhs.z0,this->z0), std::max(rhs.z1,this->z1));
+}
+AABB AABB::operator&(const AABB& rhs){
+    // intersection with a null AABB results in null
+    if( isNull(rhs) ) return AABB(0,0,0,0,0,0);
+    Interval x,y,z;
+    bool xb = Interval(x0,x1,0).intersect(Interval(rhs.x0,rhs.x1,1),x);
+    bool yb = Interval(y0,y1,0).intersect(Interval(rhs.y0,rhs.y1,1),y);
+    bool zb = Interval(z0,z1,0).intersect(Interval(rhs.z0,rhs.z1,1),z);
+    if( !xb || !yb || !zb ) return AABB(0,0,0,0,0,0); // no bounding box, it is a point
+    return AABB(x.a,x.b
+                ,y.a,y.b
+                ,z.a,z.b);
+}
+bool AABB::isNull(const AABB &aabb){
+    return !aabb.hasVolume();
+}
+Point3D AABB::origin() const{
+    return Point3D((x0+x1)/2.0
+                   ,(y0+y1)/2.0
+                   ,(z0+z1)/2.0
+                   );
+}
+bool AABB::intersects(const AABB& o)const
+{
+    bool xIntersects = (x0 >= o.x0 && x0 < o.x1)
+            || (x1 > o.x0 && x1 <= o.x1)
+            || (o.x0 >= x0 && o.x0 < x1)
+            || (o.x1 > x0 && o.x1 <= x1);
+
+    bool yIntersects = (y0 >= o.y0 && y0 < o.y1)
+            || (y1 > o.y0 && y1 <= o.y1)
+            || (o.y0 >= y0 && o.y0 < y1)
+            || (o.y1 > y0 && o.y1 <= y1);
+    bool zIntersects = (z0 >= o.z0 && z0 < o.z1)
+            || (z1 > o.z0 && z1 <= o.z1)
+            || (o.z0 >= z0 && o.z0 < z1)
+            || (o.z1 > z0 && o.z1 <= z1);
+    return (xIntersects && yIntersects && zIntersects);
+}
+
+AABB AABB::operator*(const Matrix & matrix)const{
+    Point3D x0y0z0= matrix * Point3D(x0,y0,z0);
+    Point3D x1y0z0= matrix * Point3D(x1,y0,z0);
+    Point3D x0y1z0= matrix * Point3D(x0,y1,z0);
+    Point3D x0y1z1= matrix * Point3D(x0,y1,z1);
+    Point3D x1y0z1= matrix * Point3D(x1,y0,z1);
+    Point3D x0y0z1= matrix * Point3D(x0,y0,z1);
+    Point3D x1y1z1= matrix * Point3D(x1,y1,z1);
+    Point3D x1y1z0= matrix * Point3D(x1,y1,z0);
+
+    AABB result;
+    result.x0 = std::min(x0y0z0.x
+                         ,std::min(x1y0z0.x
+                                   ,std::min(x0y1z0.x
+                                             ,std::min(x0y1z1.x
+                                                       ,std::min(x1y0z1.x
+                                                                 ,std::min(x0y0z1.x
+                                                                           ,std::min(x1y1z1.x,x1y1z0.x)))))));
+    if (std::abs(result.x0) == infinity) result.x0 = (result.x0 < 0) ? -maxReal : maxReal;
+
+    result.x1 = std::max(x0y0z0.x
+                         ,std::max(x1y0z0.x
+                                   ,std::max(x0y1z0.x
+                                             ,std::max(x0y1z1.x
+                                                       ,std::max(x1y0z1.x
+                                                                 ,std::max(x0y0z1.x
+                                                                           ,std::max(x1y1z1.x,x1y1z0.x)))))));
+    if (std::abs(result.x1) == infinity) result.x1 = (result.x1 < 0) ? -maxReal : maxReal;
+
+    result.y0 = std::min(x0y0z0.y
+                         ,std::min(x1y0z0.y
+                                   ,std::min(x0y1z0.y
+                                             ,std::min(x0y1z1.y
+                                                       ,std::min(x1y0z1.y
+                                                                 ,std::min(x0y0z1.y
+                                                                           ,std::min(x1y1z1.y,x1y1z0.y)))))));
+    if (std::abs(result.y0) == infinity) result.y0 = (result.y0 < 0) ? -maxReal : maxReal;
+
+    result.y1 = std::max(x0y0z0.y
+                         ,std::max(x1y0z0.y
+                                   ,std::max(x0y1z0.y
+                                             ,std::max(x0y1z1.y
+                                                       ,std::max(x1y0z1.y
+                                                                 ,std::max(x0y0z1.y
+                                                                           ,std::max(x1y1z1.y,x1y1z0.y)))))));
+    if (std::abs(result.y1) == infinity) result.y1 = (result.y1 < 0) ? -maxReal : maxReal;
+
+    result.z0 = std::min(x0y0z0.z
+                         ,std::min(x1y0z0.z
+                                   ,std::min(x0y1z0.z
+                                             ,std::min(x0y1z1.z
+                                                       ,std::min(x1y0z1.z
+                                                                 ,std::min(x0y0z1.z
+                                                                           ,std::min(x1y1z1.z,x1y1z0.z)))))));
+    if (std::abs(result.z0) == infinity) result.z0 = (result.z0 < 0) ? -maxReal : maxReal;
+
+    result.z1 = std::max(x0y0z0.z
+                         ,std::max(x1y0z0.z
+                                   ,std::max(x0y1z0.z
+                                             ,std::max(x0y1z1.z
+                                                       ,std::max(x1y0z1.z
+                                                                 ,std::max(x0y0z1.z
+                                                                           ,std::max(x1y1z1.z,x1y1z0.z)))))));
+    if (std::abs(result.z1) == infinity) result.z1 = (result.z1 < 0) ? -maxReal : maxReal;
+
+    return result;
+}
+
+std::vector<Point3D> AABB::vertices() const{
+    std::vector<Point3D> vertices;
+    vertices.push_back(Point3D(x0,y0,z0));
+    vertices.push_back(Point3D(x0,y0,z1));
+    vertices.push_back(Point3D(x0,y1,z0));
+    vertices.push_back(Point3D(x1,y0,z0));
+    vertices.push_back(Point3D(x1,y1,z0));
+    vertices.push_back(Point3D(x1,y0,z1));
+    vertices.push_back(Point3D(x1,y1,z1));
+    vertices.push_back(Point3D(x0,y1,z1));
+    return vertices;
+}
+}//namespace radix
+
+
+
+
+
+

radixgeometry/aabb.hh

+/*
+ * File:   aabb.hh
+ * Author: Jordan P. Lefebvre
+ *
+ * Created on August 13, 2012, 8:23 PM
+ */
+
+#ifndef RADIX_RADIXGEOMETRY_AABB_HH_
+#define	RADIX_RADIXGEOMETRY_AABB_HH_
+
+#include <vector>
+#include <string>
+
+#include "radixmath/point3d.hh"
+#include "radixmath/ray.hh"
+#include "radixmath/matrix.hh"
+#include "radixmath/constants.hh"
+namespace radix
+{
+
+/**
+     * @class AABB
+     * @brief Axis Aligned Bounding Box
+     */
+class AABB {
+public:
+    Real x0, x1, y0, y1, z0, z1;
+
+    /**
+         * Constructor
+         */
+    AABB(void);
+
+    /**
+         * Convenience Constructor
+         * @param x0 x origin
+         * @param x1
+         * @param y0 y origin
+         * @param y1
+         * @param z0 z origin
+         * @param z1
+         */
+    AABB(const Real x0, const Real x1,
+         const Real y0, const Real y1,
+         const Real z0, const Real z1);
+
+    /**
+         * Constructor
+         * @param p0
+         * @param p1
+         */
+    AABB(const Point3D p0, const Point3D p1);
+
+    /**
+         * Copy constructor
+         * @param bbox
+         */
+    AABB(const AABB& aabb);
+
+    /**
+         * Assignment operator
+         * @param rhs
+         * @return
+         */
+    AABB& operator=(const AABB& rhs);
+
+    AABB operator*(const Matrix & matrix)const;
+
+    bool operator==(const AABB& rhs)const{
+        return x0==rhs.x0
+                && x1 == rhs.x1
+                && y0 == rhs.y0
+                && y1 == rhs.y1
+                && z0 == rhs.z0
+                && z1 == rhs.z1;
+    }
+    /**
+         * @brief determine if this AABB has volume
+         * Volume is defined as anything but a point
+         * @return bool - true if and only if this AABB does not define a point
+         */
+    bool hasVolume() const{
+        return !(x0==x1 && y0==y1 && z0==z1);
+    }
+    bool isInfinite()const{
+        return x0 == -infinity || x1 == infinity
+                || y0 == -infinity || y1 == infinity
+                || z0 == -infinity || z1 == infinity;
+    }
+    bool isMaxReal()const{
+        return  x0 == -maxReal || x1 == maxReal
+                || y0 == -maxReal || y1 == maxReal
+                || z0 == -maxReal || z1 == maxReal;
+    }
+    static bool isNull(const AABB& aabb);
+    /**
+         * @brief Add two aabbs together
+         * This is used for calculating the bounding box of two geometries
+         * @param rhs
+         * @return AABB(std::min(rhs.x0,this->x0), std::max(rhs.x1,this->x1)
+                        ,std::min(rhs.y0,this->y0), std::max(rhs.y1,this->y1)
+                        std::min(rhs.z0,this->z0), std::max(rhs.z1,this->z1))
+         */
+    AABB operator+(const AABB& rhs);
+    /**
+         * @brief intersects determine if an intersection exists between this and the other aabb
+         * @param o - other aabb to determine intersection with
+         * @return true, only if these overlap, false otherwise
+         */
+    bool intersects(const AABB& o)const;
+    /**
+         * @brief creates the aabb bounding the intersection of the two
+         * @param rhs
+         * @return AABB
+         */
+    AABB operator&(const AABB& rhs);
+    /**
+         * Determine if the given AABB is contained inside this AABB
+         * Containment is determined if the given AABB is equal or
+         * smaller than this AABB
+         * @param a - the AABB to determine containment in this AABB
+         * @return true - iff the given AABB is contained
+         */
+    inline bool contains(const AABB & a)const{
+        return x0<=a.x0 && x1>=a.x1
+                && y0<=a.y0 && y1>=a.y1
+                && z0<=a.z0 && z1>=a.z1;
+    }
+    AABB& invert()
+    {
+        bool x0Edge;
+        bool x1Edge;
+        bool y0Edge;
+        bool y1Edge;
+        bool z0Edge;
+        bool z1Edge;
+        x0Edge = x0 == -maxReal || x0 == maxReal;
+        x1Edge = x1 == -maxReal || x1 == maxReal;
+        y0Edge = y0 == -maxReal || y0 == maxReal;
+        y1Edge = y1 == -maxReal || y1 == maxReal;
+        z0Edge = z0 == -maxReal || z0 == maxReal;
+        z1Edge = z1 == -maxReal || z1 == maxReal;
+
+        if( !x0Edge && x1Edge && y0Edge && y1Edge && z0Edge && z1Edge )
+        {
+            x1 = x0;
+            x0 = -maxReal;
+        }
+        else if( x0Edge && !x1Edge && y0Edge && y1Edge && z0Edge && z1Edge )
+        {
+            x0 = x1;
+            x1 = maxReal;
+        }
+        else if( x0Edge && x1Edge && !y0Edge && y1Edge && z0Edge && z1Edge )
+        {
+            y1 = y0;
+            y0 = -maxReal;
+        }
+        else if( x0Edge && x1Edge && y0Edge && !y1Edge && z0Edge && z1Edge )
+        {
+            y0 = y1;
+            y1 = maxReal;
+        }
+        else if( x0Edge && x1Edge && y0Edge && y1Edge && !z0Edge && z1Edge )
+        {
+            z1 = z0;
+            z0 = -maxReal;
+        }
+        else if( x0Edge && x1Edge && y0Edge && y1Edge && z0Edge && !z1Edge )
+        {
+            z0 = z1;
+            z1 = maxReal;
+        }else
+        {
+            x0=-maxReal;
+            x1= maxReal;
+            y0=-maxReal;
+            y1= maxReal;
+            z0=-maxReal;
+            z1= maxReal;
+        }
+        return *this;
+    }
+    /**
+         * Destructor
+         */
+    ~AABB(void);
+
+    /**
+         * Hit
+         * @brief Based on Kay and Kayjia method based on "slabs"
+
+         * @param ray ray to intersect
+         * @param t track length to populate on collision
+         * @return true on collision, false otherwise
+         */
+    virtual bool hit(const Ray& ray, Real& t) const;
+
+    /**
+         * checks if {@code p} is inside this aabb
+         * @param p Point3D
+         * @return true if inside, false otherwise
+         */
+    virtual bool inside(const Point3D& p) const;
+    virtual bool on(const Point3D& p)const;
+    std::string toString()const;
+    Point3D origin()const;
+    std::vector<Point3D> vertices()const;
+    Real height()const{ return std::abs(z0-z1);}
+    Real width()const{ return std::abs(x0-x1); }
+    Real depth()const{ return std::abs(y0-y1); }
+};
+
+}//namespace atlas
+
+
+#endif	/* ATLAS_AABB_H */
+

radixgeometry/grid.cc

@@ -38,6 +38,16 @@ std::vector<Real>& Grid::zplanes()
 }
 
 
+AABB Grid::bounds() const
+{
+    return mBounds;
+}
+
+void Grid::setBounds(const AABB &bounds)
+{
+    mBounds = bounds;
+}
+
 Grid::Grid()
     : mTranslation(0, 0, 0){
 }
@@ -78,6 +88,12 @@ Grid::Grid(Point3D lowerLeftCorner
     {
         mZPlanes[i] = mZPlanes[i-1] + deltaZ;
     }