aabb.h

Go to the documentation of this file.

/*
    This file is part of Mitsuba, a physically based rendering system.

    Copyright (c) 2007-2014 by Wenzel Jakob and others.

    Mitsuba is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License Version 3
    as published by the Free Software Foundation.

    Mitsuba is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#pragma once
#if !defined(__MITSUBA_CORE_AABB_H_)
#define __MITSUBA_CORE_AABB_H_

#include <mitsuba/core/bsphere.h>

MTS_NAMESPACE_BEGIN

/**
 * \brief Generic multi-dimensional bounding box data structure
 *
 * Maintains a component-wise minimum and maximum position and provides
 * various convenience functions to query or change them.
 *
 * \tparam T Underlying point data type (e.g. \c TPoint3<float>)
 * \ingroup libcore
 */
template <typename T> struct TAABB {
        typedef T                           PointType;
        typedef typename T::Scalar          Scalar;
        typedef typename T::VectorType      VectorType;
        typedef TRay<PointType, VectorType> RayType;

        /**
         * \brief Create a new invalid bounding box
         *
         * Initializes the components of the minimum
         * and maximum position to \f$\infty\f$ and \f$-\infty\f$,
         * respectively.
         */
        inline TAABB() {
                reset();
        }

        /// Unserialize a bounding box from a binary data stream
        inline TAABB(Stream *stream) {
                min = PointType(stream);
                max = PointType(stream);
        }

        /// Create a collapsed AABB from a single point
        inline TAABB(const PointType &p)
                : min(p), max(p) { }

        /// Create a bounding box from two positions
        inline TAABB(const PointType &min, const PointType &max)
                : min(min), max(max) {
#if defined(MTS_DEBUG)
                for (int i=0; i<PointType::dim; ++i)
                        SAssert(min[i] <= max[i]);
#endif
        }

        /// Copy constructor
        inline TAABB(const TAABB &aabb)
                : min(aabb.min), max(aabb.max) { }

        /// Equality test
        inline bool operator==(const TAABB &aabb) const {
                return min == aabb.min && max == aabb.max;
        }

        /// Inequality test
        inline bool operator!=(const TAABB &aabb) const {
                return min != aabb.min || max != aabb.max;
        }

        /// Clip to another bounding box
        inline void clip(const TAABB &aabb) {
                for (int i=0; i<PointType::dim; ++i) {
                        min[i] = std::max(min[i], aabb.min[i]);
                        max[i] = std::min(max[i], aabb.max[i]);
                }
        }

        /**
         * \brief Mark the bounding box as invalid.
         *
         * This operation sets the components of the minimum
         * and maximum position to \f$\infty\f$ and \f$-\infty\f$,
         * respectively.
         */
        inline void reset() {
                min = PointType( std::numeric_limits<Scalar>::infinity());
                max = PointType(-std::numeric_limits<Scalar>::infinity());
        }

        /// Calculate the n-dimensional volume of the bounding box
        inline Scalar getVolume() const {
                VectorType diff = max-min;
                Scalar result = diff[0];
                for (int i=1; i<PointType::dim; ++i)
                        result *= diff[i];
                return result;
        }

        /// Calculate the n-1 dimensional volume of the boundary
        inline Float getSurfaceArea() const {
                VectorType d = max - min;
                Float result = 0.0f;
                for (int i=0; i<PointType::dim; ++i) {
                        Float term = 1.0f;
                        for (int j=0; j<PointType::dim; ++j) {
                                if (i == j)
                                        continue;
                                term *= d[j];
                        }
                        result += term;
                }
                return 2.0f * result;
        }

        /// Return the center point
        inline PointType getCenter() const {
                return (max + min) * (Scalar) 0.5;
        }

        /// Return the position of one of the corners (in <tt>0..2^dim-1</tt>)
        inline PointType getCorner(int index) const {
                PointType result;
                for (int d=0; d<PointType::dim; ++d) {
                        if (index & (1 << d))
                                result[d] = max[d];
                        else
                                result[d] = min[d];
                }
                return result;
        }

        /// Return a child bounding box in a interval-, quad-, octree, etc.
        inline TAABB getChild(int index) const {
                TAABB result(getCenter());

                for (int d=0; d<PointType::dim; ++d) {
                        if (index & (1 << d))
                                result.max[d] = max[d];
                        else
                                result.min[d] = min[d];
                }

                return result;
        }

        /// Check whether a point lies on or inside the bounding box
        inline bool contains(const PointType &p) const {
                for (int i=0; i<PointType::dim; ++i)
                        if (p[i] < min[i] || p[i] > max[i])
                                return false;
                return true;
        }

        /// Check whether a given bounding box is contained within this one
        inline bool contains(const TAABB &aabb) const {
                if (!isValid())
                        return false;
                for (int i=0; i<PointType::dim; ++i)
                        if (aabb.min[i] < min[i] || aabb.max[i] > max[i])
                                return false;
                return true;
        }

        /// Axis-aligned bounding box overlap test
        inline bool overlaps(const TAABB &aabb) const {
                for (int i=0; i<PointType::dim; ++i)
                        if (max[i] < aabb.min[i] || min[i] > aabb.max[i])
                                return false;
                return true;
        }

        /// Expand the bounding box to contain another point
        inline void expandBy(const PointType &p) {
                for (int i=0; i<PointType::dim; ++i) {
                        min[i] = std::min(min[i], p[i]);
                        max[i] = std::max(max[i], p[i]);
                }
        }

        /// Expand the bounding box to contain another bounding box
        inline void expandBy(const TAABB &aabb) {
                for (int i=0; i<PointType::dim; ++i) {
                        min[i] = std::min(min[i], aabb.min[i]);
                        max[i] = std::max(max[i], aabb.max[i]);
                }
        }

        /// Calculate the squared point-AABB distance
        inline Scalar squaredDistanceTo(const PointType &p) const {
                Scalar result = 0;
                for (int i=0; i<PointType::dim; ++i) {
                        Scalar value = 0;
                        if (p[i] < min[i])
                                value = min[i] - p[i];
                        else if (p[i] > max[i])
                                value = p[i] - max[i];
                        result += value*value;
                }
                return result;
        }

        /// Calculate the point-AABB distance
        inline Scalar distanceTo(const PointType &p) const {
                return std::sqrt(squaredDistanceTo(p));
        }

        /// Calculate the minimum squared AABB-AABB distance
        inline Scalar squaredDistanceTo(const TAABB &aabb) const {
                Scalar result = 0;

                for (int i=0; i<PointType::dim; ++i) {
                        Scalar value = 0;
                        if (aabb.max[i] < min[i])
                                value = min[i] - aabb.max[i];
                        else if (aabb.min[i] > max[i])
                                value = aabb.min[i] - max[i];
                        result += value*value;
                }
                return result;
        }

        /// Calculate the minimum AABB-AABB distance
        inline Scalar distanceTo(const TAABB &aabb) const {
                return std::sqrt(squaredDistanceTo(aabb));
        }

        /// Return whether this bounding box is valid
        inline bool isValid() const {
                for (int i=0; i<PointType::dim; ++i)
                        if (max[i] < min[i])
                                return false;
                return true;
        }

        /**
         * \brief Return whether or not this bounding box
         * covers anything at all.
         *
         * A bounding box which only covers a single point
         * is considered nonempty.
         */
        inline bool isEmpty() const {
                for (int i=0; i<PointType::dim; ++i) {
                        if (max[i] > min[i])
                                return false;
                }
                return true;
        }

        /// Return the axis index with the largest associated side length
        inline int getLargestAxis() const {
                VectorType d = max - min;
                int largest = 0;

                for (int i=1; i<PointType::dim; ++i)
                        if (d[i] > d[largest])
                                largest = i;
                return largest;
        }

        /// Return the axis index with the shortest associated side length
        inline int getShortestAxis() const {
                VectorType d = max - min;
                int shortest = 0;

                for (int i=1; i<PointType::dim; ++i)
                        if (d[i] < d[shortest])
                                shortest = i;
                return shortest;
        }

        /**
         * \brief Calculate the bounding box extents
         * \return max-min
         */
        inline VectorType getExtents() const {
                return max - min;
        }

        /** \brief Calculate the near and far ray-AABB intersection
         * points (if they exist).
         *
         * The parameters \c nearT and \c farT are used to return the
         * ray distances to the intersections (including negative distances).
         * Any previously contained value is overwritten, even if there was
         * no intersection.
         *
         * \remark In the Python bindings, this function returns the
         * \c nearT and \c farT values as a tuple (or \c None, when no
         * intersection was found)
         */
        FINLINE bool rayIntersect(const RayType &ray, Float &nearT, Float &farT) const {
                nearT = -std::numeric_limits<Float>::infinity();
                farT  = std::numeric_limits<Float>::infinity();

                /* For each pair of AABB planes */
                for (int i=0; i<PointType::dim; i++) {
                        const Float origin = ray.o[i];
                        const Float minVal = min[i], maxVal = max[i];

                        if (ray.d[i] == 0) {
                                /* The ray is parallel to the planes */
                                if (origin < minVal || origin > maxVal)
                                        return false;
                        } else {
                                /* Calculate intersection distances */
                                Float t1 = (minVal - origin) * ray.dRcp[i];
                                Float t2 = (maxVal - origin) * ray.dRcp[i];

                                if (t1 > t2)
                                        std::swap(t1, t2);

                                nearT = std::max(t1, nearT);
                                farT = std::min(t2, farT);

                                if (!(nearT <= farT))
                                        return false;
                        }
                }

                return true;
        }

        /** \brief Calculate the overlap between an axis-aligned bounding box
         * and a ray segment
         *
         * This function is an extended version of the simpler \ref rayIntersect command
         * provided above. The first change is that input values passed via
         * the \c nearT and \c farT parameters are considered to specify a query interval.
         *
         * This interval is intersected against the bounding box, returning the remaining
         * interval using the \c nearT and \c farT parameters. Furthermore, the
         * interval endpoints are also returned as 3D positions via the \c near and
         * \c far parameters. Special care is taken to reduce round-off errors.
         *
         * \remark In the Python bindings, this function has the signature
         * <tt>(nearT, farT, near, far) = rayIntersect(ray, nearT, farT)</tt>.
         * It returns \c None when no intersection was found.
         */
        FINLINE bool rayIntersect(const RayType &ray, Float &nearT, Float &farT, PointType &near, PointType &far) const {
                int nearAxis = -1, farAxis = -1;

                /* For each pair of AABB planes */
                for (int i=0; i<PointType::dim; i++) {
                        const Float origin = ray.o[i];
                        const Float minVal = min[i], maxVal = max[i];

                        if (ray.d[i] == 0) {
                                /* The ray is parallel to the planes */
                                if (origin < minVal || origin > maxVal)
                                        return false;
                        } else {
                                /* Calculate intersection distances */
                                Float t1 = (minVal - origin) * ray.dRcp[i];
                                Float t2 = (maxVal - origin) * ray.dRcp[i];

                                bool flip = t1 > t2;
                                if (flip)
                                        std::swap(t1, t2);

                                if (t1 > nearT) {
                                        nearT = t1;
                                        nearAxis = flip ? (i + PointType::dim) : i;
                                }

                                if (t2 < farT) {
                                        farT = t2;
                                        farAxis = flip ? i : (i + PointType::dim);
                                }
                        }
                }

                if (!(nearT <= farT))
                        return false;

                near = ray(nearT); far = ray(farT);

                /* Avoid roundoff errors on the component where the intersection took place */
                if (nearAxis >= 0)
                        near[nearAxis % PointType::dim] = ((Float *) this)[nearAxis];

                if (farAxis >= 0)
                        far[farAxis % PointType::dim]   = ((Float *) this)[farAxis];

                return true;
        }

        /// Serialize this bounding box to a binary data stream
        inline void serialize(Stream *stream) const {
                min.serialize(stream);
                max.serialize(stream);
        }

        /// Return a string representation of the bounding box
        std::string toString() const {
                std::ostringstream oss;
                oss << "AABB" << PointType::dim << "[";
                if (!isValid()) {
                        oss << "invalid";
                } else {
                        oss << "min=" << min.toString()
                                << ", max=" << max.toString();
                }
                oss     << "]";
                return oss.str();
        }

        PointType min; ///< Component-wise minimum
        PointType max; ///< Component-wise maximum
};

/**
 * \brief Axis-aligned bounding box data structure in three dimensions
 *
 * Maintains a component-wise minimum and maximum position and provides
 * various convenience functions to query or change them.
 *
 * \ingroup libcore
 * \ingroup libpython
 */
struct AABB : public TAABB<Point> {
public:
        /**
         * \brief Create a new invalid bounding box
         *
         * Initializes the components of the minimum
         * and maximum position to \f$\infty\f$ and \f$-\infty\f$,
         * respectively.
         */
        inline AABB() : TAABB<Point>() { }

        /// Unserialize a bounding box from a binary data stream
        inline AABB(Stream *stream) : TAABB<Point>(stream) { }

        /// Create a collapsed AABB from a single point
        inline AABB(const Point &p) : TAABB<Point>(p) { }

        /// Create a bounding box from two positions
        inline AABB(const PointType &min, const PointType &max)
                : TAABB<Point>(min, max) {
        }

        /// Construct from a TAABB<Point>
        inline AABB(const TAABB<Point> &aabb)
                : TAABB<Point>(aabb) { }

        /// Calculate the surface area of the bounding box
        inline Float getSurfaceArea() const {
                Vector d = max - min;
                return (Float) 2.0 * (d.x*d.y + d.x*d.z + d.y*d.z);
        }

        /**
         * \brief Return the position of a bounding box corner
         * \param corner Requested corner index (0..7)
         */
        MTS_EXPORT_CORE Point getCorner(uint8_t corner) const;

        /**
         * \brief Bounding sphere-box overlap test
         *
         * Implements the technique proposed by Jim Arvo in
         * "A simple method for box-sphere intersection testing"
         * (Graphics Gems, 1990)
         */
        MTS_EXPORT_CORE bool overlaps(const BSphere &sphere) const;

#ifdef MTS_SSE
        /**
         * \brief Intersect against a packet of four rays.
         * \return \c false if none of the rays intersect.
         */
        FINLINE bool rayIntersectPacket(const RayPacket4 &ray, RayInterval4 &interval) const;
#endif

        /// Create a bounding sphere, which contains the axis-aligned box
        MTS_EXPORT_CORE BSphere getBSphere() const;
};

MTS_NAMESPACE_END

#endif /* __MITSUBA_CORE_AABB_H_ */