ConvexPolygon.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
//-----------------------------------------------------------------------------
// Includes
//-----------------------------------------------------------------------------
#include "MantidGeometry/Math/ConvexPolygon.h"
#include "MantidGeometry/Math/PolygonEdge.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/V2D.h"
#include <algorithm>
#include <cfloat>
#include <sstream>
#include <utility>
 
namespace Mantid::Geometry {
 
using Kernel::V2D;
 
//-----------------------------------------------------------------------------
// Public member functions
//-----------------------------------------------------------------------------
ConvexPolygon::ConvexPolygon() : m_minX(DBL_MAX), m_maxX(-DBL_MAX), m_minY(DBL_MAX), m_maxY(-DBL_MAX), m_vertices() {}
 
ConvexPolygon::ConvexPolygon(Vertices vertices) : m_vertices(std::move(vertices)) { setup(); }
 
bool ConvexPolygon::isValid() const { return (npoints() > 2); }
 
void ConvexPolygon::clear() {
  m_vertices.clear();
  m_minX = DBL_MAX;
  m_maxX = -DBL_MAX;
  m_minY = DBL_MAX;
  m_maxY = -DBL_MAX;
}
 
void ConvexPolygon::insert(const V2D &pt) {
  m_vertices.emplace_back(pt);
  // Update extrema
  m_minX = std::min(m_minX, pt.X());
  m_maxX = std::max(m_maxX, pt.X());
  m_minY = std::min(m_minY, pt.Y());
  m_maxY = std::max(m_maxY, pt.Y());
}
 
void ConvexPolygon::insert(double x, double y) { this->insert(V2D(x, y)); }
 
const V2D &ConvexPolygon::operator[](const size_t index) const { return m_vertices[index]; }
 
const Kernel::V2D &ConvexPolygon::at(const size_t index) const {
  if (index < npoints()) {
    return m_vertices[index];
  } else {
    throw Kernel::Exception::IndexError(index, npoints(), "ConvexPolygon::at()");
  }
}
 
size_t ConvexPolygon::npoints() const { return m_vertices.size(); }
 
bool ConvexPolygon::contains(const Kernel::V2D &point) const {
  for (size_t i = 0; i < npoints(); ++i) {
    PolygonEdge edge(m_vertices[i], m_vertices[(i + 1) % npoints()]);
    if (classify(point, edge) == OnLeft) {
      return false;
    }
  }
  return true;
}
 
 
bool ConvexPolygon::contains(const ConvexPolygon &poly) const {
  // Basically just have to test if each point is inside us, this could be
  // slow
  for (size_t i = 0; i < poly.npoints(); ++i) {
    if (!this->contains(poly[i]))
      return false;
  }
  return true;
}
 
double ConvexPolygon::area() const { return 0.5 * this->determinant(); }
 
double ConvexPolygon::determinant() const {
  // Arrange the points in a Nx2 matrix where N = npoints+1
  // and each row is a vertex point with the last row equal to the first.
  // Calculate the "determinant". The matrix class needs and NxN matrix
  // as the correct definition of a determinant only exists for
  // square matrices. We could fool it by putting extra zeroes but this
  // would increase the workload for no gain
 
  // Loop over all and compute the individual elements. We assume
  // that calling next() on the vertex takes us clockwise within
  // the polygon.
  double lhs(0.0), rhs(0.0);
  for (size_t i = 0; i < npoints(); ++i) {
    const V2D &v_i = m_vertices[i];
    const V2D &v_ip1 = m_vertices[(i + 1) % npoints()];
 
    lhs += v_ip1.X() * v_i.Y();
    rhs += v_i.X() * v_ip1.Y();
  }
  return lhs - rhs;
}
 
double ConvexPolygon::minX() const { return m_minX; }
 
double ConvexPolygon::maxX() const { return m_maxX; }
 
double ConvexPolygon::minY() const { return m_minY; }
 
double ConvexPolygon::maxY() const { return m_maxY; }
 
ConvexPolygon ConvexPolygon::toPoly() const { return *this; }
 
void ConvexPolygon::setup() {
  m_minX = DBL_MAX;
  m_maxX = -DBL_MAX;
  m_minY = DBL_MAX;
  m_maxY = -DBL_MAX;
 
  for (const auto &vertex : m_vertices) {
    double x{vertex.X()}, y{vertex.Y()};
    m_minX = std::min(m_minX, x);
    m_maxX = std::max(m_maxX, x);
    m_minY = std::min(m_minY, y);
    m_maxY = std::max(m_maxY, y);
  }
}
 
double ConvexPolygon::triangleArea(const V2D &a, const V2D &b, const V2D &c) const {
  return 0.5 * (b.X() - a.X()) * (c.Y() - a.Y()) - (c.X() - a.X()) * (b.Y() - a.Y());
}
 
//-----------------------------------------------------------------------------
// Non-member non-friend functions
//-----------------------------------------------------------------------------
std::ostream &operator<<(std::ostream &os, const ConvexPolygon &polygon) {
  os << "ConvexPolygon(";
  const size_t npoints(polygon.npoints());
  for (size_t i = 0; i < npoints; ++i) {
    os << polygon[i];
    if (i < npoints - 1)
      os << ",";
  }
  os << ")";
  return os;
}
 
} // namespace Mantid::Geometry