bounding_box.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2017 - 2021 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef WORLD_BUILDER_BOUNDING_BOX_H
#define WORLD_BUILDER_BOUNDING_BOX_H


#include "assert.h"
#include "point.h"

#include <limits>


namespace WorldBuilder
{


  template <unsigned int spacedim>
  class BoundingBox
  {
    public:
      BoundingBox();

      BoundingBox(const std::pair<Point<spacedim>, Point<spacedim>>
                  &boundary_points);

      template <class Container>
      BoundingBox(const Container &points);

      std::pair<Point<spacedim>, Point<spacedim>> &get_boundary_points();

      const std::pair<Point<spacedim>, Point<spacedim>> &get_boundary_points() const;

      bool
      operator==(const BoundingBox<spacedim> &box) const;

      bool
      operator!=(const BoundingBox<spacedim> &box) const;

      void
      merge_with(const BoundingBox<spacedim> &other_bbox);

      bool
      point_inside(
        const Point<spacedim> &p,
        const double tolerance = std::numeric_limits<double>::epsilon()) const;

      void
      extend(const double amount);

      double
      volume() const;

      Point<spacedim>
      center() const;

      double
      side_length(const unsigned int direction) const;

      double
      lower_bound(const unsigned int direction) const;

      double
      upper_bound(const unsigned int direction) const;

    private:
      bool
      point_inside_implementation(
        const Point<spacedim> &p,
        const double tolerance = std::numeric_limits<double>::epsilon()) const;

      std::pair<Point<spacedim>, Point<spacedim>> boundary_points;

  };


  /*------------------------ Inline functions: BoundingBox --------------------*/

#ifndef DOXYGEN


  template <>
  inline BoundingBox<2>::BoundingBox()
    :
    boundary_points ({{
      -std::numeric_limits<double>::max(),
      -std::numeric_limits<double>::max(),
      cartesian
    },
    {
      +std::numeric_limits<double>::max(),
      +std::numeric_limits<double>::max(),
      cartesian
    }
  })
  {}

  template <>
  inline BoundingBox<3>::BoundingBox()
    :
    boundary_points ({{
      -std::numeric_limits<double>::max(),
      -std::numeric_limits<double>::max(),
      -std::numeric_limits<double>::max(),
      cartesian
    },
    {
      +std::numeric_limits<double>::max(),
      +std::numeric_limits<double>::max(),
      +std::numeric_limits<double>::max(),
      cartesian
    }
  })
  {}



  template <unsigned int spacedim>
  inline BoundingBox<spacedim>::BoundingBox(
    const std::pair<Point<spacedim>, Point<spacedim>>
    &boundary_points_)
    :
    boundary_points (boundary_points_)
  {
    // We check the Bounding Box is not degenerate
    for (unsigned int i = 0; i < spacedim; ++i)
      WBAssert(boundary_points.first[i] <= boundary_points.second[i],
               "Bounding Box can't be created: the points' "
               "order should be bottom left, top right!");
  }



  template <unsigned int spacedim>
  template <class Container>
  inline BoundingBox<spacedim>::BoundingBox(const Container &points)
  {
    // Use the default constructor in case points is empty instead of setting
    // things to +oo and -oo
    if (points.size() > 0)
      {
        auto &min = boundary_points.first;
        auto &max = boundary_points.second;
        std::fill(min.begin_raw(),
                  min.end_raw(),
                  std::numeric_limits<double>::infinity());
        std::fill(max.begin_raw(),
                  max.end_raw(),
                  -std::numeric_limits<double>::infinity());

        for (const Point<spacedim> &point : points)
          for (unsigned int d = 0; d < spacedim; ++d)
            {
              min[d] = std::min(min[d], point[d]);
              max[d] = std::max(max[d], point[d]);
            }
      }
  }



  template <unsigned int spacedim>
  inline std::pair<Point<spacedim>, Point<spacedim>> &BoundingBox<spacedim>::get_boundary_points()
  {
    return this->boundary_points;
  }



  template <unsigned int spacedim>
  inline const std::pair<Point<spacedim>, Point<spacedim>> &BoundingBox<spacedim>::get_boundary_points() const
  {
    return this->boundary_points;
  }



  template <unsigned int spacedim>
  inline bool
  BoundingBox<spacedim>::
  operator==(const BoundingBox<spacedim> &box) const
  {
    return boundary_points == box.boundary_points;
  }



  template <unsigned int spacedim>
  inline bool
  BoundingBox<spacedim>::
  operator!=(const BoundingBox<spacedim> &box) const
  {
    return boundary_points != box.boundary_points;
  }



  template <unsigned int spacedim>
  inline void
  BoundingBox<spacedim>::extend(const double amount)
  {
    for (unsigned int d = 0; d < spacedim; ++d)
      {
        boundary_points.first[d] -= amount;
        boundary_points.second[d] += amount;
        WBAssert(boundary_points.first[d] <= boundary_points.second[d],
                 "Bounding Box can't be shrunk this much: the points' "
                 "order should remain bottom left, top right.");
      }
  }


  template <unsigned int spacedim>
  inline
  bool
  BoundingBox<spacedim>::point_inside(const Point<spacedim> &point,
                                      const double tolerance) const
  {
    WBAssert(boundary_points.first.get_coordinate_system() == point.get_coordinate_system(),
             "Cannot compare two points which represent different coordinate systems.");
    WBAssert(boundary_points.second.get_coordinate_system() == point.get_coordinate_system(),
             "Cannot compare two points which represent different coordinate systems.");

    if (point.get_coordinate_system() == CoordinateSystem::spherical)
      {
        Point<spacedim> other_point = point;
        if (spacedim == 2)
          {
            other_point[0] += point[0] < 0 ? 2.0 * Consts::PI : -2.0 * Consts::PI;
          }
        else
          {
            // spacedim == 3 (rad,long,lat)
            other_point[1] += point[1] < 0 ? 2.0 * Consts::PI : -2.0 * Consts::PI;
          }

        return (point_inside_implementation(point, tolerance) ||
                point_inside_implementation(other_point, tolerance));
      }

    return point_inside_implementation(point, tolerance);

  }

  template <unsigned int spacedim>
  inline
  bool
  BoundingBox<spacedim>::point_inside_implementation(const Point<spacedim> &p,
                                                     const double tolerance) const
  {
    WBAssert(boundary_points.first.get_coordinate_system() == p.get_coordinate_system(),
             "Cannot compare two points which represent different coordinate systems.");
    WBAssert(boundary_points.second.get_coordinate_system() == p.get_coordinate_system(),
             "Cannot compare two points which represent different coordinate systems.");

    for (unsigned int i = 0; i < spacedim; ++i)
      {
        // Bottom left-top right convention: the point is outside if it's smaller
        // than the first or bigger than the second boundary point The bounding
        // box is defined as a closed set
        if ((p[i] < this->boundary_points.first[i] -
             tolerance * std::abs(this->boundary_points.second[i] -
                                  this->boundary_points.first[i])) ||
            (p[i] > this->boundary_points.second[i] +
             tolerance * std::abs(this->boundary_points.second[i] -
                                  this->boundary_points.first[i])))
          return false;
      }
    return true;
  }



  template <unsigned int spacedim>
  inline
  void
  BoundingBox<spacedim>::merge_with(
    const BoundingBox<spacedim> &other_bbox)
  {
    for (unsigned int i = 0; i < spacedim; ++i)
      {
        this->boundary_points.first[i] =
          std::min(this->boundary_points.first[i],
                   other_bbox.boundary_points.first[i]);
        this->boundary_points.second[i] =
          std::max(this->boundary_points.second[i],
                   other_bbox.boundary_points.second[i]);
      }
  }



  template <unsigned int spacedim>
  double
  BoundingBox<spacedim>::volume() const
  {
    double vol = 1.0;
    for (unsigned int i = 0; i < spacedim; ++i)
      vol *= (this->boundary_points.second[i] - this->boundary_points.first[i]);
    return vol;
  }



  template <unsigned int spacedim>
  inline
  double
  BoundingBox<spacedim>::lower_bound(const unsigned int direction) const
  {
    WBAssert(direction, spacedim);

    return boundary_points.first[direction];
  }



  template <unsigned int spacedim>
  inline
  double
  BoundingBox<spacedim>::upper_bound(const unsigned int direction) const
  {
    WBAssert(direction, spacedim);

    return boundary_points.second[direction];
  }



  template <unsigned int spacedim>
  Point<spacedim>
  BoundingBox<spacedim>::center() const
  {
    Point<spacedim> point = boundary_points.first; // initialize to inherit coordinate system
    for (unsigned int i = 0; i < spacedim; ++i)
      point[i] = .5 * (boundary_points.first[i] + boundary_points.second[i]);

    return point;
  }



  template <unsigned int spacedim>
  inline
  double
  BoundingBox<spacedim>::side_length(const unsigned int direction) const
  {
    WBAssert(direction < spacedim, "Invalid index");

    return boundary_points.second[direction] - boundary_points.first[direction];
  }


#endif // DOXYGEN

}


#endif