JConvexHull2D.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
 
#include "TROOT.h"
#include "TRandom3.h"
#include "TH2D.h"
#include "TGraph.h"
#include "TMarker.h"
#include "TEllipse.h"
#include "TCanvas.h"
#include "TApplication.h"
 
#include "JGeometry2D/JVector2D.hh"
#include "JGeometry2D/JCircle2D.hh"
#include "JGeometry2D/JGeometry2DToolkit.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Example program to test convex hull and enclosing circle.
 * \author mdejong
 */
int main(int argc, char* argv[])
{
  using namespace std;
 
  string  inputFile;
  string  outputFile;
  int     numberOfPoints;
  ULong_t seed;
  int     debug;
 
  try {
 
    JParser<> zap("Example program to test convex hull and enclosing circle.");
 
    zap['f'] = make_field(inputFile)           = "";
    zap['o'] = make_field(outputFile)          = "";
    zap['n'] = make_field(numberOfPoints)      =  0;
    zap['S'] = make_field(seed)                =  0;
    zap['d'] = make_field(debug)               =  0;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  gRandom->SetSeed(seed);
 
  using namespace JPP;
 
 
  vector<JVector2D> buffer;
 
  typedef vector<JVector2D>::const_iterator const_iterator;
  typedef vector<JVector2D>::iterator       iterator;
 
 
  if (inputFile != "") {
 
    ifstream in(inputFile.c_str());
 
    for (double x, y; in >> x >> y; ) {
      buffer.push_back(JVector2D(x,y));
    }
 
    in.close();
 
  } else if (numberOfPoints > 0) {
 
    NOTICE("Seed: " << gRandom->GetSeed() << endl);
 
    for (int i = 0; i != numberOfPoints; ++i) {
 
      buffer.push_back(JVector2D(gRandom->Uniform(-1.0, +1.0),
                                 gRandom->Uniform(-1.0, +1.0)));
    }
 
    if (outputFile != "") {
      
      ofstream out(outputFile.c_str(), ios::out);
 
      for (const_iterator i = buffer.begin(); i != buffer.end(); ++i)
        out << FIXED(7,3) << i->getX() << ' ' 
            << FIXED(7,3) << i->getY() << endl;
      
      out.close();
    }
 
  } else {
 
    FATAL("No points to draw." << endl);
  }
 
  for (const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
    DEBUG(FIXED(7,3) << i->getX() << ' ' << FIXED(7,3) << i->getY() << endl);
  }
 
  const JCircle2D circle(buffer.begin(), buffer.end());
 
  pair<iterator,iterator> hull = getConvexHull2D(buffer.begin(), 
                                                 buffer.end());
 
 
  DEBUG("upper: " 
        << 0                                    << ' ' 
        << distance(buffer.begin(),hull.first)  << endl);
 
  DEBUG("lower: " 
        << distance(buffer.begin(),hull.first)  << ' '
        << distance(buffer.begin(),hull.second) << endl);
 
  DEBUG("circle: " << circle.getX() << ' ' << circle.getY() << ' ' << circle.getRadius() << endl);
 
  double area = getArea2D(buffer.begin(), hull.second);
  
  NOTICE("Area: " << area << endl);
 
 
  { 
    int n[] = { 0, 0 };
 
    for (const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
      if (!inside2D(buffer.begin(), hull.second,             *i)) { ++n[0]; }
      if (!inside2D(buffer.begin(), hull.first, hull.second, *i)) { ++n[1]; }
 
      if (!inside2D(buffer.begin(), hull.first, hull.second, *i)) { 
        
      }
    }
    
    for (int i = 0; i != sizeof(n)/sizeof(n[0]); ++i) {
      NOTICE("Number of points outside [" << i << "] " << n[i] << endl);
    }
  }
 
  TApplication* tp = new TApplication("user", NULL, NULL);
 
  TCanvas cv("cv", "", 400, 400);
 
  cv.SetFillStyle(4000);
  cv.SetFillColor(0);
 
  cv.Divide(1, 1);
  cv.cd(1);
 
 
  TEllipse  ellipse(circle.getX(), circle.getY(), circle.getRadius());
 
 
  const Int_t MAX_BUFFER_SIZE = buffer.size() + 1;
 
  Double_t  x[MAX_BUFFER_SIZE];
  Double_t  y[MAX_BUFFER_SIZE];
 
  int N = 0;
 
  for (const_iterator i = buffer.begin(); i != buffer.end(); ++i, ++N) {
    x[N] = i->getX(); 
    y[N] = i->getY(); 
  }
 
 
  Double_t xmin = -2.0;
  Double_t xmax = +2.0;
  Double_t ymin = -2.0;
  Double_t ymax = +2.0;
 
  TH2D h2("h2", "", 1, xmin, xmax, 1, ymin, ymax);
 
  h2.SetStats(kFALSE);
  h2.Draw();
 
 
  ellipse.SetLineWidth(2);
  ellipse.Draw();
 
 
  TGraph g(N, x, y);
 
  g.SetMarkerStyle(20);
  g.SetMarkerColor(kBlack);
  g.SetMarkerSize(0.7);
  g.Draw("P");
 
 
  // lower hull
  
  N = distance(buffer.begin(), hull.first);
 
  TGraph g1(N, x, y);
 
  g1.SetLineColor(kBlack);
  g1.SetLineWidth(2);
  g1.Draw("L");
 
 
  // upper hull
 
  const int i = --N;
 
  N = distance(buffer.begin(), hull.second);
 
  x[N] = x[0];
  y[N] = y[0];
 
  N = N + 1 - i;
 
  TGraph g2(N, &x[i], &y[i]);
 
  g2.SetLineColor(kRed);
  g2.SetLineWidth(2);
  g2.Draw("L");
 
 
  cv.Update();
 
  tp->Run();
}