/* Last edited on 2017-06-30 01:13:59 by stolfilocal */
#include <operacoes_b.h>
#include <assert.h>
#include <math.h>

#include <float_image_from_uint16_image.h>  
#include <float_image_to_uint16_image.h>  


//PROTOTIPOS

float decodeLo(sample_uint_t iv,sample_uint_t maxval,sample_uint_t *imin,sample_uint_t *imax,float *vmin,float *vmax);
float decodeHi(sample_uint_t iv,sample_uint_t maxval,sample_uint_t *imin,sample_uint_t *imax,float *vmin,float *vmax);
float decodeMdLo(sample_uint_t iv,sample_uint_t maxval,sample_uint_t *imin,sample_uint_t *imax,float *vmin,float *vmax);
float decodeMdHi(sample_uint_t iv,sample_uint_t maxval,sample_uint_t *imin,sample_uint_t *imax,float *vmin,float *vmax);
float decodeSdLo(sample_uint_t iv,sample_uint_t maxval,sample_uint_t *imin,sample_uint_t *imax,float *vmin,float *vmax);
float decodeSdHi(sample_uint_t iv,sample_uint_t maxval,sample_uint_t *imin,sample_uint_t *imax,float *vmin,float *vmax);
sample_uint_t codeLo(float fv,sample_uint_t maxval,float *vmin,float *vmax,int *clo,int *chi,sample_uint_t *imin,  sample_uint_t *imax);
sample_uint_t codeHi(float fv,sample_uint_t maxval,float *vmin,float *vmax,int *clo,int *chi,sample_uint_t *imin,  sample_uint_t *imax);
sample_uint_t codeMdInterval(float fv,sample_uint_t maxval,float *vmin,float *vmax,int *clo,int *chi,sample_uint_t *imin,sample_uint_t *imax);
sample_uint_t codeSdInterval(float fv,sample_uint_t maxval,float *vmin,float *vmax,int *clo,int *chi,sample_uint_t *imin,sample_uint_t *imax);

void leImagemCond
  ( char *bandir,
    char *nome,
    char *tipo,
    int res,
    char *ext,
    int round,
    float_image_t **IMG
  );
  /* Se {IMG} == NULL nao le. */ 

void leImagemIntervaloCond
  ( char *bandir,
    char *nome,
    char *tipo,
    int res,
    char *ext,
    int roundLo,
    float_image_t **LO,
    int roundHi,
    float_image_t **HI
  );
  /* Se {LO} ou {HI} == NULL nao le. */ 

//IMPLEMENTACOES

float_image_t *ex_read_image(FILE *rd, char *name, int opt)
  { bool_t close_it = FALSE;
    if (rd == NULL) { rd = open_read(name, FALSE); close_it = TRUE; }
    uint16_image_t *pim = uint16_image_read_pnm_file(rd);
    float_image_t *fim = float_image_from_uint16_image_interval(pim, opt, FALSE);
    uint16_image_free(pim);
    if (close_it) { fclose(rd); }
    return fim;
  }

void ex_write_image(FILE *wr, char *name, float_image_t *img,int opt)
  { bool_t close_it = FALSE;
    if (wr == NULL) { wr = open_write(name, FALSE); close_it = TRUE; }
    sample_uint_t maxval=255;
 /*   if (img->sz[0]>1)
    {  maxval=255; }
    else { maxval=1; }*/
    uint16_image_t *pim = float_image_to_uint16_image_interval(img, maxval ,opt,FALSE);
    bool_t forceplain = FALSE;
    bool_t verbose = FALSE;
    uint16_image_write_pnm_file(wr, pim, forceplain, verbose);
    if (close_it) { fclose(wr); }
    uint16_image_free(pim);
  }

float_image_t *getImageOpenRGB(char *name, int opt)
{
   if (strlen(name)>=1)
   {
      FILE *arq;
      arq = fopen(name,"rw");
      if (arq == NULL) { fprintf(stderr, "** falhou abertura de '%s'\n", name); assert(0); }

      float_image_t *I = ex_read_image(arq, name, opt);

      fclose(arq);
      if (I==NULL)
      {
         fprintf(stderr, "Não foi possível abrir o arquivo.\n");
      }
      return I;
   }
   return NULL;
}

void leImagens
  ( char *bandir,
    char *nome,
    int res,
    char *ext,
    float_image_t **LO,
    float_image_t **HI,
    float_image_t **MDLO,
    float_image_t **MDHI,
    float_image_t **SDLO,
    float_image_t **SDHI
  )
{
  if (res>0)
    {
      leImagemCond(bandir,nome,"Lo",res,ext,0,LO);
      leImagemCond(bandir,nome,"Hi",res,ext,1,HI);
      leImagemIntervaloCond(bandir,nome,"Md",res,ext,2,MDLO,3,MDHI);
      leImagemIntervaloCond(bandir,nome,"Sd",res,ext,4,SDLO,5,SDHI);     

    }
  else
    {
      float_image_t *IMG;
      leImagemCond(bandir,nome,"Md",res,ext,2,&IMG);
      if (LO != NULL)   { *LO = IMG; }
      if (HI != NULL)   { *HI = IMG; }
      if (MDLO != NULL) { *MDLO = IMG; }
      if (MDHI != NULL) { *MDHI = IMG; }
      if (SDLO != NULL) { *SDLO = NULL; }
      if (SDHI != NULL) { *SDHI = NULL; }
    }
}

void leImagemCond
  ( char *bandir,
    char *nome,
    char *tipo,
    int res,
    char *ext,
    int round,
    float_image_t **IMG
  )
  {  
    if (IMG != NULL) 
      { char *filename=NULL;
        char *filename = jsprintf("%s/%s/%s/R%02d.%s",bandir,nome,tipo,res,ext);
        *IMG =  getImageOpenRGB(filename,round); assert(*IMG !=NULL);
        free(filename);
      }
  }
  
void leImagemIntervaloCond
  ( char *bandir,
    char *nome,
    char *tipo,
    int res,
    char *ext,
    int roundLo,
    float_image_t **LO,
    int roundHi,
    float_image_t **HI
  )
  {  
    if ((LO != NULL) || (HI != NULL))
      { char *filename=NULL;
        char *filename = jsprintf("%s/%s/%s/R%02d.%s",bandir,nome,tipo,res,ext);
        if (LO != NULL) { *LO = getImageOpenRGB(filename,roundLo); assert(*LO !=NULL); }
        if (HI != NULL) { *HI = getImageOpenRGB(filename,roundHi); assert(*HI !=NULL); }
      }
  }

void liberaImagens(float_image_t **LO,float_image_t **HI,float_image_t **MDLO,float_image_t **MDHI,float_image_t **SDLO,float_image_t **SDHI)
{
    if ((*HI!=NULL)&&(*HI!=*MDLO)){float_image_free(*HI);}
    if ((*LO!=NULL)&&(*LO!=*MDLO)){float_image_free(*LO);}
    if ((*MDHI!=NULL)&&(*MDLO!=*MDHI)){float_image_free(*MDHI);}
    if (*MDLO!=NULL){float_image_free(*MDLO);}
    if (*SDHI!=NULL){float_image_free(*SDHI);}
    if (*SDLO!=NULL){float_image_free(*SDLO);}
    *HI=NULL;
    *LO=NULL;
    *MDHI=NULL;
    *MDLO=NULL;
    *SDHI=NULL;
    *SDLO=NULL;
}

//operações de gravação e leitura

float decodeLo
(  sample_uint_t iv, 
   sample_uint_t maxval,
   sample_uint_t *imin, 
   sample_uint_t *imax, 
   float *vmin,
   float *vmax
)
  /* Convert integer intensity to float, kep stats: */
{
  if ((imin != NULL) && (iv < (*imin))) { (*imin) = iv; }
  if ((imax != NULL) && (iv > (*imax))) { (*imax) = iv; }
  /* Convert integer {iv} to double {rv} accounting for quantization error in writing: */
  double rv;
  if (iv == 0)
    { rv = 0.0; }
  else
    { rv = ((double)iv - 1)/((double)maxval); }
  /* Convert {rv} to float {fv} and adjust to account for arithmetic rounding errors: */
  float fv = rv - 1.0e-6;
  if ((iv != 0) && (fv <= 0)) { fv = 1.0e-100; }
  if (fv > 1.0) { fv = 1.0; }
  if ((vmin != NULL) && (fv < (*vmin))) { (*vmin) = fv; }
  if ((vmax != NULL) && (fv > (*vmax))) { (*vmax) = fv; }
  return fv;
}

float decodeHi
(  sample_uint_t iv, 
   sample_uint_t maxval,
   sample_uint_t *imin, 
   sample_uint_t *imax, 
   float *vmin,
   float *vmax
)
  /* Convert integer intensity to float, kep stats: */
{
  if ((imin != NULL) && (iv < (*imin))) { (*imin) = iv; }
  if ((imax != NULL) && (iv > (*imax))) { (*imax) = iv; }
  /* Convert integer {iv} to double {rv} accounting for quantization error in writing: */
  double rv;
  if (iv == 0)
    { rv = 0.0; }
  else
    { rv = ((double)iv)/((double)maxval); }
  /* Convert {rv} to float {fv} and adjust to account for arithmetic rounding errors: */
  float fv = rv + 1.0e-6;
  if ((iv != 0) && (fv <= 0)) { fv = 1.0e-100; }
  if (fv > 1.0) { fv = 1.0; }
  if ((vmin != NULL) && (fv < (*vmin))) { (*vmin) = fv; }
  if ((vmax != NULL) && (fv > (*vmax))) { (*vmax) = fv; }
  return fv;
}

float decodeMdLo
( sample_uint_t iv, 
   sample_uint_t maxval,
   sample_uint_t *imin, 
   sample_uint_t *imax, 
   float *vmin,
   float *vmax
)
/* Convert integer intensity to float, kep stats: */
{ 
  if ((imin != NULL) && (iv < (*imin))) { (*imin) = iv; }
  if ((imax != NULL) && (iv > (*imax))) { (*imax) = iv; }
  /* Convert integer {iv} to double {rv} accounting for quantization error in writing: */
  double rv;
  if (iv == 0)
    { 
      rv = 0.0;
    }
  else
    { 
      rv = ((double)iv)/((double)maxval);  
    }
  /* Convert {rv} to float {fv} and adjust to account for arithmetic rounding errors: */
  float fv = rv;
  if ((iv != 0) && (fv <= 0)) { fv = 1.0e-100; }
  if (fv > 1.0) { fv = 1.0; }
  if ((vmin != NULL) && (fv < (*vmin))) { (*vmin) = fv; }
  if ((vmax != NULL) && (fv > (*vmax))) { (*vmax) = fv; }
  return fv;
}

float decodeMdHi
( sample_uint_t iv, 
   sample_uint_t maxval,
   sample_uint_t *imin, 
   sample_uint_t *imax, 
   float *vmin,
   float *vmax
)
/* Convert integer intensity to float, kep stats: */
{ 
  if ((imin != NULL) && (iv < (*imin))) { (*imin) = iv; }
  if ((imax != NULL) && (iv > (*imax))) { (*imax) = iv; }
  /* Convert integer {iv} to double {rv} accounting for quantization error in writing: */
  double rv;
  if (iv == 0)
    { 
      rv = 0.0;
    }
  else
    { 
      rv = ((double)iv + 1)/((double)maxval); 
    }
  /* Convert {rv} to float {fv} and adjust to account for arithmetic rounding errors: */
  float fv = rv;
  if ((iv != 0) && (fv <= 0)) { fv = 1.0e-100; }
  if (fv > 1.0) { fv = 1.0; }
  if ((vmin != NULL) && (fv < (*vmin))) { (*vmin) = fv; }
  if ((vmax != NULL) && (fv > (*vmax))) { (*vmax) = fv; }
  return fv;
}

float decodeSdLo
( sample_uint_t iv, 
  sample_uint_t maxval,
  sample_uint_t *imin, 
  sample_uint_t *imax, 
  float *vmin,
  float *vmax
)
/* Convert integer intensity to float, kep stats: */
{ 
  if ((imin != NULL) && (iv < (*imin))) { (*imin) = iv; }
  if ((imax != NULL) && (iv > (*imax))) { (*imax) = iv; }
  /* Convert integer {iv} to double {rv} accounting for quantization error in writing: */
  double rv;
  if (iv == 0)
    { 
      rv = 0.0; 
    }
  else
    { 
      rv = (((double)iv)/((double)maxval))*0.5;
    }
  /* Convert {rv} to float {fv} and adjust to account for arithmetic rounding errors: */
  float fv = rv;
  if ((iv != 0) && (fv <= 0)) { fv = 1.0e-100; }
  if (fv > 0.5) { fv = 0.5; }
  if ((vmin != NULL) && (fv < (*vmin))) { (*vmin) = fv; }
  if ((vmax != NULL) && (fv > (*vmax))) { (*vmax) = fv; }
  return fv;
}

float decodeSdHi
( sample_uint_t iv, 
  sample_uint_t maxval,
  sample_uint_t *imin, 
  sample_uint_t *imax, 
  float *vmin,
  float *vmax
)
/* Convert integer intensity to float, kep stats: */
{ 
  if ((imin != NULL) && (iv < (*imin))) { (*imin) = iv; }
  if ((imax != NULL) && (iv > (*imax))) { (*imax) = iv; }
  /* Convert integer {iv} to double {rv} accounting for quantization error in writing: */
  double rv;
  if (iv == 0)
    { 
      rv = 0.0;
    }
  else
    { 
      rv = (((double)iv + 1)/((double)maxval))*0.5;
    }
  /* Convert {rv} to float {fv} and adjust to account for arithmetic rounding errors: */
  float fv = rv;

  if ((iv != 0) && (fv <= 0)) { fv = 1.0e-100; }
  if (fv > 0.5) { fv = 0.5; }
  if ((vmin != NULL) && (fv < (*vmin))) { (*vmin) = fv; }
  if ((vmax != NULL) && (fv > (*vmax))) { (*vmax) = fv; }
  return fv;
}

sample_uint_t codeLo
( float fv, 
  sample_uint_t maxval,
  float *vmin,
  float *vmax, 
  int *clo,
  int *chi,
  sample_uint_t *imin, 
  sample_uint_t *imax
)
{ 
  demand(! isnan(fv), "{fv} is NaN"); 
  if ((vmin != NULL) && (fv < (*vmin))) { (*vmin) = fv; }
  if ((vmax != NULL) && (fv > (*vmax))) { (*vmax) = fv; }

  if (fv < 0.0) { fv = 0.0;  if (clo != NULL) { (*clo)++; } }
  if (fv > 1.0) { fv = 1.0;  if (chi != NULL) { (*chi)++; } }
  sample_uint_t zv = (int)ceil(fv*maxval);
  demand((zv >= 0) && (zv <= (int)maxval), "bad {zv}");
  sample_uint_t iv = zv;
  if ((imin != NULL) && (iv < (*imin))) { (*imin) = iv; }
  if ((imax != NULL) && (iv > (*imax))) { (*imax) = iv; }
  return iv;
}

sample_uint_t codeHi
( float fv, 
  sample_uint_t maxval,
  float *vmin,
  float *vmax, 
  int *clo,
  int *chi,
  sample_uint_t *imin, 
  sample_uint_t *imax
)
{ 
  demand(! isnan(fv), "{fv} is NaN"); 
  if ((vmin != NULL) && (fv < (*vmin))) { (*vmin) = fv; }
  if ((vmax != NULL) && (fv > (*vmax))) { (*vmax) = fv; }
 
  if (fv < 0.0) { fv = 0.0;  if (clo != NULL) { (*clo)++; } }
  if (fv > 1.0) { fv = 1.0;  if (chi != NULL) { (*chi)++; } }
  sample_uint_t zv = (int)ceil(fv*maxval);
  demand((zv >= 0) && (zv <= (int)maxval), "bad {zv}");
  sample_uint_t iv = zv;
  if ((imin != NULL) && (iv < (*imin))) { (*imin) = iv; }
  if ((imax != NULL) && (iv > (*imax))) { (*imax) = iv; }
  return iv;
}

sample_uint_t codeMdInterval
  ( float fv, 
    sample_uint_t maxval,
    float *vmin,
    float *vmax, 
    int *clo,
    int *chi,
    sample_uint_t *imin, 
    sample_uint_t *imax
  )
  { demand(! isnan(fv), "{fv} is NaN"); 
    if ((vmin != NULL) && (fv < (*vmin))) { (*vmin) = fv; }
    if ((vmax != NULL) && (fv > (*vmax))) { (*vmax) = fv; }

    if (fv < 0.0) { fv = 0.0;  if (clo != NULL) { (*clo)++; } }
    if (fv > 1.0) { fv = 1.0;  if (chi != NULL) { (*chi)++; } }
    sample_uint_t zv = (int)ceil(((fv==1)?maxval-1:maxval)*fv);
    demand((zv >= 0) && (zv <= (int)maxval), "bad {zv}");
    sample_uint_t iv = zv;
    if ((imin != NULL) && (iv < (*imin))) { (*imin) = iv; }
    if ((imax != NULL) && (iv > (*imax))) { (*imax) = iv; }
    return iv;
  }

sample_uint_t codeSdInterval
  ( float fv, 
    sample_uint_t maxval,
    float *vmin,
    float *vmax, 
    int *clo,
    int *chi,
    sample_uint_t *imin, 
    sample_uint_t *imax
  )
  { demand(! isnan(fv), "{fv} is NaN"); 
    if ((vmin != NULL) && (fv < (*vmin))) { (*vmin) = fv; }
    if ((vmax != NULL) && (fv > (*vmax))) { (*vmax) = fv; }

    if (fv < 0.0) { fv = 0.0;  if (clo != NULL) { (*clo)++; } }
    if (fv > 0.5) { fv = 0.5;  if (chi != NULL) { (*chi)++; } }
  //  sample_uint_t zv = (int)ceil(((fv==0.5)?maxval-1:maxval)*(fv/0.5));
    sample_uint_t zv = (int)ceil(maxval*(fv/0.5));
    demand((zv >= 0) && (zv <= (int)maxval), "bad {zv}");
    sample_uint_t iv = zv;
    if ((imin != NULL) && (iv < (*imin))) { (*imin) = iv; }
    if ((imax != NULL) && (iv > (*imax))) { (*imax) = iv; }
    return iv;
  }



float_image_t *float_image_from_uint16_image_interval
  ( uint16_image_t *iim, 
    sign_t opt, 
    bool_t verbose
  )
  { 
    /* Get image dimensions: */
    int NX = iim->cols;
    int NY = iim->rows;
    
    /* Channel counts: */
    int chns = iim->chns; /* Num of channels. */
    
    /* Allocate float image: */
    float_image_t *fim = float_image_new(chns, NX, NY);
    
    /* Max sample value in integer image: */
    uint16_t maxval = iim->maxval;
    
    /* Input and output range registers: */
    sample_uint_t imin[chns], imax[chns]; /* Input range registers. */ 
    float vmin[chns], vmax[chns];         /* Output range registers. */ 
    int c; /* Channel index. */
    for (c = 0; c < chns; c++) 
      { imin[c] = maxval;
        imax[c] = 0;
        vmin[c] = +INF;
        vmax[c] = -INF;
      }
    
    /* Convert pixels, keep statistics: */
    int x, y;
    for(y = 0; y < NY; y++)
      { int pgmy = NY - 1 - y;
        uint16_t *prow = iim->smp[pgmy];
        for(x = 0; x < NX; x++)
          { for (c = 0; c < chns; c++)
              { /* Convert int sample {*prow} to float {v}, store, keep stats: */
                uint16_t ismp = (*prow);
                float fsmp = 0;
                if (opt==0)
                   fsmp = decodeLo(ismp,maxval,&(imin[c]),&(imax[c]),&(vmin[c]),&(vmax[c]));
                else if (opt==1)
                   fsmp = decodeHi(ismp,maxval,&(imin[c]),&(imax[c]),&(vmin[c]),&(vmax[c]));
                else if (opt==2)
                   fsmp = decodeMdLo(ismp,maxval,&(imin[c]),&(imax[c]),&(vmin[c]),&(vmax[c]));
                else if (opt==3)
                   fsmp = decodeMdHi(ismp,maxval,&(imin[c]),&(imax[c]),&(vmin[c]),&(vmax[c]));
                else if (opt==4)
                   fsmp = decodeSdLo(ismp,maxval,&(imin[c]),&(imax[c]),&(vmin[c]),&(vmax[c]));
                else if (opt==5)
                   fsmp = decodeSdHi(ismp,maxval,&(imin[c]),&(imax[c]),&(vmin[c]),&(vmax[c]));

                float_image_set_sample(fim, c, x, y, fsmp);
                prow++;
              }
          }
      }
    
    if (verbose) 
      { /* Print statistics: */
        long int NPIX = ((long int)NX)*((long int)NY);
        fprintf(stderr, "  %ld pixels in PNM image\n", NPIX);
        if (NPIX > 0)
          { for (c = 0; c < chns; c++)
              { double loc = 0.0;
                double hic = 1.0;
                sample_conv_print_floatize_stats
                  ( c, c, imin[c], imax[c], maxval, loc, hic, vmin[c], vmax[c] );
              }
          }
      }
    
    return fim;
  }

uint16_image_t *float_image_to_uint16_image_interval
  ( float_image_t *fim, 
    sample_uint_t maxval,
    int opt,
    bool_t verbose
  )
  { /* Get image dimensions: */
    int NC = fim->sz[0];
    int NX = fim->sz[1];
    int NY = fim->sz[2];
    
    /* Channel counts: */
    int fchns = NC; /* Num channels in float image. */
    int ichns = NC; /* Num channels in integer image. */
    int chns = NC;
    /* Allocate PGM/PPM image: */
    uint16_image_t *iim = uint16_image_new(NX, NY, ichns);
    
    /* Set max sample value in integer image: */
    iim->maxval = maxval;
    
    /* Channel indexing variables: */
    int k; /* Channel of integer image. */
    int c; /* Channel of float image. */
    
    /* Input and output range registers: */
    float vmin[ichns], vmax[ichns];         /* Float pixel range. */
    sample_uint_t imin[ichns], imax[ichns]; /* Int pixel range. */
    int clo[ichns], chi[ichns];             /* Counts of lo-clipped and hi-clipped pixels. */
    for (k = 0; k < ichns; k++) 
      { clo[k] = chi[k] = 0;
        vmin[k] = +INF;
        vmax[k] = -INF; 
        imin[k] = maxval;
        imax[k] = 0;  
      }
    
    /* Convert pixels, store in {iim}, keep statistics: */
    int x, y;
    for(y = 0; y < NY; y++)
      { int ppmy = NY-1-y;
        uint16_t *prow = iim->smp[ppmy];
        for(x = 0; x < NX; x++)
          { /* Convert float pixel {fpxy[c..c+2]} to integer pixel {ipxy[0..2]}, keep stats: */
            for (k = 0; k < ichns; k++)
              { c = k;
                float v = ((c < 0) || (c >= fchns) ? 0.0 : float_image_get_sample(fim, c, x, y));

                if (opt==0)
                   (*prow) = codeLo(v,maxval,&(vmin[k]),&(vmax[k]),&(clo[k]),&(chi[k]),&(imin[k]),&(imax[k]));
                else if (opt==1)
                   (*prow) = codeHi(v,maxval,&(vmin[k]),&(vmax[k]),&(clo[k]),&(chi[k]),&(imin[k]),&(imax[k]));
                else if ((opt==2) || (opt==3))
                   (*prow) = codeMdInterval(v,maxval,&(vmin[k]),&(vmax[k]),&(clo[k]),&(chi[k]),&(imin[k]),&(imax[k]));              
                else if ((opt==4) || (opt==5))
                   (*prow) = codeSdInterval(v,maxval,&(vmin[k]),&(vmax[k]),&(clo[k]),&(chi[k]),&(imin[k]),&(imax[k]));               
             
                prow++;
              }
          }
      }
    
    if (verbose)
      { /* Print statistics: */
        long int NPIX = ((long int)NX)*((long int)NY);
        fprintf(stderr, "  %ld pixels in float image\n", NPIX);
        if (NPIX > 0)
          { for (k = 0; k < chns; k++)
              { double lok = 0.0;
                double hik = 1.0;
                c = k;
                sample_conv_print_quantize_stats
                  ( c, k, vmin[k], vmax[k], lok, hik, clo[k], chi[k], maxval, imin[k], imax[k]);
              }
          }
      }
    
    return iim;
  }