#include <math.h>
#include "pisurfdr.h"

/* --Methode-- */
PISurfaceDrawer::PISurfaceDrawer(P2DArrayAdapter* arr, bool asxy, bool asz, bool ad)
{
  mArr = arr;  mAda = ad;
  mScx = mScy = mScz = 1.;
  mZmin = 0.; mZmax = 1.;
  mLut = NULL;
  if (arr == NULL)  return;
  arr->XYfromxy(0,0,mXmin, mYmin);
  arr->XYfromxy(arr->XSize(),arr->YSize(),mXmax, mYmax);
  if (asxy)  mScy = (mXmax-mXmin)/(mYmax-mYmin);
  double scm =  ((mXmax-mXmin) > (mYmax-mYmin)) ? (mXmax-mXmin)  : (mYmax-mYmin);
  int i,j;   double cv;
  mZmin = mZmax = arr->Value(0,0);
  for(i=0; i<arr->XSize(); i++) 
    for(j=0; j<arr->YSize(); j++) {
      cv  = arr->Value(i,j);
      if (cv < mZmin) mZmin = cv;
      if (cv > mZmax) mZmax = cv;
    }
  if ((mZmax-mZmin) <= 1.e-19)  mZmax = mZmin + 1.e-19;
  if (asz) mScz = scm/(mZmax-mZmin);
//  printf("PISurfaceDrawer::PISurfaceDrawer Min/Max: %g,%g  %g,%g  %g,%g  Sc: %g %g %g \n",
//          mXmin, mXmax, mYmin, mYmax, mZmin, mZmax, mScx, mScy, mScz);
  double zmin = mZmin; 
  double zmax = mZmax;
  double dzz = (mZmax-mZmin)*0.05;
  mLut = new LUT(zmin-dzz*2, zmax-dzz, 256);
}

/* --Methode-- */
PISurfaceDrawer::~PISurfaceDrawer()
{
  if (mAda) delete mArr;
  if (mLut) delete mLut;
}

/* --Methode-- */
void PISurfaceDrawer::UpdateLimits()
{
  if (mArr == NULL)  return;
  double xmin, xmax, ymin, ymax, zmin, zmax;
  xmin = mXmin*mScx;  xmax = mXmax*mScx;
  ymin = mYmin*mScy;  ymax = mYmax*mScy;
  zmin = mZmin*mScz;  zmax = mZmax*mScz;
// Centre du champ en C = (xmin+xmax)/2., (ymin+ymax)/2  (zmin+zmax)*0.5
// Distance D = Max(xmax-xmin,ymin-ymax)*2
// Observateur en O = X+D, Yc+2*D 
  double D = xmax-xmin;
  if (D < (ymax-ymin))  D = ymax-ymin;
  D *= 1.4;
  
  Set3DView((xmin+xmax)/2., (ymin+ymax)/2, zmin+(zmax-zmin)*0.40, 
            (xmin+xmax)/2.+D , (ymin+ymax)/2.-2.5*D , zmin+(zmax-zmin)*0.85, 0.25, 0.25);  
  
  x3Min = xmin;   x3Max = xmax;
  y3Min = ymin;   y3Max = ymax;
  z3Min = zmin-0.05*(zmax-zmin);   z3Max = zmax+0.05*(zmax-zmin);

//  printf("PISurfaceDrawer::UpdateLimits() : %g-%g  %g-%g  %g-%g (%g) \n", xmin,xmax,ymin,ymax,zmin,zmax,D);
//  printf("PISurfaceDrawer::UpdateLimits() :  %g %g %g <<- %g %g %g \n", 
//         (xmin+xmax)/2., (ymin+ymax)/2, zmin+(zmax-zmin)*0.40, 
//         (xmin+xmax)/2.+D , (ymin+ymax)/2.+2.*D , zmin+(zmax-zmin)*0.60 );

}


typedef struct   /*  structures pour trier en ZP les facettes */
  { double dp;  double zh;  int ij;  } ZPSORT;
/* Nouvelle-Fonction */
static int ijzpsort(const void * sp1, const void * sp2)
/*  Fonctions de tri des etoiles en flux ds un starlist */
{
  ZPSORT* p1 = (ZPSORT*)sp1;   
  ZPSORT* p2 = (ZPSORT*)sp2;  
  if (p1->dp > p2->dp)  return(-1);
  else if (p1->dp < p2->dp)  return(1);
  return(0);
}

/* --Methode-- */
void PISurfaceDrawer::DrawAxes(PIGraphicUC* g)
{
  if (axesFlags == kAxesNone) return;
#if defined(__DECCXX) || defined(__KCC__)
  PIGraphic3D* g3 = dynamic_cast<PIGraphic3D*>(g);
#else
  PIGraphic3D* g3 = (PIGraphic3D*)(g);
#endif
  Draw3DBoxe(g3, XMin3(), XMax3(), YMin3(), YMax3(), ZMin3(), ZMax3(),
                 mXmin, mXmax, mYmin, mYmax, mZmin-0.05*(mZmax-mZmin), mZmax+0.05*(mZmax-mZmin));
}

/* --Methode-- */
void PISurfaceDrawer::Draw(PIGraphicUC* g, double xmin, double ymin, double xmax, double ymax)
{
// On trace les axes - En attendant de faire mieux - Reza 8/12/98
  if (axesFlags != kAxesNone)  DrawAxes(g);

  if (mArr == NULL)  return;

#if defined(__DECCXX) || defined(__KCC__)
  PIGraphic3D* g3 = dynamic_cast<PIGraphic3D*>(g);
#else
  PIGraphic3D* g3 = (PIGraphic3D*)(g);
#endif

  g3->SelBackground(PI_White);
//  g3->Erase(xmin, ymin, xmax-xmin, ymax-ymin);

  int xsz = mArr->XSize();
  int ysz = mArr->YSize();

  ZPSORT* ijzp = new ZPSORT[xsz*ysz];
  double* xw = new double[(xsz+1)*(ysz+1)];
  double* yw = new double[(xsz+1)*(ysz+1)];
  double* zw = new double[(xsz+1)*(ysz+1)];

  double dx,dy;   // taille du bin 
  dx = (mXmax-mXmin)*mScx/xsz; 
  dy = (mYmax-mYmin)*mScy/ysz;
 
  int i1,j1;
  int i,j,k,l;
  double x,y,z,xp,yp,zp;

  k = 0;
  for(i=0; i<xsz; i++) 
    for(j=0; j<ysz; j++) { 
      ijzp[k].ij = j*xsz+i; 
      mArr->XYfromxy(i,j,x,y);
      x *= mScx; x += dx/2.;
      y *= mScy; y += dy/2.;
      ijzp[k].zh = mArr->Value(i,j);
      z = ijzp[k].zh * mScz;
      g3->C3DC2ObsCS(x,y,z,xp,yp,zp);
      ijzp[k].dp = zp; 
      k++;

// Calcul de la projection du coin de la facette
      x -= dx/2.;   y -= dy/2.;
      i1 = (i>0) ? i-1 : 0;
      j1 = (j>0) ? j-1 : 0;
      z = mArr->Value(i,j) + mArr->Value(i1,j) + mArr->Value(i1,j1) +  mArr->Value(i,j1);
      z = z/4.*mScz;        
//      g3->Proj3DC2GrC(x,y,z, xw[j*(xsz+1)+i], yw[j*(xsz+1)+i]);
      xw[j*(xsz+1)+i] = x;      yw[j*(xsz+1)+i] = y;    zw[j*(xsz+1)+i] = z;    
    }

  j = ysz;
  for(i=0; i<xsz; i++) {
     mArr->XYfromxy(i,j-1,x,y);
     x *= mScx; 
     y *= mScy; y += dy;
     i1 = (i>0) ? i-1 : 0;
     j1 = j-1;
     z = mArr->Value(i,j-1) + mArr->Value(i1,j-1) + mArr->Value(i1,j1) +  mArr->Value(i,j1);
     z = z/4.*mScz;
//     g3->Proj3DC2GrC(x,y,z, xw[j*(xsz+1)+i], yw[j*(xsz+1)+i]);     
     xw[j*(xsz+1)+i] = x;      yw[j*(xsz+1)+i] = y;   zw[j*(xsz+1)+i] = z;   
  }
  i = xsz;
  for(j=0; j<ysz; j++) {
     mArr->XYfromxy(i-1,j,x,y);
     x *= mScx; x += dx;
     y *= mScy; 
     i1 = i-1;
     j1 = (j>0) ? j-1 : 0;
     z = mArr->Value(i-1,j) + mArr->Value(i1,j) + mArr->Value(i1,j1) +  mArr->Value(i-1,j1);
     z = z/4.*mScz;
//     g3->Proj3DC2GrC(x,y,z, xw[j*(xsz+1)+i], yw[j*(xsz+1)+i]);     
     xw[j*(xsz+1)+i] = x;      yw[j*(xsz+1)+i] = y;   zw[j*(xsz+1)+i] = z;   
  }
  i = xsz;  j = ysz;
  mArr->XYfromxy(i-1,j-1,x,y);
  i1 = i-1;   j1 = j-1;
  x *= mScx; x += dx;
  y *= mScy; y += dy;
  z = mArr->Value(i-1,j-1) + mArr->Value(i1,j-1) + mArr->Value(i1,j1) +  mArr->Value(i-1,j1);
  z = z/4.*mScz;
//  g3->Proj3DC2GrC(x,y,z, xw[j*(xsz+1)+i], yw[j*(xsz+1)+i]);     
  xw[j*(xsz+1)+i] = x;      yw[j*(xsz+1)+i] = y;    zw[j*(xsz+1)+i] = z;   

// On trie par ordre de distance zp de projection / observateur, du plus loin au plus proche
  qsort(ijzp,(size_t) k,sizeof(ZPSORT),ijzpsort);

  PIGrCoord xg[5],yg[5],zg[5];
  PIColors fgc; 


  if (mCmapid == CMAP_OTHER) {   // Trace avec des lignes
    if (mLAtt == PI_NotDefLineAtt)    g3->SelLine(PI_ThinLine);
    fgc = ( mFCol == PI_NotDefColor ) ? PI_Black : mFCol;
    for(l=0;l<k;l++) {
      i =  ijzp[l].ij%xsz;
      j =  ijzp[l].ij/xsz;
      xg[4] = xg[0] = xw[i+j*(xsz+1)];    yg[4] = yg[0] = yw[i+j*(xsz+1)];  
      zg[4] = zg[0] =   zw[i+j*(xsz+1)];
      xg[1] = xw[i+1+j*(xsz+1)];          yg[1] = yw[i+1+j*(xsz+1)];   
      zg[1] = zw[i+1+j*(xsz+1)];   
      xg[2] = xw[i+1+(j+1)*(xsz+1)];      yg[2] = yw[i+1+(j+1)*(xsz+1)];
      zg[2] = zw[i+1+(j+1)*(xsz+1)];
      xg[3] = xw[i+(j+1)*(xsz+1)];        yg[3] = yw[i+(j+1)*(xsz+1)];  
      zg[3] = zw[i+(j+1)*(xsz+1)]; 
// On efface derriere 
      g3->SelForeground(PI_White);
//      g3->BaseGraphic()->DrawFPolygon(xg,yg,5);
      g3->DrawFPolygon3D(xg,yg,zg,5);
// On trace 
      g3->SelForeground(fgc);
//      g3->BaseGraphic()->DrawPolygon(xg,yg,5);
      g3->DrawPolygon3D(xg,yg,zg,5);
    }
  }
  else {   // Trace avec niveaux de couleur
    PIColorMap* cmap = new PIColorMap(mCmapid);
    PIColors fgc; 
    if ( (mCmapid == CMAP_GREY32) || (mCmapid == CMAP_GREYINV32) || 
         (mCmapid == CMAP_GREY128)  || (mCmapid == CMAP_GREYINV128) )
      fgc = ( mFCol == PI_NotDefColor ) ? PI_Red : mFCol;
    else  fgc = ( mFCol == PI_NotDefColor ) ? PI_Black : mFCol;
    for(l=0;l<k;l++) {
      i =  ijzp[l].ij%xsz;
      j =  ijzp[l].ij/xsz;
      xg[4] = xg[0] = xw[i+j*(xsz+1)];    yg[4] = yg[0] = yw[i+j*(xsz+1)];  
      zg[4] = zg[0] =   zw[i+j*(xsz+1)];
      xg[1] = xw[i+1+j*(xsz+1)];          yg[1] = yw[i+1+j*(xsz+1)];   
      zg[1] = zw[i+1+j*(xsz+1)];   
      xg[2] = xw[i+1+(j+1)*(xsz+1)];      yg[2] = yw[i+1+(j+1)*(xsz+1)];
      zg[2] = zw[i+1+(j+1)*(xsz+1)];
      xg[3] = xw[i+(j+1)*(xsz+1)];        yg[3] = yw[i+(j+1)*(xsz+1)];  
      zg[3] = zw[i+(j+1)*(xsz+1)]; 
// On efface derriere 
      g3->SelForeground(*cmap, mLut->ApplyFast(ijzp[l].zh));
      g3->DrawFPolygon3D(xg,yg,zg,5);
// On trace 
      if (mLAtt == PI_NotDefLineAtt)   continue;
      g3->SelForeground(fgc);
      g3->DrawPolygon3D(xg,yg,zg,5);
    }
      
  }

  delete[] ijzp;
  delete[] xw; 
  delete[] yw;
  delete[] zw;
}
