gel/vsol/vsol_polygon_2d.cxx

// This is gel/vsol/vsol_polygon_2d.cxx
#include <iostream>
#include <cmath>
#include "vsol_polygon_2d.h"
//:
// \file
#include "vsl/vsl_vector_io.h"
#include <cassert>
#ifdef _MSC_VER
#  include "vcl_msvc_warnings.h"
#endif
#include <vsol/vsol_point_2d.h>
#include "vgl/vgl_vector_2d.h"

//***************************************************************************
// Initialization
//***************************************************************************

//---------------------------------------------------------------------------
//: Constructor from a std::vector (not a geometric vector but a list of points)
// Require: new_vertices.size()>=3
//---------------------------------------------------------------------------
vsol_polygon_2d::vsol_polygon_2d(const std::vector<vsol_point_2d_sptr> &new_vertices)
  : vsol_region_2d()
{
  // require
  assert(new_vertices.size()>=3);

  storage_=new std::vector<vsol_point_2d_sptr>(new_vertices);
}

//---------------------------------------------------------------------------
// Copy constructor
//---------------------------------------------------------------------------
vsol_polygon_2d::vsol_polygon_2d(const vsol_polygon_2d &other)
  : vsol_region_2d(other)
{
  //vsol_point_2d_sptr p;
  storage_=new std::vector<vsol_point_2d_sptr>(*other.storage_);
  for (unsigned int i=0;i<storage_->size();++i)
    (*storage_)[i]=new vsol_point_2d(*((*other.storage_)[i]));
}

//---------------------------------------------------------------------------
// Destructor
//---------------------------------------------------------------------------
vsol_polygon_2d::~vsol_polygon_2d()
{
  for (auto & i : *storage_)
   i = nullptr;
  delete storage_;
}

//---------------------------------------------------------------------------
//: Clone `this': creation of a new object and initialization
// See Prototype pattern
//---------------------------------------------------------------------------
vsol_spatial_object_2d *vsol_polygon_2d::clone() const {
  return new vsol_polygon_2d(*this);
}

//***************************************************************************
// Safe casting
//***************************************************************************

vsol_polygon_2d *vsol_polygon_2d::cast_to_polygon() {
  if (!cast_to_triangle()||!cast_to_rectangle())
    return this;
  else
    return nullptr;
}

const vsol_polygon_2d *vsol_polygon_2d::cast_to_polygon() const {
  if (!cast_to_triangle()||!cast_to_rectangle())
    return this;
  else
    return nullptr;
}

vsol_triangle_2d *vsol_polygon_2d::cast_to_triangle() { return nullptr; }
const vsol_triangle_2d *vsol_polygon_2d::cast_to_triangle() const {
  return nullptr;
}

vsol_rectangle_2d *vsol_polygon_2d::cast_to_rectangle() { return nullptr; }
const vsol_rectangle_2d *vsol_polygon_2d::cast_to_rectangle() const {
  return nullptr;
}

//***************************************************************************
// Access
//***************************************************************************

//---------------------------------------------------------------------------
//: Return vertex `i'
// Require: valid_index(i)
//---------------------------------------------------------------------------
vsol_point_2d_sptr vsol_polygon_2d::vertex(const int i) const
{
  // require
  assert(valid_index(i));

  return (*storage_)[i];
}

//***************************************************************************
// Comparison
//***************************************************************************

//---------------------------------------------------------------------------
//: Has `this' the same points than `other' in the same order ?
//---------------------------------------------------------------------------
bool vsol_polygon_2d::operator==(const vsol_polygon_2d &other) const
{
  bool result = (this==&other);

  if (!result)
  {
    result = (storage_->size()==other.storage_->size());
    if (result)
    {
      vsol_point_2d_sptr p=(*storage_)[0];

      unsigned int i=0;
      for (result=false;i<storage_->size()&&!result;++i)
        result = (*p==*(*other.storage_)[i]);
      if (result)
      {
        for (unsigned int j=1;j<size()&&result;++i,++j)
        {
          if (i>=storage_->size()) i=0;
          result = ((*storage_)[i]==(*storage_)[j]);
        }
      }
    }
  }
  return result;
}

//: spatial object equality

bool vsol_polygon_2d::operator==(const vsol_spatial_object_2d& obj) const
{
  return
    obj.cast_to_region() && obj.cast_to_region()->cast_to_polygon() &&
    *this == *obj.cast_to_region()->cast_to_polygon();
}


//---------------------------------------------------------------------------
//: Compute the bounding box of `this'
//---------------------------------------------------------------------------
void vsol_polygon_2d::compute_bounding_box() const {
  set_bounding_box((*storage_)[0]->x(), (*storage_)[0]->y());
  for (unsigned int i=1;i<storage_->size();++i)
    add_to_bounding_box((*storage_)[i]->x(), (*storage_)[i]->y());
}

//---------------------------------------------------------------------------
//: Return the area of `this'
//---------------------------------------------------------------------------
double vsol_polygon_2d::area() const {
  double area = 0.0;
  unsigned int last = storage_->size()-1;

  for (unsigned int i=0; i<last; ++i)
    area += ((*storage_)[i]->x() * (*storage_)[i+1]->y())
          - ((*storage_)[i+1]->x() * (*storage_)[i]->y());

  area += ((*storage_)[last]->x() * (*storage_)[0]->y())
        - ((*storage_)[0]->x() * (*storage_)[last]->y());

  return std::abs(area / 2.0);
}

//---------------------------------------------------------------------------
//: Return the centroid of `this'
//---------------------------------------------------------------------------
// The centroid is computed by using Green's theorem to compute the
// area-weighted 1st moments of the polygon.
//  Green's theorem relates the surface integral to the line integral around
//  the boundary as:
//     Int(surface) x dxdy = 0.5 * Int(boundary) x*x dy
//     Int(surface) y dxdy = 0.5 * Int(boundary) y*y dx
//  The centroid is given by
//     xc = Int(surface) x dxdy / Int(surface) dxdy  = Int(surface) x dxdy/area
//     yc = Int(surface) y dxdy / Int(surface) dxdy  = Int(surface) y dxdy/area
//
//  For a polygon: with vertices x[i], y[i]
//   0.5 * Int(boundary) x*x dy =
//   1/6 * Sum(i)( x[i+1] + x[i] ) * ( x[i] * y[i+1] - x[i+1] * y[i] )
//
//   0.5 * Int(boundary) y*y dx =
//   1/6 * Sum(i)( y[i+1] + y[i] ) * ( x[i] * y[i+1] - x[i+1] * y[i] )
//
//  In the case of degenerate polygons, where area == 0, return the average of
//  the vertex locations.
//
vsol_point_2d_sptr vsol_polygon_2d::centroid() const {
  unsigned int nverts = storage_->size();
  assert(nverts>0);
  double sx = 0, sy = 0;
  double area = 0;
   //fill in edge from last point to first point
  vsol_point_2d_sptr pi = (*storage_)[nverts-1];
  vsol_point_2d_sptr pi1 = (*storage_)[0];
  double xi = pi->x(), yi = pi->y();
  double xi1 = pi1->x(), yi1 = pi1->y();

  //for degenerate polygons
  double xsum = xi, ysum = yi;

  double temp = xi*yi1 - xi1*yi;
  area += temp;
  sx += temp*(xi1 + xi);
  sy += temp*(yi1 + yi);
  for (unsigned int i=0; i<nverts-1; ++i)
    {
      pi = (*storage_)[i];
      pi1 = (*storage_)[i+1];
      xi = pi->x(), yi = pi->y();
      xi1 = pi1->x(), yi1 = pi1->y();
      temp = xi*yi1 - xi1*yi;
      area += temp;
      sx += temp*(xi1 + xi);
      sy += temp*(yi1 + yi);
      xsum += xi; ysum += yi;
    }
  area /= 2.0;
  if(std::fabs(area)!=0.0)
    {
      double xc = sx/(6.0*area), yc = sy/(6.0*area);
      return new vsol_point_2d(xc, yc);
    }
  return new vsol_point_2d(xsum/nverts, ysum/nverts);
}

//---------------------------------------------------------------------------
//: Is `this' convex ?
// A polygon is convex if the direction of "turning" at every vertex is
// the same.  This is checked by calculating the cross product of two
// consecutive edges and verifying that these all have the same sign.
//---------------------------------------------------------------------------
bool vsol_polygon_2d::is_convex() const {
  if (storage_->size() == 3)
    return true; // A triangle is always convex

  // First find a non-zero cross product.  This is certainly present,
  // unless the polygon collapses to a line segment.
  // Note that cross-product=0 means that two edges are parallel, which
  // is perfectly valid, but the other "turnings" should still all be in
  // the same direction.  An earlier implementation allowed for turning
  // in the other direction after a cross-product=0.

  double n = 0.0;
  for (unsigned int i = 0; i < storage_->size(); ++i) {
    int j = (i > 1) ? i - 2 : i - 2 + storage_->size();
    int k = (i > 0) ? i - 1 : i - 1 + storage_->size();
    vsol_point_2d_sptr p0 = (*storage_)[k];
    vsol_point_2d_sptr p1 = (*storage_)[j];
    vsol_point_2d_sptr p2 = (*storage_)[i];
    vgl_vector_2d<double> v1 = p0->to_vector(*p1);
    vgl_vector_2d<double> v2 = p1->to_vector(*p2);
    n = cross_product(v1, v2);
    if (n != 0.0)
      break;
  }
  if (n == 0.0)
    return true;

  for (unsigned int i = 0; i < storage_->size(); ++i) {
    int j = (i > 1) ? i - 2 : i - 2 + storage_->size();
    int k = (i > 0) ? i - 1 : i - 1 + storage_->size();
    vsol_point_2d_sptr p0 = (*storage_)[k];
    vsol_point_2d_sptr p1 = (*storage_)[j];
    vsol_point_2d_sptr p2 = (*storage_)[i];
    vgl_vector_2d<double> v1 = p0->to_vector(*p1);
    vgl_vector_2d<double> v2 = p1->to_vector(*p2);
    double n2 = cross_product(v1, v2);
    if (n2 * n < 0)
      return false; // turns in the other direction
  }
  return true;
}

//----------------------------------------------------------------
// ================   Binary I/O Methods ========================
//----------------------------------------------------------------

//: Binary save self to stream.
void vsol_polygon_2d::b_write(vsl_b_ostream &os) const
{
  vsl_b_write(os, version());
  vsol_spatial_object_2d::b_write(os);
  if (!storage_)
    vsl_b_write(os, false); // Indicate null pointer stored
  else
  {
    vsl_b_write(os, true); // Indicate non-null pointer stored
    vsl_b_write(os, *storage_);
  }
}

//: Binary load self from stream (not typically used)
void vsol_polygon_2d::b_read(vsl_b_istream &is)
{
  short ver;
  vsl_b_read(is, ver);
  switch (ver)
  {
   case 1:
    vsol_spatial_object_2d::b_read(is);

    delete storage_;
    storage_ = new std::vector<vsol_point_2d_sptr>();
    bool null_ptr;
    vsl_b_read(is, null_ptr);
    if (!null_ptr)
      return;
    vsl_b_read(is, *storage_);
    break;
   default:
    std::cerr << "vsol_polygon_2d: unknown I/O version " << ver << '\n';
  }
}

//: Return IO version number;
short vsol_polygon_2d::version() const
{
  return 1;
}

//: Print an ascii summary to the stream
void vsol_polygon_2d::print_summary(std::ostream &os) const
{
  os << *this;
}

//***************************************************************************
// Implementation
//***************************************************************************

//---------------------------------------------------------------------------
//: Default constructor.
//---------------------------------------------------------------------------
vsol_polygon_2d::vsol_polygon_2d() {
  storage_=new std::vector<vsol_point_2d_sptr>();
}

bool vsol_polygon_2d::valid_vertices(const std::vector<vsol_point_2d_sptr> ) const
{
  return true;
}


//: Binary save vsol_polygon_2d to stream.
void
vsl_b_write(vsl_b_ostream &os, const vsol_polygon_2d* p)
{
  if (p==nullptr) {
    vsl_b_write(os, false); // Indicate null pointer stored
  }
  else{
    vsl_b_write(os,true); // Indicate non-null pointer stored
    p->b_write(os);
  }
}


//: Binary load vsol_polygon_2d from stream.
void
vsl_b_read(vsl_b_istream &is, vsol_polygon_2d* &p)
{
  delete p;
  bool not_null_ptr;
  vsl_b_read(is, not_null_ptr);
  if (not_null_ptr) {
    p = new vsol_polygon_2d();
    p->b_read(is);
  }
  else
    p = nullptr;
}


inline void vsol_polygon_2d::describe(std::ostream &strm, int blanking) const
{
  if (blanking < 0) blanking = 0; while (blanking--) strm << ' ';
  if (size() == 0)
    strm << "[null]";
  else {
    strm << "[Nverts=" << size()
         << " Area=" << area();
    for (unsigned int i=0; i<size(); ++i)
      strm << " p" << i << ':' << *(vertex(i));
    strm << ']';
  }
  strm << std::endl;
}