//     Classes ArrayAdapter           R. Ansari  05/98
// LAL (Orsay) / IN2P3-CNRS  DAPNIA/SPP (Saclay) / CEA

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

// #include "ctimer.h"   Pour mesure TCPU 
#include "parradapter.h"

/* --Methode-- */
P1DArrayAdapter::P1DArrayAdapter(int sz)
{
mSize=sz;
DefineXCoordinate(0, 1.);
}

/* --Methode-- */
P1DArrayAdapter::~P1DArrayAdapter()
{
}

/* --Methode-- */
void P1DArrayAdapter::DefineXCoordinate(double x0, double dx)
{
mOx = x0;  mDx = dx;
}

/* --Methode-- */
double P1DArrayAdapter::X(int i)
{
return(mOx+mDx*i);
}

/* --Methode-- */
P2DArrayAdapter::P2DArrayAdapter(int sx, int sy)
{
sX = sx; sY = sy; 
sgnX = sgnY = 1;
offX = offY = 0;
invX = invY = eXY = false;
DefineXYCoordinates(0.,0.,1.,1.);
}

/* --Methode-- */
P2DArrayAdapter::~P2DArrayAdapter()
{
}

/* --Methode-- */
void P2DArrayAdapter::ConfigureAxes(bool invx, bool invy, bool exy, int sx, int sy)
{
invX = invx;  invY = invy;  eXY = exy; 
if ( (sx > 0) && (sy > 0) )   { sX = sx; sY = sy; }
if (invX) { sgnX = -1;  offX = sX-1; } 
else { sgnX = 1; offX = 0; }
if (invY) { sgnY = -1;  offY = sY-1; } 
else { sgnY = 1; offY = 0; }
}

/* --Methode-- */
void P2DArrayAdapter::DefineXYCoordinates(double x0, double y0, double dx, double dy)
{
mOx = x0;   mOy = y0;
mDx = dx;   mDy = dy;
}

/* --Methode-- */
void P2DArrayAdapter::XYfromxy(int ix, int iy, double& x, double& y)
{
x = mOx+ix*mDx;
y = mOy+iy*mDy; 
}

/* --Methode-- */
int P2DArrayAdapter::CheckDyn(double& min, double& max, double& moy, double& sig, int& nbnul, int& nbsat, int vit)
{
// TIMEF ; 
double m = 0.;
double s = 0.;
double dv; 

int n9 = 0;
int n0 = 0;

if (min >= max) { min = -9.e19; max = 9.e19; }
double minpix = max;
double maxpix = min; 
double pixv;
// printf("DEBUG_CHECKDYN-IN %d - %g %g ", vit, min, max);

int n = 0;
if (vit < 1) vit = 1;
if (sX*sY/vit < 1000) {
  vit = sX*sY/1000;
  if (vit < 1) vit = 1;
}

// printf(" -> vit=%d \n", vit);
for(int j=0; j<sY; j++) 
  for(int i=0; i<sX; i++)   {
    if (++n % vit != 0)  continue;
    pixv = Value(i,j);
    if (pixv <= max && pixv >= min) {  
      if (pixv < minpix)  minpix = pixv;
      if (pixv > maxpix)  maxpix = pixv; 
      dv = (double)pixv;
      m += dv;  s += (dv*dv);
      }
    else  {  
      if (pixv > max)  n9++;
      if (pixv < min)  n0++; 
    }
  }

nbnul = n0*vit;  nbsat = n9*vit;
min = minpix;  max = maxpix; 

n -= (n0+n9);
if (n > 0)
  { dv = (double)n;  m /= dv;  s = s/dv - m*m;
  if (s>0.)  s = sqrt(s);  else s = 0.;
  moy = m;  sig = s; }
else { moy = 0.;  sig = -1.; }

// printf("DEBUG_CHECKDYN: %d %g %g - %g %g \n", n, min, max, moy, sig);
return(n);
}


/* --Methode-- */
void P2DArrayAdapter::ComputeLut_PicHisto(double& min, double& max, int nbsig, double fracmax)
{
// TIMEF ; 
int MINNBIN = 64;
int MAXNBIN = 1024;

int nbin = sX*sY/10;
nbin = (nbin < MINNBIN) ? MINNBIN : ((nbin > MAXNBIN) ? MAXNBIN : nbin );
int vitesse = sX*sY / 20000;
vitesse = (vitesse < 1) ? 1 : ((vitesse > 25) ? 25 : vitesse );


//  Variable decrivant l'histogramme 
double minhis, maxhis, binwidth;
int under,over,tothis;
int *phis;
int *q;

double MINDYNAMIC = 1.e-8;          /* Gamme dynamique minimale (MaxAff-MinAff)*/
double FRACPIC = 0.6065;            /*  Valeur de exp(-0.5)   */



if ( min >= (max-MINDYNAMIC) ) { 
  int nbnul, nbsat, nok;
  double moy, sig; 
  nok = CheckDyn(min, max, moy, sig, nbnul, nbsat, vitesse);
  if (nok < 1) { min=0.; max = 1.;  return; }
  }
if ( min >= (max-MINDYNAMIC) )  max = min+MINDYNAMIC;
minhis = min;
maxhis = max;
binwidth = (maxhis-minhis) / nbin;
if (binwidth < MINDYNAMIC/MINNBIN) { 
  binwidth = MINDYNAMIC/MINNBIN;
  maxhis = nbin*binwidth + minhis;
  }
// printf("DEBUG nbin minhis maxhis binwidth = %d %g %g %g \n", nbin, minhis, maxhis, binwidth);

/* allocation dynamique de memoire et mise a zero de l'histogramme */
phis = new int_4[nbin];
for ( q=phis ; q<phis+nbin ; *q++ = 0 )  ;
under = over = tothis = 0;

//  Boucle sur les pixels   
//  Pour aller plus vite, on prend un pixel sur vitesse   
int kk = 0;
int bin;
for(int j=0; j<sY; j++)
  for(int i=0; i<sX; i++) {
    if ( (++kk % vitesse) != 0) continue;
   /* remplissage de l'histogramme */
    bin = (int)((Value(i, j) - minhis) / binwidth) ;
    if ( bin < 0 )  under++;
    else  {
      if (bin < nbin) { (*(phis+bin))++; tothis++; }
      else   over++; 
      }
    }     /*  Fin de for : Boucle sur les pixels   */

/*  On verifie qu'il y a eu une entree ds l'histo */
if (tothis <= 0)
  { delete[] phis;  return; }
  
/* recherche du maximum */
int binmax = -1 ; 
int hmax = -1 ; 
q = phis;
for (kk = 0 ; kk < nbin ;  kk++, q++ )
  if ( *(q) > hmax ) {hmax = *(q) ; binmax = kk ;}

if ( (binmax < 0) || (binmax >= nbin) )
  { delete[] phis;  return; }


/*   On va chercher le sigma du fond  SigmaLeft, SigmaRight */
int hmin = (int) ( (double)(hmax) * FRACPIC ) ;

int sigr = 0;
q = phis+binmax;
while ( ((*q) > hmin) && (sigr < (nbin-binmax)) )  { sigr++;  q++; }
int sigl = 0;
q = phis+binmax;
while ( ((*q) > hmin) && (sigl < binmax) )  { sigl++; q--; }


if (nbsig == 0) nbsig = 1;
if (nbsig < 0)  {  /*   On considere le cas ou NbSigmas est negatif   */
  if ( (bin = binmax+nbsig*sigl) < 0 )  bin = 0;
  min = minhis + (double)bin*binwidth;
  if ( (bin = binmax-nbsig*sigr) >= nbin ) bin = nbin-1;
  max = minhis + (double)bin*binwidth;
  if (max < (min+MINDYNAMIC))  max = min+MINDYNAMIC;
  delete[] phis;  return; 
  }

// NbSigma positif 
if ( (bin = binmax+nbsig*sigl) >= nbin ) {
  min = minhis + binmax*binwidth;
  max = maxhis;
  if (max < (min+MINDYNAMIC))  max = min+MINDYNAMIC;
  delete[] phis;  return;
  }

min = minhis + bin*binwidth;

/*   Calcul du nb de pixel entre MinAff et ValSat  */
int nbpix = 0;
for( q = phis+bin ; q < phis+nbin ; q++)  nbpix += (*q);
nbpix = (int_4)((double)nbpix * fracmax);
int nbpix2 = 0; q = phis+bin; 
int k = bin;
while( (nbpix2 < nbpix) && (k < nbin) )  
  { nbpix2 += (*q); k++; q++; }

max = minhis + (double)k*binwidth;
delete[] phis;
if (max < (min+MINDYNAMIC))  max = min+MINDYNAMIC;

return;
}


/* --Methode-- */
unsigned char * P2DArrayAdapter::ComputePixmap(LUT* lut, int ofx, int ofy, int zm, 
                int xwsz, int ywsz, unsigned char * opix, int * oxsp, int * oysp)
{
// TIMEF ; 
unsigned char *pix, *pp, ucp;
int zmm=1;
int npx, npy, nppx, nppy, nszp, oszp;
int i,j,k,l;
double vpx;


if ( (XSize() <= 0) || (YSize() <= 0) || (lut == NULL) ) 
  { *oxsp = *oysp = 0; return(NULL); }

npx = XSize()-ofx;
npy = YSize()-ofy;

if ((zm == 0) || (zm == -1))  zm = 1;

if (zm >= 1)
  {
  if (npx > (xwsz/zm))  npx = xwsz/zm; 
  if (npy > (ywsz/zm))  npy = ywsz/zm;
  nppx = npx*zm;   nppy = npy*zm;
  }
else
  {
  zmm = -zm;
  nppx = npx/zmm;   nppy = npy/zmm;
  if (nppx > xwsz)  nppx = xwsz; 
  if (nppy > ywsz)  nppy = ywsz;
  npx = nppx*zmm;   npy = nppy*zmm; 
  }

oszp = (*oxsp)*(*oysp);
nszp = nppx*nppy;
if ( oszp != nszp )
  {
  if (opix)  delete[] opix;
  pix = new unsigned char[nppx*nppy];
  if (pix == NULL)  { *oxsp = *oysp = 0; ; return(NULL); }
  }
else  pix = opix;
*oxsp = nppx; *oysp = nppy;

/*
printf("Debug_ComputePixmap %d %d %d (%d-%d %d-%d)\n", 
      ofx, ofy, zm, npx,nppx, npy,nppy);
*/

pp = pix;
if (eXY) { // Echange Axe X,Y 
  if (zm == 1)  {  // Pas de zoom  
    for(j=ofy; j<ofy+npy; j++)
      for(i=ofx; i<ofx+npx; i++) 
        { *pp = (unsigned char) lut->ApplyFast( 
                this->Value(j*sgnX+offX, i*sgnY+offY));  
           pp++; }
    }
  else if (zm < -1) {       // Compression 
    double fv = (double)(zmm*zmm); 
    for(j=ofy; j<ofy+npy; j+=zmm)
      for(i=ofx; i<ofx+npx; i+=zmm) {
        vpx = 0;
        for(l=0; l<zmm; l++)
          for(k=0; k<zmm; k++)
            vpx += this->Value( (j+l)*sgnX+offX, (i+k)*sgnY+offY);
        *pp = (unsigned char) lut->ApplyFast(vpx/fv); 
        pp++;
        }
    }
  else  {      // Agrandissement
    for(j=ofy; j<ofy+npy; j++)
      for(i=ofx; i<ofx+npx; i++)  {
        ucp = (unsigned char) lut->ApplyFast(this->Value(j*sgnX+offX, i*sgnY+offY));
        for(l=0; l<zm; l++)  {
          pp = pix+((j-ofy)*zm+l)*nppx+((i-ofx)*zm);
          for(k=0; k<zm; k++)  { *pp = ucp; pp++; }  
	  }
        }
    }
  }     // Fin du cas avec echange d'axe X/Y 
else {
  if (zm == 1)  {  // Pas de zoom  
    for(j=ofy; j<ofy+npy; j++)
      for(i=ofx; i<ofx+npx; i++) 
        { *pp = (unsigned char) lut->ApplyFast( 
                this->Value(i*sgnX+offX, j*sgnY+offY));  
           pp++; }
    }
  else if (zm < -1) {       // Compression 
    double fv = (double)(zmm*zmm); 
    for(j=ofy; j<ofy+npy; j+=zmm)
      for(i=ofx; i<ofx+npx; i+=zmm) {
        vpx = 0;
        for(l=0; l<zmm; l++)
          for(k=0; k<zmm; k++)
            vpx += this->Value( (i+k)*sgnX+offX, (j+l)*sgnY+offY);
        *pp = (unsigned char) lut->ApplyFast(vpx/fv); 
        pp++;
        }
    }
  else  {      // Agrandissement
    for(j=ofy; j<ofy+npy; j++)
      for(i=ofx; i<ofx+npx; i++)  {
        ucp = (unsigned char) lut->ApplyFast(this->Value(i*sgnX+offX, j*sgnY+offY));
        for(l=0; l<zm; l++)  {
          pp = pix+((j-ofy)*zm+l)*nppx+((i-ofx)*zm);
          for(k=0; k<zm; k++)  { *pp = ucp; pp++; }  
	  }
        }
    }
  }
return(pix);
}

