GammaLib

#include <iostream>

#include <string>

#include <sstream>

#include <TPad.h>

#include <TBox.h>

#include <TLine.h>

#include <TRandom.h>

#include <TGeoPolygon.h>

#include <sash/HESSArray.hh>

#include <sash/TimeDiff.hh>

#include <sash/Time.hh>

#include <stash/constants.h>

#include "../include/Region.hh"

#include "../include/Polygon.hh"

#define DEBUG 1

#include <utilities/debugging.hh>

using namespace std;

/**

 * 

 * \param vertices: list of vertex coordinates in absolute coordinates

 * (must have at least one entry). These must be defined clockwise. 

 *

 * The center of the polygon (used for e.g. bounding boxes, etc) will

 * be calculated using GetCentreOfGravity().

 */

Tools::Polygon::

Polygon(vector<Stash::Coordinate> vertices, string name) 

  : Region( vertices[0], name, "POLYGON", Stash::Lambda(0) ),

    fPoly(NULL), fX(NULL), fY(NULL)

{

  DEBUG_OUT << "Initializing vertices: " << vertices.size() << endl;

  // construct the TGeoPolygon:

  fX = new TArrayD(vertices.size());

  fY = new TArrayD(vertices.size());

  for (size_t ii=0; ii < vertices.size(); ii++) {

    (*fX)[ii] = vertices[ii].GetLambda().GetDegrees();

    (*fY)[ii] = vertices[ii].GetBeta().GetDegrees();

  }

  fPoly = new TGeoPolygon( vertices.size() );

  fPoly->SetXY( fX->GetArray(), fY->GetArray() );

  fPoly->FinishPolygon();

  // move the center to the COG:

  fCentre = GetCentreOfGravity();

}

Tools::Polygon::

Polygon( const Tools::Polygon &cpy) 

  : Region(cpy)

{

  fPoly = new TGeoPolygon( cpy.fPoly->GetNvert() );

  fX = new TArrayD( *(cpy.fX) );

  fY = new TArrayD( *(cpy.fY) );

  fPoly->SetXY(fX->GetArray(),fY->GetArray());

  fPoly->FinishPolygon();

}

Tools::Polygon::

Polygon() 

      : fPoly(NULL),fX(NULL),fY(NULL) {

  Sash::HESSArray::GetHESSArray();

}

Tools::Polygon::~Polygon() {

  delete fX;

  delete fY;

  delete fPoly;

}

double 

Tools::Polygon::GetArea() const {

  if (!fPoly) return 0;

  double centreX = fCentre.GetLambda().GetDegrees();

  double centreY = fCentre.GetBeta().GetDegrees();

  unsigned int nvert = fPoly->GetNvert();

  Double_t *px = fPoly->GetX();

  Double_t *py = fPoly->GetY();

  Double_t dx[nvert];

  Double_t dy[nvert];

  double rafactor;

  double area=0;

  for (unsigned int ii = 0; ii < nvert; ii++) {

    rafactor = cos(py[ii]*M_PI/180.0);

    dx[ii] = (px[ii] - centreX) * rafactor;

    dy[ii] = py[ii] - centreY;

  }

  for (unsigned int ii = 0; ii < nvert; ii++) {

    unsigned int jj = (ii+1)%nvert;

    area += dx[ii]*dy[jj] - dx[jj]*dy[ii];

  }

  return fabs(area)/2.0;

}

bool

Tools::Polygon::IsInside( const Stash::Coordinate& dir ) const {

  double point[2];

  point[0] = dir.GetLambda().GetDegrees();

  point[1] = dir.GetBeta().GetDegrees();

  return fPoly->Contains( point );

}

Tools::Region::VertexList 

Tools::Polygon::AsVertexList() const {

  Tools::Region::VertexList vlist;

  if (!fPoly) return vlist;

  int nvert = fPoly->GetNvert();

  Double_t *px = fPoly->GetX();

  Double_t *py = fPoly->GetY();

  vlist.resize(1);

  for (int ii=0; ii < nvert; ii++) {

    Stash::Lambda lambda(px[ii], Stash::Angle::Degrees);

    Stash::Beta beta(py[ii], Stash::Angle::Degrees);

    Stash::Coordinate cor(lambda,

                          beta,

                          GetCentre().GetSystem());

    vlist[0].push_back(InverseTransform(cor));

  }

  vlist[0].push_back( vlist[0][0] ); // close the curve

  return vlist;

}

string

Tools::Polygon::

AsFITSRegionString() const{

  if (!fPoly) return "";

  ostringstream ost;

  int nvert = fPoly->GetNvert();

  Double_t *px = fPoly->GetX();

  Double_t *py = fPoly->GetY();   

  ost << GetFITSSystemName()<<"; polygon(";

  for (int ii=0; ii < nvert; ii++) {

    Stash::Lambda lambda(px[ii], Stash::Angle::Degrees);

    Stash::Beta beta(py[ii], Stash::Angle::Degrees);

    Stash::Coordinate cor(lambda,

                          beta,

                          GetCentre().GetSystem());

    double xx = InverseTransform(cor).GetLambda().GetDegrees();

    double yy = InverseTransform(cor).GetBeta().GetDegrees();

    ost<<xx<<","<<yy;

    if (ii!= nvert-1) ost<<",";

  }  

  ost << ") # text={" << GetName() <<"} ";

  return ost.str();

}

/**

 * Grows the boundary of the region out from the center point. 

 *

 * \bug: this doesn't actually work. Need to properly get the polar coordinates.

 */

Tools::Polygon* 

Tools::Polygon::

GrowBoundary(const double degrees, const std::string name) const{

  // implement using GetOffsetOfPointFromCenter and GetPointAtOffsetFromCenter

  if (!fPoly) return NULL;

  int nvert = fPoly->GetNvert();

  Double_t *px = fPoly->GetX();

  Double_t *py = fPoly->GetY();

  vector<Stash::Coordinate> vertices;

  for (int ii=0; ii < nvert; ii++) {

    Stash::Lambda lambda(px[ii], Stash::Angle::Degrees);

    Stash::Beta beta(py[ii], Stash::Angle::Degrees);

    Stash::Coordinate cor(lambda,

                          beta,

                          fCentre.GetSystem());

    // calculate radial coordinate (R,phi) for each point centered on

    // the region center:

    Stash::DirectionVector cdv = fCentre.GetDirectionVector();

    Stash::DirectionVector pdv = cor.GetDirectionVector();

    Stash::DirectionVector ddv = pdv.RotateToZenith(cdv);

    Stash::Lambda phi( ddv.GetTheta().GetDegrees(), Stash::Angle::Degrees );

    double rad = cor.GetAngularDistance(fCentre).GetDegrees();

    // increase the radius, and create the new point:

    Stash::Beta newrad(degrees + rad, Stash::Angle::Degrees);

    DEBUG_OUT << ii  <<" PHI "<< phi.GetDegrees()

              << " RAD : " << rad << " NEW : " << newrad.GetDegrees() << endl;

    Stash::Coordinate newcor = GetPointAtRadius( newrad, phi );

    vertices.push_back(newcor);

  }

  return new Tools::Polygon( vertices, string(GetName())+name );

};

double

Tools::Polygon::

GetBound1() const {

  if (!fPoly) return 0;

  int nvert = fPoly->GetNvert();

  Double_t *px = fPoly->GetX();

  Double_t *py = fPoly->GetY();

  double maxdx = -1000;

  for (int ii=0; ii < nvert; ii++) {

    Stash::Lambda lambda(px[ii], Stash::Angle::Degrees);

    Stash::Beta beta(py[ii], Stash::Angle::Degrees);

    Stash::Coordinate cor(lambda,

                          beta,

                          GetCentre().GetSystem());

    double dx,dy;

    GetOffsetOfPointFromCenter( InverseTransform(cor), dx, dy, true );

    if (fabs(dx) > maxdx) maxdx = fabs(dx);

  }

  return 2.0*maxdx;

}

double

Tools::Polygon::

GetBound2() const {

  if (!fPoly) return 0;

  int nvert = fPoly->GetNvert();

  Double_t *px = fPoly->GetX();

  Double_t *py = fPoly->GetY();

  double maxdy = -1000;

  Stash::Coordinate cen = GetCentreOfGravity();

  for (int ii=0; ii < nvert; ii++) {

    Stash::Lambda lambda(px[ii], Stash::Angle::Degrees);

    Stash::Beta beta(py[ii], Stash::Angle::Degrees);

    Stash::Coordinate cor(lambda,

                          beta,

                          cen.GetSystem());

    double dx,dy;

    GetOffsetOfPointFromCenter( InverseTransform(cor), dx, dy,true );

    if (fabs(dy) > maxdy) maxdy = fabs(dy);

  }

  return 2.0*maxdy;

}

double

Tools::Polygon::

GetRadialBound() const {

  if (!fPoly) return 0;

  int nvert = fPoly->GetNvert();

  Double_t *px = fPoly->GetX();

  Double_t *py = fPoly->GetY();

  double maxdist = -1000;

  for (int ii=0; ii < nvert; ii++) {

    Stash::Lambda lambda(px[ii], Stash::Angle::Degrees);

    Stash::Beta beta(py[ii], Stash::Angle::Degrees);

    Stash::Coordinate cor(lambda,

                          beta,

                          GetCentre().GetSystem());

    double dist = cor.GetAngularDistance( fCentre ).GetDegrees();

    if (dist > maxdist) maxdist=dist;

  }

  return maxdist;

}

void 

Tools::Polygon::

MoveTo( Stash::Coordinate newpos ) {

  if (!fPoly) return;

  int nvert = fPoly->GetNvert();

  Double_t *px = fPoly->GetX();

  Double_t *py = fPoly->GetY();

  double cx = fCentre.GetLambda().GetDegrees();

  double cy = fCentre.GetBeta().GetDegrees();

  double nx = newpos.GetLambda().GetDegrees();

  double ny = newpos.GetBeta().GetDegrees();

  for (int ii=0; ii < nvert; ii++) {

    px[ii] += cx-nx;

    py[ii] += cy-ny;

  }

  fCentre = newpos;

}

/**

 * Custom streamer since TGeoPolygon doesn't automatically get updated when streamed...

 */

void Tools::Polygon::Streamer(TBuffer &R__b)

{

   // Stream an object of class Tools::Polygon.

   //This works around a msvc bug and should be harmless on other platforms

   typedef ::Tools::Polygon thisClass;

   UInt_t R__s, R__c;

   if (R__b.IsReading()) {

      Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v) { }

      //This works around a msvc bug and should be harmless on other platforms

      typedef Tools::Region baseClass0;

      baseClass0::Streamer(R__b);

      R__b >> fPoly;

      R__b >> fX;

      R__b >> fY;

      R__b.CheckByteCount(R__s, R__c, thisClass::IsA());

      // Initialize Polygon:

      fPoly->SetXY( fX->GetArray(), fY->GetArray() );

      fPoly->FinishPolygon();

   } else {

      R__c = R__b.WriteVersion(thisClass::IsA(), kTRUE);

      //This works around a msvc bug and should be harmless on other platforms

      typedef Tools::Region baseClass0;

      baseClass0::Streamer(R__b);

      R__b << fPoly;

      R__b << fX;

      R__b << fY;

      R__b.SetByteCount(R__c, kTRUE);

   }

}