/******************************************************************************

   This file is a part of the

             COMPLEX PSEUDO-RAW SAR IMAGE PROCESSING TOOLBOX

   Copyright (C) 2017 Rémy Abergel (Remy.Abergel AT parisdescartes.fr), Loïc
   Denis (Loic.Denis AT univ-st-etienne.fr), Saïd Ladjal (Said.Ladjal AT
   telecom-paristech.fr) and Florence Tupin (Florence.Tupin AT
   telecom-paristech.fr).

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful, but
   WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>.

******************************************************************************/

/*

   Compilation (using gcc):
   ========================

   gcc -O3 modulus.c tiffread.c tiffwrite.c -lfftw3 -lm -ltiff -o modulus


   Description:
   ============

   This module can be used to compute the image modulus of a complex valued
   image.

   The output image can be saved in TIFF 8-bits (default) of in TIFF-float
   (option -f).

   Use option '-p p (in [0,100])' to saturate a ratio of 'p' percents pixels
   with minimal values.

   Use option '-q q (in [0,100])' to saturate a ratio of 'q' percents pixels
   with maximal values.

   Use option '-d d (in [0,100])' to saturate a ratio of 'd/2' percents pixels
   with minimal values and 'd/2' percents pixels with maximal values.

 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "fftw3.h"

#define ASSERT_ALLOC(cmd) if(NULL == (cmd)) { printf("Not enough memory\n"); exit(EXIT_FAILURE); }

extern int load_2pages_tiff_into_fftw_complex(fftw_complex**,int*,int*,char*,char); // see source file 'tiffread.c'
extern int save_monopage_striped_tiff_from_double(char*,double*,int,int,char,int); // see source file 'tiffwrite.c'



/* comparison function for qsort **increasing order** */
static int compare(const void *a, const void *b)
{
  double v;
  v = *((double *)a)-*((double *)b);
  return (v>0.?1:(v<0.?-1:0));
}


static void display_usage()
{
  printf("\nCompute the modulus of a complex-valued image: out = |in|.\n\n");
  printf("Usage: modulus [-p p] [-q q] [-d d] [-c] [-f] in out\n\n");
  printf("   -p p (in [0,100]) : saturate a ratio of 'p' percents pixels with minimal values\n");
  printf("   -q q (in [0,100]) : saturate a ratio of 'q' percents pixels with maximal values\n");
  printf("   -d d (in [0,100]) : saturate a ratio of 'd' percents pixels with extremal values\n");
  printf("   -c  : set output min/max values equal to 0 and 255 (affine contrast change)\n");
  printf("   -f  : save output as a tiff-float (instead of unsigned char tiff)\n");
  printf("   in  : input complex image, (stack of 2 TIFF images, first image=real part\n");
  printf("         of 'in', second image=imaginary part of 'in')\n");
  printf("   out : output image (TIFF)\n\n");
}


int main(int argc, char **argv)
{
  char *filenamein = NULL,*filenameout = NULL,*pvalue = NULL,*qvalue = NULL,*dvalue = NULL,vflag,cflag,fflag;
  double *v = NULL,*w = NULL;
  fftw_complex *u;
  double min,max,val,p,q,d,c1,c2;
  int x,y,adr,i,nx,ny,err;
  vflag = cflag = fflag = 0;

  /******************************************/
  /*            arguments parser            */
  /******************************************/

  /* main parser */
  for(i=1;i<argc;i++) {
    if (strcmp(argv[i],"-f") == 0) fflag = 1;
    else if (strcmp(argv[i],"-c") == 0) cflag = 1;
    else if (strcmp(argv[i],"-p") == 0) {
      if(i==argc-1) { display_usage(); printf("Error: option -p requires an argument\n\n"); return EXIT_FAILURE; }
      else {
	pvalue = argv[i+1]; err = sscanf(pvalue,"%lf",&p); i++;
	if(err != 1) { display_usage(); printf("Error: could not retrieve properly the argument of option -p\n\n"); return EXIT_FAILURE; }
      }
    }
    else if (strcmp(argv[i],"-q") == 0) {
      if(i==argc-1) { display_usage(); printf("Error: option -q requires an argument\n\n"); return EXIT_FAILURE; }
      else {
	qvalue = argv[i+1]; err = sscanf(qvalue,"%lf",&q); i++;
	if(err != 1) { display_usage(); printf("Error: could not retrieve properly the argument of option -q\n\n"); return EXIT_FAILURE; }
      }
    }
    else if (strcmp(argv[i],"-d") == 0) {
      if(i==argc-1) { display_usage(); printf("Error: option -d requires an argument\n\n"); return EXIT_FAILURE; }
      else {
	dvalue = argv[i+1]; err = sscanf(dvalue,"%lf",&d); i++;
	if(err != 1) { display_usage(); printf("Error: could not retrieve properly the argument of option -d\n\n"); return EXIT_FAILURE; }
      }
    }
    else if (NULL == filenamein) filenamein = argv[i];
    else if (NULL == filenameout) filenameout = argv[i];
    else { display_usage(); printf("Error: could not retrieve properly input arguments.\n"); return EXIT_FAILURE; }
  }

  /* check consistency */
  if ((NULL == filenamein) || (NULL == filenameout)) {
    display_usage();
    if (NULL == filenamein) printf("Error: input 'in' is missing\n\n");
    else if(NULL == filenameout) printf("Error: output 'out' is missing\n\n");
    return EXIT_FAILURE;
  }
  if(pvalue && dvalue) {
    display_usage();
    printf("Error: options -p and -d cannot be used together.\n");
  }
  if(qvalue && dvalue) {
    display_usage();
    printf("Error: options -q and -d cannot be used together.\n");
  }
  if(p<0||p>100){
    display_usage();
    printf("Error: you must select 'p' in [0,100].\n");
  }
  if(q<0||q>100){
    display_usage();
    printf("Error: you must select 'q' in [0,100].\n");
  }
  if(d<0||d>100){
    display_usage();
    printf("Error: you must select 'd' in [0,100].\n");
  }
  if(pvalue && qvalue && p+q > 100) {
    display_usage();
    printf("Error: when options -p and -q are used together, p+q must be smaller than 100.\n");
  }

  /*******************/
  /* Open TIFF Image */
  /*******************/

  if(EXIT_FAILURE == load_2pages_tiff_into_fftw_complex(&u,&nx,&ny,filenamein,vflag)) {
    printf("Error: could not open input '%s'\n",filenamein);
   return EXIT_FAILURE;
  }

  /*********************/
  /* Memory allocation */
  /*********************/
  ASSERT_ALLOC(v = (double*) malloc (nx*ny * sizeof(double)));
  if(pvalue||qvalue||dvalue) ASSERT_ALLOC(w = (double*) malloc (nx*ny * sizeof(double)));

  /**********/
  /* KERNEL */
  /**********/
  min = 0.; max = -1;
  // set v = |u|
  for(adr=y=0;y<ny;y++)
    for(x=0;x<nx;x++,adr++) {
      val = hypot(u[adr][0],u[adr][1]);
      v[adr] = val;
      if(w) w[adr] = val;
      if(min > val) min = val;
      if(max < val) max = val;
    }

  // Deal with saturation options
  if(w) {
    qsort((void*)w,nx*ny,sizeof(double),compare);
    if(pvalue) {
      i = (int)round(0.01*p*(double)nx*(double)ny);
      min = w[i<0?0:i>=(nx*ny)?nx*ny-1:i];
    }
    if(qvalue) {
      i = (int)round(0.01*(100.-q)*(double)nx*(double)ny);
      max = w[i<0?0:i>=(nx*ny)?nx*ny-1:i];
    }
    if(dvalue) {
      i = (int)round(0.01*0.5*d*(double)nx*(double)ny); min = w[i<0?0:i>=(nx*ny)?nx*ny-1:i];
      i = (int)round(0.01*(100.-0.5*d)*(double)nx*(double)ny); max = w[i<0?0:i>=(nx*ny)?nx*ny-1:i];
    }
    for(adr=0;adr<nx*ny;adr++) {
      val = v[adr];
      if(val > max) val = max;
      if(val < min) val = min;
      v[adr] = val;
    }
  }

  // apply contrast change (if needed)
  if(cflag && (min < max))
    for(c1=255./(max-min),c2=-c1*min,adr=0;adr<nx*ny;adr++)
      v[adr] = c1*v[adr]+c2;

  /*********************************************/
  /*               END OF KERNEL               */
  /*********************************************/

  /********************************/
  /* Write output as a TIFF Image */
  /********************************/
  if(EXIT_FAILURE == save_monopage_striped_tiff_from_double(filenameout,v,nx,ny,vflag,fflag?0:8)) {
    printf("Error: could not save output '%s'\n",filenameout);
    free(v); free(w);
    return EXIT_FAILURE;
  }

  /* free memory */
  free(v);
  if (w) free(w);

  return EXIT_SUCCESS;
}
