// Test Delaunay Triangularization in ROOT

#include <iostream>

#include "TCanvas.h"
#include "TRandom3.h"
#include "TFile.h"
#include "TGraph2D.h"
#include "TGraphDelaunay.h"

void test_delaunay(){

  //
  // PART 1: create Delaunay triangulation and write object to file
  //

  TCanvas* c = new TCanvas("c","Delaunay test", 0, 0, 700, 600);
  TRandom3* myrnd = new TRandom3( 0 );

  const int n = 100; // number of points that will be interpolated

  double x[n];
  double y[n];
  double z[n];

  TFile* fout = new TFile( "test_delaunay.root", "RECREATE");

  TGraph2D* gr = new TGraph2D();
  gr->SetName("MyGraph2D");
  // create random set of 3D points
  for( unsigned int i=0; i<n; i++){
    x[i] = myrnd->Uniform();
    y[i] = myrnd->Uniform();
    z[i] = myrnd->Uniform();
    gr->SetPoint( i, x[i], y[i], z[i] );
  }
  
  gr->SetMarkerStyle( 24 );
  gr->SetMarkerSize( 1 );
  gr->Draw("TRI1");

  // create Delaunay interpolant from TGraph
  TGraphDelaunay* grdt = new TGraphDelaunay( gr );
  grdt->FindAllTriangles();
  double xq = 0.5;
  double yq = 0.5;
  double zq = grdt->ComputeZ( xq, yq );
  printf("Z value for (%f, %f) = %f\n",xq, yq, zq);
  
  // write the TGraphDelaunay to a ROOT file
  grdt->SetName("MyDelaunay");
  grdt->Write();
  gr->Write();
  fout->Close();


  //
  // PART 2: re-read the file
  //

  TFile *fin  = new TFile("test_delaunay.root", "READ");

  TGraph2D* gr2 = (TGraph2D*)fin->Get( "MyGraph2D" );
  if ( gr2 ){
    TCanvas* c2 = new TCanvas("c2","Delaunay test - read", 100, 100, 700, 600);
    gr2->SetMarkerStyle( 24 );
    gr2->SetMarkerSize( 1 );
    gr2->Draw("TRI1");
  }

  TGraphDelaunay* d = (TGraphDelaunay*)fin->Get( "MyDelaunay" );
  zq = 0.0;
  if ( d ){
    printf("Found object %s\n",d->GetName());
    // d->FindAllTriangles();
    xq = 0.5;
    yq = 0.5;
    zq = d->ComputeZ( xq, yq );
    printf("Z value for (%f, %f) = %f\n",xq, yq, zq);
  }

  gr2->SetDirectory(0);
  fin->Close();



}