// Module PI : Peida Interactive     PISurfaceDrawer
// Traceur de surface 3D  - R.Ansari, C.Magneville 98 
// LAL (Orsay) / IN2P3-CNRS  DAPNIA/SPP (Saclay) / CEA

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

#define PERC_GARDE 0.05

//++
// Class	PISurfaceDrawer
// Lib		PI
// include	pisurfdr.h
//
//	Classe de traceur d'une surface Z = f(X,Y)
//--
//++
// Links	Parents
// PIDrawer3D
//--
//++
// Links	Voir aussi
// P2DArrayAdapter
//--
//++
// Titre	Constructeur
//--
//++
//  PISurfaceDrawer(P2DArrayAdapter* arr, bool asxy=false, bool asz=true, bool ad = true)
//	Constructeur. Tracé d'une surface définie par le tableau 2D "arr".
//	Si asxy == true : Application d'un facteur d'échelle pour égaliser 
//	le rapport d'aspect des deux axes X et Y
//	Si asz == true : Application d'un facteur d'échelle pour égaliser 
//	le rapport d'aspect de l'axe Z par rapport a XY.
//--

/* --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);

  aScxy = asxy;
  aScz = asz;

  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;
//  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)*PERC_GARDE;
  mLut = new LUT(zmin-dzz*2, zmax-dzz, 256);

  SetName("SurfDrw");
}

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

/* --Methode-- */
void PISurfaceDrawer::UpdateLimits()
{
  if (mArr == NULL)  return;
  ComputeScaleFactor(mXmin, mXmax, mYmin, mYmax, mZmin, mZmax, 
		     aScxy, aScz, mScy, mScz);
  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);  

  Set3DBox(xmin, xmax, ymin, ymax, zmin-PERC_GARDE*(zmax-zmin), 
	   zmax+PERC_GARDE*(zmax-zmin), mScx, mScy, mScz);

//  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(CC_HAS_RTTI_SUPPORT)
  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-PERC_GARDE*(mZmax-mZmin),mZmax+PERC_GARDE*(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(CC_HAS_RTTI_SUPPORT)
  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 = ScaleX(mXmax-mXmin)/xsz; 
  dy = ScaleY(mYmax-mYmin)/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 = ScaleX(x); x += dx/2.;
      y = ScaleY(y); y += dy/2.;
      ijzp[k].zh = mArr->Value(i,j);
      z = ScaleZ(ijzp[k].zh);
      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 = ScaleZ(z/4.);        
//      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 = ScaleX(x); 
     y = ScaleY(y); 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 = ScaleZ(z/4.);
//     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 = ScaleX(x); x += dx;
     y = ScaleY(y); 
     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 = ScaleZ(z/4.);
//     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 = ScaleX(x); x += dx;
  y = ScaleY(y); y += dy;
  z = mArr->Value(i-1,j-1) + mArr->Value(i1,j-1) + mArr->Value(i1,j1) +  mArr->Value(i-1,j1);
  z = ScaleZ(z/4.);
//  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 mfcol = GetGraphicAtt().GetFgColor();
  PIColors fgc = mfcol;


  if (GetGraphicAtt().GetColMapId() == CMAP_OTHER) {   // Trace avec des lignes
    if (GetGraphicAtt().GetLineAtt() == 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
    bool revcm = false;
    PIColorMap* cmap = new PIColorMap(GetGraphicAtt().GetColMapId(revcm));
    cmap->ReverseColorIndex(revcm);
    CMapId mcmapid = GetGraphicAtt().GetColMapId();
    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 (GetGraphicAtt().GetLineAtt() == PI_NotDefLineAtt)   continue;
      g3->SelForeground(fgc);
      g3->DrawPolygon3D(xg,yg,zg,5);
    }
      
  }

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

