JGetSmallestDistance2D.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 "TCanvas.h"
#include "TApplication.h"

#include "JGeometry2D/JVector2D.hh"
#include "JGeometry2D/JGeometry2DToolkit.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JTimer.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"


/**
 * \file
 *
 * Example program to find smallest distance between two points.
 * \author mdejong
 */
int main(int argc, char* argv[])
{
  using namespace std;

  string  inputFile;
  string  outputFile;
  int     numberOfPoints;
  UInt_t  seed;
  int     debug;

  try {

    JParser<> zap("Example program to find smallest distance between two points.");

    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;


  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 > 1) {

    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 << setw(7) << i->getX() << ' ' 
            << setw(7) << i->getY() << endl;
      
      out.close();
    }
  }

  if (buffer.size() < 2) {
    FATAL("Not enough points." << endl);
  }

  for (const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
    DEBUG(i->getX() << ' ' << i->getY() << endl);
  }


  {
    JTimer timer("classic");
  
    timer.start();

    double dmin = numeric_limits<double>::max();
    
    for (const_iterator i = buffer.begin(); i != buffer.end(); ++i) {
      for (const_iterator j = i; ++j != buffer.end(); ) {
        
        const double d = i->getDistance(*j);
        
        if (d < dmin) {
          dmin = d;
        }
      }
    }
    
    timer.stop();

    cout << "Minimal distance " << SCIENTIFIC(12,5) << dmin << endl;
    timer.print(cout);
  }


  JTimer timer("O(n log(n))");
  
  timer.start();
  
  const double dmin = getSmallestDistance2D(buffer.begin(), buffer.end());

  timer.stop();

  cout << "Minimal distance " << SCIENTIFIC(12,5) << dmin << endl;
  timer.print(cout);
  

  pair<const_iterator, const_iterator> result = JSmallestDistance2D::getPair(buffer.begin(), buffer.end(), dmin);

  cout << "pair " 
       << "(" << FIXED(7,5) << result.first ->getX() << "," << FIXED(7,5) << result.first ->getY() << ")" << ' ' 
       << "(" << FIXED(7,5) << result.second->getX() << "," << FIXED(7,5) << result.second->getY() << ")" << endl; 
  cout << "Distance         " << SCIENTIFIC(12,5) << result.first->getDistance(*result.second) << 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);


  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 = -1.1;
  Double_t xmax = +1.1;
  Double_t ymin = -1.1;
  Double_t ymax = +1.1;

  TH2D h2("h2", "", 1000, xmin, xmax, 1000, ymin, ymax);

  h2.SetStats(kFALSE);
  h2.Draw();


  TGraph g(N, x, y);

  g.SetMarkerStyle(20);
  g.SetMarkerColor(kBlack);
  g.SetMarkerSize(0.7);
  g.Draw("P");
  
  TMarker m1(result.first ->getX(), result.first ->getY(), 20);
  TMarker m2(result.second->getX(), result.second->getY(), 20);
  
  m1.SetMarkerColor(kRed);
  m2.SetMarkerColor(kRed);
  
  m1.SetMarkerSize(0.7);
  m2.SetMarkerSize(0.7);

  m1.Draw();
  m2.Draw();

  cv.Update();

  tp->Run();
}