#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <math.h>

#include <typeinfo>

#include <vector>
#include <string>

#include "ctimer.h"
#include "piacmd.h"
#include "nobjmgr.h"
#include "pistdimgapp.h"

#include "pidrawer.h"
#include "nomgadapter.h"

#include "histos.h"
#include "histos2.h"
#include "hisprof.h"
#include "ntuple.h"

#include "servnobjm.h"
#include <nbtri.h> 

#include "picntools.h"
#include "pigncont.h"
//A-Virer #include "pidrwtools.h"
//A-Virer #include "piaxestools.h"

extern "C" { 
#include "contour.h"
#include "alloc.h"
#include "set.h"
}

TBOOLEAN multiplot;
double zero =0;
TBOOLEAN polar;
int xleft, xright, ybot, ytop;
float xsize; 
float ysize; 
float xoffset, yoffset;


GNUPlotContour::GNUPlotContour(struct iso_curve * iso){
  int k=0;
  iso1 = iso;
  struct iso_curve *iso_cur = iso ;
  struct iso_curve *iso_old;
  _xmin = 1.e37;
  _ymin = 1.e37;
  _xmax = -1.e37;
  _ymax = -1.e37;
  _myiso = NULL;
  _zmin = 1.e37;
  _zmax = -1.e37;
  while(iso_cur){
    k++;
    iso_cur = iso_cur->next;
    _ny = iso_cur->p_count;
    for(int k=0 ; k<_ny ; k++){
      double xx=iso->points[k].x;
      double yy=iso->points[k].y;
      double zz=iso->points[k].z;
      if(xx<_xmin) _xmin = xx;
      if(xx>_xmax) _xmax = xx;
      if(yy<_ymin) _ymin = yy;
      if(yy>_ymax) _ymax = yy;
      if(zz<_zmin) _zmin = zz;
      if(zz>_zmax) _zmax = zz;

    }
  }

  _nx = k;
  _contours = NULL;
  My_Levels = NULL;
  
  _npolys=-1;  
  _Dxp =_Dyp =1.;
  _Exy =_Invy =_Invy = FALSE ;

}


GNUPlotContour::GNUPlotContour(NTupleInterface *nt,int *nz=NULL,double *z=NULL,double *dz=NULL){
  _xmin = 1.e37;
  _ymin = 1.e37;
  _xmax = -1.e37;
  _ymax = -1.e37;
  _zmin = 1.e37;
  _zmax = -1.e37;
  char *nom[4]={"x","y","z","niso"};
  _myiso = new NTuple(4,nom);  
  //_myiso->Show();
  float xnt[4];
  double mn,mx;

  nt->GetMinMax(0,_xmin,_xmax);
  nt->GetMinMax(1,_ymin,_ymax);
  nt->GetMinMax(2,_zmin,_zmax);

  
  struct iso_curve *istmp = NULL ;
  struct iso_curve *iso = NULL ;
  struct iso_curve *iso_old = NULL ;
  double x,y;
  int nen = nt->NbLines();
  
  // histo des y
  for(int k= 0 ; k<1000 ; k++){
    istmp = iso_alloc(100);
    iso_free(istmp);
  }

#define NUMPTMIN 25
#define MAXISO 100
#define MINISO 5

  Histo prep(_ymin,_ymax,nen);  
  for(int ix=0 ; ix<nen ; ix++) prep.Add(nt->GetCell(ix,1),1.);
  Histo Intg = prep;
  Intg.HInteg(0);
  //prep.Print();
  //Intg.Print();
  float *bin;
  int *cnt;

  int niso = sqrt( (double)nen);
  bin = new float[niso+1];
  cnt = new int[niso+1];

  for(int i=0 ; i<niso ; i++){    
    bin[i] = prep.BinHighEdge(prep.BinPercent( ((double)i)/(double)niso));
    cnt[i] = Intg(prep.BinPercent( ((double)i)/(double)niso)); 
  }

  bin[niso]=_ymax;
  
  cnt[niso] = prep.NEntries();
  int numx = 0;
  int numi = 9999;
  for(int i=0;i<niso ; i++){
    int numy = cnt[i+1]-cnt[i];
    
    if(numy>0){
      numx = numx>numy ? numx : numy ;
      numi = numy>numi ? numi : numy ;
    }
  }

  
  if (numi<=0) {
    
    return;
  }
  _nx = 0;
  
  for(int i=0;i<niso ; i++){
    double ymi = bin[i];
    double ymx = bin[i+1];
    int numy = cnt[i+1]-cnt[i];
    
    if(numy<=0) continue;
    double *xp;
    xp = new double[numy];
    
    int *indx;
    indx = new int[numy];

    istmp = iso_alloc(numy) ; // structure temporaire 
    
    int ip=0;
    
    for(int ix=0 ; ix<nen ; ix++){
      
      if(nt->GetCell(ix,1)>=ymi&&nt->GetCell(ix,1)<ymx){
        

	istmp->points[ip].type = INRANGE;
	istmp->points[ip].x = nt->GetCell(ix,0);
	istmp->points[ip].y = nt->GetCell(ix,1);
	istmp->points[ip].z = nt->GetCell(ix,2);
	indx[ip] = ip;
	xp[ip++] = nt->GetCell(ix,0);
      }      
    }
    
    // tri des points de l'iso_curve 
    //
    istmp->p_count = ip;
    
    tri_double(xp,indx,ip);
    
    // remplissage structure finale 
    
    iso = iso_alloc(numx) ;
    if( _nx == 0 ) iso1 = iso; //premiere structure remplie 
    if (iso_old) iso_old->next = iso;
    for(int ix=0 ; ix<numy   ; ix++){
      int jx = indx[ix];
    
      
      iso->points[ix].type = INRANGE ;
      iso->points[ix].x = istmp->points[jx].x ;
      iso->points[ix].y = istmp->points[jx].y ;
      iso->points[ix].z = istmp->points[jx].z ;
      xnt[0] = iso->points[ix].x;
      xnt[1] = iso->points[ix].y;
      xnt[2] = iso->points[ix].z;
      xnt[3] = _nx;
      //cout << " fill avec "<<xnt[0]<<"  "<<xnt[1] <<" "<<xnt[2] << " "<<xnt[3]<<endl;
      // _myiso->Show();
 
      _myiso->Fill(xnt);
    }
    //cout << " completion "<<ip<<" =>> "<<numx<<endl;
    // on complete 
    for( int ix=ip; ix<numx ; ix++ ){
      int jmx = indx[ip-1];
      iso->points[ix].type = INRANGE ;
      iso->points[ix].x = istmp->points[jmx].x ;
      iso->points[ix].y = istmp->points[jmx].y ;
      iso->points[ix].z = istmp->points[jmx].z ;
      xnt[0] = iso->points[ix].x;
      xnt[1] = iso->points[ix].y;
      xnt[2] = iso->points[ix].z;
      xnt[3] = _nx;
      //cout << " fill avec "<<xnt[0]<<"  "<<xnt[1] <<" "<<xnt[2] << " "<<xnt[3]<<endl;
      _myiso->Fill(xnt);
      //_myiso->Show();
    }
    
    iso->p_count = numx;

    iso_old = iso;
    _nx++;
    //cout << " destruction structure tempo "<<_nx <<endl;
    if(istmp){
      iso_free( istmp );// destruction de la structure temporaire 
      istmp = NULL;
    }
    if(xp){delete[] xp ; xp=NULL;}
    if(indx) {delete[] indx ; indx = NULL;}

  }
  //if(_myiso!=NULL)_myiso->Write("myiso.ppf"); 
  //cout << " fin du remplissage des iso_curves xmin,max "<< _xmin<<","<<_xmax<<" ymin,max "<<_ymin<<","<<_ymax<<endl;
  _contours = NULL;
  My_Levels = NULL;
  _npolys=-1;  
  Contour_Levels_kind = LEVELS_AUTO ;
  Contour_kind = CONTOUR_KIND_LINEAR ;
  set_contour_kind(Contour_kind);
  Contour_Levels = 5;
  set_contour_levels(Contour_Levels);
  set_contour_levels_kind(Contour_Levels_kind);
  
 
  if(dz!=NULL){
    My_Levels = new double[2];
    Contour_Levels = *nz;
    My_Levels[0] = *z ;
    My_Levels[1] = *dz;
    Contour_Levels_kind = LEVELS_INCREMENTAL ;
    set_contour_levels(Contour_Levels);
    set_contour_levels_kind(Contour_Levels_kind);
    set_contour_levels_list(My_Levels);//,2);

  }else  if(nz!=NULL){
    
    Contour_Levels = *nz;
    Contour_Levels_kind = LEVELS_NUM ;
    set_contour_levels(Contour_Levels);
    set_contour_levels_kind(Contour_Levels_kind);
  }
  
}

GNUPlotContour::GNUPlotContour(P2DArrayAdapter*arr,int *nz=NULL,double *z=NULL,double *dz=NULL){
  _xmin = 1.e37;
  _ymin = 1.e37;
  _xmax = -1.e37;
  _ymax = -1.e37;
  _zmin = 1.e37;
  _zmax = -1.e37;
  // _nx = arr->XSize()-1;
  // _ny = arr->YSize()-1;
   _nx = arr->XSize()  ;
   _ny = arr->YSize()  ;
   //cout << " size x,y "<<_nx<<" , "<<_ny <<endl;
   char *nom[4]={"x","y","z","niso"};
   
   _myiso = new NTuple(4,nom);  
   float xnt[4];
   struct iso_curve *iso =NULL ;
   struct iso_curve *iso_old =NULL ;
   double x,y;
   // Calcul de la taille en X et en Y d'un pixel
   double dxp, dyp;
   double x1,y1;
   arr->XYfromxy(0,0,x,y);
   arr->XYfromxy(1,1,x1,y1);
   _Dxp = (x1-x);
   _Dyp = (y1-y);
   double zz;
   for(int ix=0 ; ix<_nx ; ix++){
     iso = iso_alloc(_ny) ;
     if(ix==0) iso1 = iso;
     if (iso_old) iso_old->next = iso;
     for(int iy = 0 ; iy<_ny ; iy++){
       /*i = ix*100+iy;*/
       arr->XYfromxy(ix,iy,x,y);

       iso->points[iy].type = INRANGE;
       iso->points[iy].x = x + _Dxp/2.;
       iso->points[iy].y = y + _Dyp/2.;
       zz=arr->Value(ix,iy);
       iso->points[iy].z = zz;
       //iso->points[iy].z = ((x-25)*(x-25) + (y-25)*(y-25))*.05 ;
       if(x <_xmin) _xmin = x;
       if(x >_xmax) _xmax = x;
       if(y <_ymin) _ymin = y;
       if(y >_ymax) _ymax = y;
       if(zz <_zmin) _zmin = zz;
       if(zz >_zmax) _zmax = zz;
       //printf("ix=%d iy=%d x=%f y=%f z=%f \n",ix,iy,iso->points[iy].x,iso->points[iy].y,iso->points[iy].z);
       xnt[0] = iso->points[iy].x;
       xnt[1] = iso->points[iy].y;
       xnt[2] = iso->points[iy].z;
       xnt[3] = _nx;

       _myiso->Fill(xnt);
     }
     iso->p_count = _ny;
     iso_old = iso;
     /*printf(" remplissage %lx old %lx old next %lx \n",iso,iso_old,iso->next,iso_old->next);*/
     /*iso->next = iso ;*/
   }
   //cout << " fin du remplissage des iso_curves xmin,max "<< _xmin<<","<<_xmax<<" ymin,max "<<_ymin<<","<<_ymax<<endl;
   _contours = NULL;
   My_Levels = NULL;
  _npolys=-1;  
  Contour_kind = CONTOUR_KIND_LINEAR ;
  set_contour_kind(Contour_kind);
  Contour_Levels_kind = LEVELS_AUTO ;
  Contour_Levels = 5;
  set_contour_levels(Contour_Levels);
  set_contour_levels_kind(Contour_Levels_kind);
  

  if(dz!=NULL){
    My_Levels = new double[2];
    Contour_Levels = *nz;
    My_Levels[0] = *z ;
    My_Levels[1] = *dz;
    Contour_Levels_kind = LEVELS_INCREMENTAL ;
    set_contour_levels(Contour_Levels);
    set_contour_levels_kind(Contour_Levels_kind);
    set_contour_levels_list(My_Levels);//,2);

  }else  if(nz!=NULL){
    Contour_Levels = *nz;
    Contour_Levels_kind = LEVELS_NUM ;
    set_contour_levels(Contour_Levels);
    set_contour_levels_kind(Contour_Levels_kind);
  }
  
}

GNUPlotContour::~GNUPlotContour(){
 cout << " destructeur de GNUPlotContour "<<endl;
  
 struct iso_curve *iso_cur;
 struct iso_curve *iso_nx;
 iso_cur = iso1;
 
  while(iso_cur){
    // cout << " ~GNUPlotContour() : destruction de  iso1 "<< iso_cur << endl;
   iso_nx = iso_cur->next;  
   iso_free(iso_cur);iso_cur = NULL;
   
   iso_cur = iso_nx;
   
   }
 
 iso1 = NULL;
 
 struct gnuplot_contours *cntcur = _contours;
 struct gnuplot_contours *cntold;

 while(cntcur) {
   //cout << " ~GNUPlotContour() : destruction de _contours" << _contours <<endl; 
   cntold = cntcur;
   cntcur = cntold->next;
   gp_free(cntold);
   cntold=NULL;
 } 
 
 _contours = NULL;

 if(My_Levels) {
   //cout << " ~GNUPlotContour() : destruction de MyLevels "<<My_Levels <<endl;
   delete[] My_Levels; My_Levels=NULL;}

 if(_myiso!=NULL){ delete _myiso; _myiso=NULL;}

} 

void GNUPlotContour::CalcContour(){
 cout << " GNUPlotContour::CalcContour(): determination des contours "<<endl;

 set_contour_kind(Contour_kind);
 if(Contour_Levels_kind == LEVELS_INCREMENTAL){
    
    set_contour_levels(Contour_Levels);
    set_contour_levels_kind(Contour_Levels_kind);
    //free_contour_levels_list();
    set_contour_levels_list(My_Levels);//,2);

  }else  if(Contour_Levels_kind == LEVELS_DISCRETE ){
    set_contour_levels(Contour_Levels);
    set_contour_levels_kind(Contour_Levels_kind);
    //free_contour_levels_list();
    set_contour_levels_list(My_Levels);//,Contour_Levels);
    
  }else  if(Contour_Levels_kind == LEVELS_NUM ){
    
    set_contour_levels(Contour_Levels);
    set_contour_levels_kind(Contour_Levels_kind);
  }else if(Contour_Levels_kind == LEVELS_AUTO){
    Contour_Levels = 5;
    set_contour_levels(Contour_Levels);
    set_contour_levels_kind(Contour_Levels_kind);
  }

 if(_contours) {
   //cout << " GNUPlotContour::CalcContour(): destruction des contours "<<endl;
   struct gnuplot_contours *cntcur = _contours;
   struct gnuplot_contours *cntold;
   while(cntcur) {
     cntold = cntcur;
     cntcur = cntold->next;
     gp_free(cntold);
   }
   _contours = NULL;
   
 }

   
 
 //struct gnuplot_contours *cntcur;
 //struct gnuplot_contours *cntold;
 
 _contours = contour (_nx,iso1);
 //free_contour_levels_list();

 
}
//_+_ Methode 
void 
GNUPlotContour::SetMyLevels(double *ptr,int k){
  
  if(My_Levels!=NULL) delete My_Levels;
  My_Levels = new double[k];
  for(int i=0 ; i<k ; i++)My_Levels[i] = ptr[i];
}

// Class PIContourDrawer
/* --Methode-- */
PIContourDrawer::PIContourDrawer(P2DArrayAdapter* arr,bool autodel,int *nz,double *z0,double *dz )
  :GNUPlotContour(arr,nz,z0,dz)

{
  _arr = arr;
  _nti = NULL;
  _nz=5;
  if(nz!=NULL)_nz = *nz;
  _autodel = autodel;
  InitAtts();
}

PIContourDrawer::PIContourDrawer(NTupleInterface* nti,bool autodel,int *nz,double *z0,double *dz )
  :GNUPlotContour(nti,nz,z0,dz)

{
  _arr = NULL;
  _nti = nti;
  _nz=5;
  if(nz!=NULL)_nz = *nz;
  _autodel = autodel;
  InitAtts();
}



/* --Methode-- */
PIContourDrawer::~PIContourDrawer()
{
  if(_autodel){ 
    if(_arr!=NULL){delete _arr; _arr=NULL;}
    if(_nti!=NULL){delete _nti ; _nti=NULL;}
  }
  //if (mColorMap!=NULL){delete  mColorMap;  mColorMap=NULL;}
}

/* --Methode-- */
void PIContourDrawer::Draw(PIGraphicUC* g, double xmin, double ymin, double xmax, double ymax)
{
  // double xmin, double ymin, double xmax, double ymax LIMITES DE LA ZONE A REFRESHER
  //PIGrCoord * xx = new PIGrCoord[10];
  //PIGrCoord * yy = new PIGrCoord[10];
  PIGrCoord * xx =NULL;
  PIGrCoord * yy =NULL;
  PIGrCoord xa;
  PIGrCoord ya;
  char buff[64];
  // On met les bons attributs graphiques
  //SelGraAtt(g);
  g->SaveGraphicAtt();
  struct gnuplot_contours *cntcur;
  struct gnuplot_contours *cntold;
  cntcur = _contours;

  double a,b;
  int tot=0;
  bool Invx, Invy, Exy;
  PIColorMap * cmap = NULL;
  if(mCmapid !=CMAP_OTHER ){
    cmap = new PIColorMap(mCmapid);
  }
  int numdr = 0;
  while(cntcur){
    int npts = cntcur->num_pts ;
    if(npts<0) continue;

    xx = new PIGrCoord[npts];
    yy = new PIGrCoord[npts];
    double lev = (cntcur->coords)[0].z ;
    for(int l=0 ; l< npts ; l++){
      a = (cntcur->coords)[l].x;
      b = (cntcur->coords)[l].y;
      
      xx[l] = a ;
      yy[l] = b ;

    }

    // choix de la couleur 
    
    g->SelForeground(mFCol);
    if(mCmapid !=CMAP_OTHER ){
      if(_zmax-_zmin!=0){
	int kc = ((lev - _zmin)/(_zmax-_zmin))*cmap->NCol();	
	
	g->SelForeground(*cmap,kc);
      }
    }    
    // traitement des  des markers 
    if(IsMarkOn()==true){
      
      g->SelMarker(mMSz,mMrk);
      g->DrawMarkers(xx,yy,npts);
    }
    // traitement des lignes 
    if(mLineOn==true){
      
      g->SelLine(mLAtt);
      g->DrawPolygon(xx,yy,npts,false);
    }
    // traitement des labels     
    // SIMPLISTE POUR L'INSTANT : UN AFFICHAGE 
    if(mLabelOn==true){
      // 
      
      char strg[10];
      sprintf(strg,"%g",lev);
      PIFont myfont(mFName);
      myfont.SetFontAtt(mFAtt);
      myfont.SetFontSz(mFSz);
      
      g->SelFont(myfont);
      double px,py;
      px = (cntcur->coords)[0].x+2*(Xmax()-Xmin())/100.;
      py = (cntcur->coords)[0].y;
      g->DrawString(px,py,strg, PI_HorizontalCenter&&PI_VerticalCenter );
      
    }
    numdr++;
    delete[] xx;
    xx=NULL;
    delete[] yy;
    yy=NULL;
    cntcur = cntcur->next;
    tot++;
  }    
  
  //cout << "PIContourDrawer::Draw  fin de trace des "<<tot<< " polygones "<<endl; 
 
  // 
  if(cmap !=NULL)delete cmap;
  
  g->RestoreGraphicAtt();
}

/* --Methode-- */
void PIContourDrawer::UpdateLimits()
{
  // Doit calculer les limites

  double xmn = _xmin;
  double xmx = _xmax;  
  double ymn = _ymin;
  double ymx = _ymax;
  
  SetLimits(xmn, xmx, ymn, ymx);
}


/* --Methode-- */
void PIContourDrawer::ShowControlWindow(PIBaseWdgGen* wdg)
{
  PICnTools::SetCurrentBaseWdg(wdg);
  PICnTools::SetCurrentCnDrw(this);
  PICnTools::ShowPICnTools();
}


/* --Methode-- */
bool PIContourDrawer::IsLabelOn(){
  return (mLabelOn);
}

bool PIContourDrawer::IsLineOn(){
  return (mLineOn);
}

bool PIContourDrawer::IsMarkOn(){
  return (mMarkOn);
}

void PIContourDrawer::SetLabelOn(bool state){
  mLabelOn = state;
}

void PIContourDrawer::SetMarkOn(bool state){
  mMarkOn = state;
}

void PIContourDrawer::SetLineOn(bool state){
  mLineOn = state;
}

void PIContourDrawer::InitAtts(){
  mLineOn = false;
  mMarkOn = true;
  mLabelOn = false;
}