GammaLib

#include <iostream>

#include <string>

#include <TPad.h>

#include <TBox.h>

#include <TLine.h>

#include <TPRegexp.h>

#include <TObjString.h>

#include <sash/HESSArray.hh>

#include <sash/TimeDiff.hh>

#include <sash/Time.hh>

#include <stash/constants.h>

#include <utilities/StringTools.hh>

#include <crash/EquatorSystem.hh>

#include <crash/RotatedSystem.hh>

#include <stdtools/strtools.hh>

#include "Region.hh"

#include "Polygon.hh"

#include "RegionTools.hh"

#define DEBUG 0

#include <utilities/debugging.hh>

using namespace std;

Tools::Region* Tools::StringToRegion (std::string input,

                                      bool ds9, bool inverted){

  Tools::Region* reg = 0;

  if (ds9) 

    reg = Tools::DS9ToRegion(input, inverted);

  else 

    reg = Tools::HESSToRegion(input);

  return reg;

}

Tools::Region* 

Tools::HESSToRegion (std::string input){

  // remove comments:

  size_t i = input.find("#");

  if (i != std::string::npos) input = input.substr(0,i);

  // Ignore empty lines (and comments)

  if(input.empty()) return 0;

  std::vector<std::string> info = split(input.c_str()," \t");

#ifdef DEBUG

  if (DEBUG > 1) {

    DEBUG_OUT_L(2) << " Read: ";

    for (unsigned int i = 0; i < info.size(); i++) 

      std::cout << info[i] << " ; ";

    std::cout << std::endl;

  }

#endif

  if (info.size() < 5) {

    throw std::runtime_error("Parse error in regions file: "+input);

  }

  // Construct region

  Tools::Region* reg = 0;

  std::string type = info[0];

  bool radec = (info[2] == "RADec" || info[2] == "RADEC");

  double xpos = atof(info[3].c_str());

  double ypos = atof(info[4].c_str());

  std::string name = info[5];

  const Stash::SystemRef *sys;

  if (radec) 

    sys = &Crash::EquatorSystem::GetRADecJ2000System();

  else

    sys = &Crash::RotatedSystem::GetGalacticSystem();

  Stash::Coordinate centre 

    = Stash::Coordinate (Stash::Lambda(xpos, Stash::Angle::Degrees),

                         Stash::Beta(ypos, Stash::Angle::Degrees), *sys);

  if (type == "BOX") {

    DEBUG_OUT_L(2) << "Found Box, columns: " << info.size() << std::endl;

    if (info.size() < 8) return 0;

    double boxX = atof(info[6].c_str());

    double boxY = atof(info[7].c_str());

    double phi = 0.0;

    if (info.size() > 8) phi = atof(info[8].c_str());

    reg = new Tools::Box(centre,boxX,boxY,

                         Stash::Lambda(phi,Stash::Angle::Degrees),name);

  } 

  else if (type == "SEGMENT") {

    DEBUG_OUT_L(2) << "Found Segment, columns: " << info.size() << std::endl;

    if (info.size() < 10) return 0;

    double r1 = atof(info[6].c_str());

    double r2 = atof(info[7].c_str());

    double phi1 = atof(info[8].c_str());

    double phi2 = atof(info[9].c_str());

    // simplify the region if the bounds are over the full range:

    if (fabs(phi1)<1e-10 && fabs(phi2-360.0) <1e-10)

      if (fabs(r1) < 1e-10)

        reg = new Tools::Circle(centre,r2,name);

      else

        reg = new Tools::Ring(centre,r1,r2,name);

    else 

      reg = new Tools::CircleSegment(centre,r1,r2,phi1,phi2,name);

  } 

  else if (type == "POLYGON" || type == "POLYCORNERS") {

    DEBUG_OUT_L(2) << "Found Polygon, columns: " << info.size() << std::endl;

    if (info.size() < 13) return 0;

    bool relative = (type == "POLYGON");

    unsigned int ptNb = atoi(info[6].c_str());

    std::vector<Stash::Coordinate> corners;

    for (unsigned int c = 7; c < 7+2*ptNb; c+=2) {

      Stash::Lambda lambda(info[c]);

      Stash::Beta beta(info[c+1]);

      Stash::Coordinate pos(lambda,

                            beta,

                            centre.GetSystem());

      if (relative)

        pos += centre; 

      corners.push_back(pos);

    }

    reg = new Tools::Polygon(corners, name); 

  } 

  else {

    WARN_OUT << "Unknown region type: " << type << std::endl;

    return 0; 

  }

  return reg;

}

/**

 * \todo: document this? what does inverted do?

 */

Tools::Region* Tools::DS9ToRegion (std::string input,

                                   bool inverted)

{

  // Ignore empty lines (and comments)

  if(input.empty()) return 0;

  if(TString(input).BeginsWith("#")) return 0;

  if(TString(input).BeginsWith("global")) return 0;

  std::vector<std::string> infos = split(input.c_str(),"={}()#;, \t");

  std::vector<std::string>::iterator info;

  if (DEBUG) {

    DEBUG_OUT << " Read: ";

    for (unsigned int i = 0; i < infos.size(); i++) 

      std::cout << infos[i] << " ; ";

    std::cout << std::endl;

  }

  // Coordinate system

  bool system = false;

  bool radec = true;

  info = infos.begin();

  while (info != infos.end()) {

    if ( *info == "fk5" || *info == "FK5") {

      system = true;

      break;

    }

    if (*info == "galactic" || *info == "GALACTIC") {

      system = true;

      radec = false;

      break;

    }

    info++;

  }

  if (system) {

    info = infos.erase(info);

  } else {

    if (radec)

      WARN_OUT << "No system found: assume RADec" << std::endl;

    else

      WARN_OUT << "No system found: assume galactic" << std::endl;

  }

  if (infos.size() == 0) return 0;

  // Region name

  std::string name = "";   

  bool text = false;

  info = infos.begin();

  while (info != infos.end()) {

    if ( *info == "text" || *info == "TEXT") {

      text = true;

      break;

    }

    info++;

  }

  if (text) {

    info = infos.erase(info);

    name = *info;

    info = infos.erase(info);

    text = true;

  }

  // Region type

  info = infos.begin();

  bool fullCS = true;

  bool disk = true;

  unsigned int paramNb = 3;

  std::string type = "SEGMENT";

  if (*info == "circle" || *info == "CIRCLE") {

    type = "SEGMENT";

    fullCS = true;

    paramNb = 3; // x,y,r

  } else if (*info == "annulus" || *info == "ANNULUS") {

    type = "SEGMENT";

    fullCS = true;

    disk = false;

    paramNb = 4; //x,y,r1,r2

  } else if (*info == "panda" || *info == "PANDA") {

    type = "SEGMENT";

    fullCS = false;

    disk = false;

    paramNb = 8; //x,y,a1,a2,na,r1,r2,nr

  } else if (*info == "box" || *info == "BOX") {

    type = "BOX";

    paramNb = 5; //x,y,d1,d2,nd

  } else if (*info == "polygon" || *info == "POLYGON") {

    //       type = "POLYGON";

    //       paramNb = (int) floor( ((int) info.size() -idx -3) / 2)+2;

    std::cout << "Polygon support not yet implemented" << std::endl;

    return 0;

  } else {

    throw std::runtime_error("Parse error in regions file: "+input);

  }

  info = infos.erase(info);

  if (infos.size() < paramNb) {

    throw std::runtime_error("Parse error in regions file: "+input);

  }

  double xpos = atof((*info).c_str());

  info = infos.erase(info);

  double ypos = atof((*info).c_str());

  info = infos.erase(info);

  if (!text) {

    name = "Reg_"+strtools::num2str(xpos)

      +"_"+strtools::num2str(ypos);

    WARN_OUT << "No name found: assume " << name << std::endl;

  }

  // Construct region    

  const Stash::SystemRef *sys;

  if (radec) 

    sys = &Crash::EquatorSystem::GetRADecJ2000System();

  else

    sys = &Crash::RotatedSystem::GetGalacticSystem();

  Stash::Coordinate centre 

    = Stash::Coordinate (Stash::Lambda(xpos, Stash::Angle::Degrees),

                         Stash::Beta(ypos, Stash::Angle::Degrees), *sys);

  DEBUG_OUT << "Found Centre : " << xpos << " , " << ypos << std::endl; 

  Tools::Region* reg = 0;

  if (type == "BOX") {

    double boxX = atof((*info).c_str());

    info = infos.erase(info);

    double boxY = atof((*info).c_str());

    info = infos.erase(info);

    double phi = 0.;

    if (infos.size() ) {

      phi = atof((*info).c_str());

      if (inverted) {

        phi -= 180.;

        phi *= -1;

      }

      info = infos.erase(info);

    }

    DEBUG_OUT << "Found Box : " << boxX 

              << " , " << boxY 

              << " , " << phi << std::endl;

    reg = new Tools::Box(centre,boxX,boxY,

                         Stash::Lambda(phi,Stash::Angle::Degrees),name);

  } else if (type == "SEGMENT") {

    double phi1 = 0.;

    double phi2 = 360.;

    int Nphi = 1;

    double r1 = 0.;

    double r2 = 0.;

    int Nr = 1;

    if (!fullCS) {

      phi1 = atof((*info).c_str());

      info = infos.erase(info);

      phi2 = atof((*info).c_str());

      info = infos.erase(info);

      Nphi = atoi((*info).c_str());

      info = infos.erase(info);

      r1   = Tools::ParseAngle((*info));

      info = infos.erase(info);

      r2   = Tools::ParseAngle((*info));

      info = infos.erase(info);

      Nr   = atoi((*info).c_str());

      info = infos.erase(info);

    } else {

      if (disk) {

        r2 = Tools::ParseAngle((*info));

        info = infos.erase(info);

      } else {

        r1   = Tools::ParseAngle((*info));

        info = infos.erase(info);

        r2   = Tools::ParseAngle((*info));

        info = infos.erase(info);

        Nr   = atoi((*info).c_str());

        info = infos.erase(info);

      }

    }

    if(!fullCS){

      if (inverted) {

        double phi = phi2;

        phi2 = 180. - phi1;

        phi1 = 180. - phi;

      }

    }

    DEBUG_OUT << "Found Segment: " 

              << phi1 << " " << phi2 << " " << Nphi << " , "

              << r1 << " " << r2 << " " << Nr 

              << std::endl;

    Tools::CircleSegment *seg = new Tools::CircleSegment

      (centre,r1,r2,phi1,phi2,name);

    reg=seg;

    // TODO: reimplement segmented regions:

    // if (Nr > 1) {

    //   Tools::Collection *sliceR = new Tools::Collection(name);

    //   sliceR->MakeSlicesFrom(seg, true, Nr);

    //   if (Nphi > 1) {

    //     Tools::Collection *slicePhi = new Tools::Collection(name);

    //     for (int r = 0; r < Nr; r++) {

    //       Tools::Collection *slice = new Tools::Collection();

    //       Tools::CircleSegment *cs 

    //         = dynamic_cast<Tools::CircleSegment*> (sliceR->GetOneRegion(r));

    //       slice->MakeSlicesFrom(cs, false, Nphi);

    //       slicePhi->AddRegion(slice);

    //       delete slice;

    //     }

    //     reg = slicePhi;

    //     delete sliceR;

    //   } else {

    //     reg = sliceR;

    //   }

    //   delete seg;

    // } else if (Nphi > 1) {

    //   Tools::Collection *slicePhi = new Tools::Collection(name);

    //   slicePhi->MakeSlicesFrom(seg, false, Nphi);

    //   delete seg;

    // } else {

    //   reg = seg;

    // }

  } 

  else if (type == "POLYGON") {

    WARN_OUT << "Polygon region not currently supported" << endl;

    return NULL;

    DEBUG_OUT << "Found Polygon: " << std::endl;

    if (infos.size() < 13) return 0;

    unsigned int ptNb = infos.size() / 2;

    std::vector<Stash::Coordinate> corners;

    for (unsigned int c = 0; c < 2*ptNb; c+=2) {

      double x = atof((*info).c_str());

      info = infos.erase(info);

      double y = atof((*info).c_str());

      info = infos.erase(info);

      Stash::Coordinate pos( Stash::Lambda(x, Stash::Angle::Degrees),

                             Stash::Beta(y, Stash::Angle::Degrees),

                             centre.GetSystem() );

      corners.push_back(pos);

    }

    reg = new Tools::Polygon (corners, name); 

  }

  return reg;

}

// Parse an angle+unit string according to ds9 format:

// ( from ds9 reference manual:

// size arguments 

//    [num]                   # context-dependent (see below)  

//    [num]"                  # arc sec 

//    [num]'                  # arc min 

//    [num]d                  # degrees 

//    [num]r                  # radians 

//    [num]p                  # physical pixels 

//    [num]i                  # image pixels )

//

//  output: angle in degrees

//

Double_t Tools::ParseAngle(std::string input){

  DEBUG_OUT << "Parsing \"" << input << "\"" << std::endl;

  if( TPRegexp("[^\\.\\d\"'drpi]").Match(input.c_str()) > 0 ){

    WARN_OUT << "ERROR: Illegal string \"" << input << "\"" << std::endl;

    return 0.0;

//     throw std::runtime_error("Parse error in ParseAngle, input: "+input);

  }

  TObjArray *subStrL = TPRegexp("([\\.\\d]*)([\"'drpi])").MatchS(input.c_str());

  for (Int_t i = 0; i < subStrL->GetLast()+1; i++) {

    const TString subStr = ((TObjString *)subStrL->At(i))->GetString();

    DEBUG_OUT_L(2) << "substr " << i << ": \"" << subStr << "\"" << std::endl;

  }

  if( subStrL->GetLast() >= 3 ){

    WARN_OUT << "ERROR parsing string \"" << input << "\"" << std::endl;

    throw std::runtime_error("Parse error in ParseAngle, input: "+input);

  }

  if( subStrL->GetLast() < 1 ){

    // input string is just a floating point number, so there is no match, but it's a perfectly valid string

    return TString(input).Atof();

  }    

  TString numberStr = ((TObjString *)subStrL->At(1))->GetString();

  Double_t number = numberStr.Atof();

  if( subStrL->GetLast() == 2 ){

    TString unitStr   = ((TObjString *)subStrL->At(2))->GetString();

    if( unitStr == "\"" ) number /= 3600.;

    else if( unitStr == "'" ) number /= 60.;

    else if( unitStr == "d" ) number /= 1.;

    else if( unitStr == "r" ) number *= TMath::RadToDeg();

    else{

      WARN_OUT << "ERROR: illegal unit string \"" << unitStr << "\" in input \"" << input << "\"" << std::endl;

      throw std::runtime_error("Parse error in ParseAngle, input: "+input);

    }

  }

  return number;

}