+#ifndef MSP_GEOMETRY_BOUNDINGBOX_H_
+#define MSP_GEOMETRY_BOUNDINGBOX_H_
+
+#include <algorithm>
+#include <stdexcept>
+#include <msp/linal/vector.h>
+
+namespace Msp {
+namespace Geometry {
+
+template<typename T, unsigned D>
+class BoundingBox
+{
+private:
+ LinAl::Vector<T, D> min_pt;
+ LinAl::Vector<T, D> max_pt;
+ bool empty;
+ bool negated;
+
+public:
+ BoundingBox();
+ BoundingBox(const LinAl::Vector<T, D> &, const LinAl::Vector<T, D> &);
+
+ static BoundingBox negate(const BoundingBox &);
+
+ const LinAl::Vector<T, D> &get_minimum_point() const { return min_pt; }
+ T get_minimum_coordinate(unsigned i) const { return min_pt[i]; }
+ const LinAl::Vector<T, D> &get_maximum_point() const { return max_pt; }
+ T get_maximum_coordinate(unsigned i) const { return max_pt[i]; }
+ bool is_empty() const { return empty; }
+ bool is_negated() const { return negated; }
+};
+
+template<typename T, unsigned D>
+inline BoundingBox<T, D>::BoundingBox():
+ empty(true),
+ negated(false)
+{ }
+
+template<typename T, unsigned D>
+inline BoundingBox<T, D>::BoundingBox(const LinAl::Vector<T, D> &n, const LinAl::Vector<T, D> &x):
+ min_pt(n),
+ max_pt(x),
+ empty(false),
+ negated(false)
+{
+ for(unsigned i=0; i<D; ++i)
+ if(min_pt[i]>max_pt[i])
+ throw std::invalid_argument("BoundingBox::BoundingBox");
+}
+
+template<typename T, unsigned D>
+inline BoundingBox<T, D> BoundingBox<T, D>::negate(const BoundingBox<T, D> &bb)
+{
+ BoundingBox<T, D> result = bb;
+ result.negated = !bb.negated;
+ return result;
+}
+
+template<typename T, unsigned D>
+inline BoundingBox<T, D> operator&(const BoundingBox<T, D> &bb1, const BoundingBox<T, D> &bb2)
+{
+ if(bb1.is_empty() || bb2.is_empty())
+ return BoundingBox<T, D>();
+
+ if(bb1.is_negated())
+ {
+ if(bb2.is_negated())
+ return ~((~bb1)|(~bb2));
+ else
+ return bb2&bb1;
+ }
+
+ LinAl::Vector<T, D> result_min;
+ LinAl::Vector<T, D> result_max;
+
+ if(bb2.is_negated())
+ {
+ // This is effectively subtraction of (non-negated) bb2 from bb1
+ int uncovered_axis = -1;
+ for(unsigned i=0; i<D; ++i)
+ if(bb1.get_minimum_coordinate(i)<bb2.get_minimum_coordinate(i) || bb1.get_maximum_coordinate(i)>bb2.get_maximum_coordinate(i))
+ {
+ if(uncovered_axis!=-1)
+ return bb1;
+ uncovered_axis = i;
+ }
+
+ if(uncovered_axis==-1)
+ return BoundingBox<T, D>();
+
+ result_min = bb1.get_minimum_point();
+ result_max = bb1.get_maximum_point();
+ if(bb2.get_minimum_coordinate(uncovered_axis)<bb1.get_minimum_coordinate(uncovered_axis))
+ result_min[uncovered_axis] = bb2.get_maximum_coordinate(uncovered_axis);
+ else
+ result_max[uncovered_axis] = bb2.get_minimum_coordinate(uncovered_axis);
+ }
+ else
+ {
+ for(unsigned i=0; i<D; ++i)
+ {
+ result_min[i] = std::max(bb1.get_minimum_coordinate(i), bb1.get_minimum_coordinate(i));
+ result_max[i] = std::min(bb1.get_minimum_coordinate(i), bb1.get_minimum_coordinate(i));
+ if(result_min[i]>result_max[i])
+ return BoundingBox<T, D>();
+ }
+ }
+
+ return BoundingBox<T, D>(result_min, result_max);
+}
+
+template<typename T, unsigned D>
+inline BoundingBox<T, D> operator|(const BoundingBox<T, D> &bb1, const BoundingBox<T, D> &bb2)
+{
+ if(bb1.is_empty())
+ return bb2;
+ if(bb2.is_empty())
+ return bb1;
+
+ if(bb1.is_negated())
+ return ~((~bb1)&(~bb2));
+ else if(bb2.is_negated())
+ return ~((~bb2)&(~bb1));
+
+ LinAl::Vector<T, D> result_min;
+ LinAl::Vector<T, D> result_max;
+ for(unsigned i=0; i<D; ++i)
+ {
+ result_min[i] = std::min(bb1.get_minimum_coordinate(i), bb1.get_minimum_coordinate(i));
+ result_max[i] = std::max(bb1.get_minimum_coordinate(i), bb1.get_minimum_coordinate(i));
+ }
+
+ return BoundingBox<T, D>(result_min, result_max);
+}
+
+template<typename T, unsigned D>
+inline BoundingBox<T, D> operator~(const BoundingBox<T, D> &bb)
+{
+ return BoundingBox<T, D>::negate(bb);
+}
+
+} // namespace Geometry
+} // namespace Msp
+
+#endif