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

   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 sardecomp.c sardecomp_kernel.c tools.c tiffread.c tiffwrite.c -lfftw3 -lm -ltiff -o sardecomp


   Source file description:
   ========================

   This module contains a main function (with TIFF Input/Output support)
   associated to the 'SPECKLE PLUS TARGETS DECOMPOSITION' algorithm coming from
   the research paper [1] (see Algorithm 4 therein for a pseudo-code
   description)

   [1] R. Abergel, L. Denis, S. Ladjal and F. Tupin, ``Sub-pixellic Methods for
       Sidelobes Suppression and Strong Targets Extraction in Single Look
       Complex SAR Images'', preprint HAL-IMT, submitted in July 2017.

   Note: the work [1] corresponds to the revision of the manuscript entitled
   ``Sidelobes Suppression and Strong Targets Extraction in Single Look Complex
   SAR Images with Several Applications'' that we submitted in July 2017 (we
   changed the title according to one Reviewer's recommandation).

   The core functions of the algorithm are available in the source file
   'sardecomp_kernel.c'.


   Documentation of the sardecomp module:
   ======================================

   This algorithm can be used to decompose a single look complex pseudo raw
   image (u0) into two components: a speckle-dominated image (w0) and the
   point-like targets (a set of 3-uple of the type (x,y,A), such as (x,y)
   denotes the sub-pixellic coordinates of the target center, and A its complex
   amplitude).

   The algorithm can also produce the images s0 and d0. Image s0 is such as

                                u0 = w0 + s0,

   or equivalently, s0 = u0 -w0. Image s0 corresponds to a linear combination of
   sampled cardinal sine functions (each cardinal sine represents one extracted
   target). The image d0 corresponds to the linear combination of discrete
   Diracs (each discrete Dirac correspond to one extracted target whose
   sub-pixellic position was projected on the initial sampling grid), so that
   the image R(u0), defined by

                              R(u0) = w0 + d0,

   is free of sidelobe effects.

   See more detailed explanations in the source file 'sardecomp_kernel.c', as
   well as in the plain-text file 'README.org' embedded at the root directory of
   this toolbox.

   Please refer to [1] for more details and explanations about the algorithm.

*/

#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 void sardecomp(fftw_complex*,fftw_complex*,fftw_complex*,FILE*,double*,int,int,double,int,int,int,char,char); // see source file 'sardecomp_kernel.c'
extern int load_2pages_tiff_into_fftw_complex(fftw_complex**,int*,int*,char*,char); // see source file 'tiffread.c'
extern int save_2pages_striped_tiff_float_from_fftw_complex(char*,fftw_complex*,int,int,char); // see source file 'tiffwrite.c'


/* display_usage: display usage message in the user terminal */
static void display_usage()
{
  printf("\nSpeckle plus targets decomposition.\n\n");
  printf("Usage: sardecomp [-N N] [-K K] [-e e] [-n n] [-s s] [-s0 s0] [-d0 d0] [--nfaorder] [-nfaorder] [-v] [-vv] u0 w0 targets\n\n");
  printf("   -N N      : (default = 20) number of tested translations\n");
  printf("   -K K      : (default = 25) locality parameter\n");
  printf("   -e e      : (default = 1) maximal number of allowed false detections\n");
  printf("   -n n      : limit to n the number of target extractions\n");
  printf("   -s s      : sigma parameter of the NFA model (by default sigma is learnt)\n");
  printf("   -s0 s0    : output recombination of extracted cardinal sines s0 = u0 - w0\n");
  printf("   -d0 d0    : output recombination into discrete Diracs of the extracted targets\n");
  printf("   -nfaorder : force target extraction in decreasing NFA order (time consuming)\n");
  printf("   -v        : light verbose mode (display algorithm progression)\n");
  printf("   -vv       : heavy verbose mode (also display information about TIFF I/O)\n");
  printf("   u0        : input complex pseudo-raw image (a stack of 2 TIFF images)\n");
  printf("   w0        : output speckle component (a stack of 2 TIFF images)\n");
  printf("   targets   : output file containing the targets features\n\n");
  printf("Inputs/Outputs description:\n\n");
  printf("   This function is an implementation of the ''SPECKLE PLUS TARGETS\n");
  printf("   DECOMPOSITION'' algorithm coming from the research paper [1] (see Algorithm 4\n");
  printf("   therein for a pseudo-code description)\n\n");
  printf("   [1] R. Abergel, L. Denis, S. Ladjal and F. Tupin, ``Sidelobes Suppression and\n");
  printf("       Strong Targets Extraction in Single Look Complex SAR Images with Several\n");
  printf("       Applications'', preprint HAL-IMT, submitted in July 2017.\n\n");
  printf("   This algorithm can be used to decompose a single look complex pseudo raw\n");
  printf("   image (u0) into two components: a speckle-dominated image (w0) and the\n");
  printf("   point-like targets (a set of 3-uple of the type (x,y,A), such as (x,y)\n");
  printf("   denotes the sub-pixellic coordinates of the target center, and A its complex\n");
  printf("   amplitude).\n\n");
  printf("   The algorithm can also produce the images s0 and d0. Image s0 is such as\n\n");
  printf("                                u0 = w0 + s0,\n\n");
  printf("   or equivalently, s0 = u0 -w0. Image s0 corresponds to a linear combination of\n");
  printf("   sampled cardinal sine functions (each cardinal sine represents one extracted\n");
  printf("   target). The image d0 corresponds to the linear combination of discrete\n");
  printf("   Diracs (each discrete Dirac correspond to one extracted target whose\n");
  printf("   sub-pixellic position was projected on the initial sampling grid), so that\n");
  printf("   the image R(u0), defined by\n\n");
  printf("                              R(u0) = w0 + d0,\n\n");
  printf("   is free of sidelobe effects. See more details in [1].\n\n");
}


/* main function, including a simple hand-made parser */
int main(int argc, char **argv)
{

  char *filename_u0 = NULL,*filename_w0 = NULL,*filename_targets = NULL,*filename_s0 = NULL,*filename_d0 = NULL,*Nvalue = NULL,*Kvalue = NULL,*evalue = NULL,*nvalue = NULL,*svalue = NULL,nfaorder,vflag,vvflag,verbose;
  fftw_complex *w0 = NULL, *s0 = NULL, *d0 = NULL;
  double s,epsilon;
  int Ntransl,K,i,err,nx,ny,n;
  FILE *targets;
  nfaorder = vvflag = vflag = 0;

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

  /* main parser */
  for(i=1;i<argc;i++) {
    if (strcmp(argv[i],"-v") == 0) vflag = 1;
    else if (strcmp(argv[i],"-vv") == 0) vvflag = 1;
    else if (strcmp(argv[i],"-nfaorder") == 0) nfaorder = 1;
    else if (strcmp(argv[i],"-N") == 0) {
      if(i==argc-1) { display_usage(); printf("Error: option -N requires an argument\n\n"); return EXIT_FAILURE; }
      else {
	Nvalue = argv[i+1]; err = sscanf(Nvalue,"%d",&Ntransl); i++;
	if(err != 1) { display_usage(); printf("Error: could not retrieve properly the argument of option -N\n\n"); return EXIT_FAILURE; }
      }
    }
    else if (strcmp(argv[i],"-K") == 0) {
      if(i==argc-1) { display_usage(); printf("Error: option -K requires an argument\n\n"); return EXIT_FAILURE; }
      else {
	Kvalue = argv[i+1]; err = sscanf(Kvalue,"%d",&K); i++;
	if(err != 1) { display_usage(); printf("Error: could not retrieve properly the argument of option -K\n\n"); return EXIT_FAILURE; }
      }
    }
    else if (strcmp(argv[i],"-e") == 0) {
      if(i==argc-1) { display_usage(); printf("Error: option -e requires an argument\n\n"); return EXIT_FAILURE; }
      else {
	evalue = argv[i+1]; err = sscanf(evalue,"%lf",&epsilon); i++;
	if(err != 1) { display_usage(); printf("Error: could not retrieve properly the argument of option -e\n\n"); return EXIT_FAILURE; }
      }
    }
    else if (strcmp(argv[i],"-s") == 0) {
      if(i==argc-1) { display_usage(); printf("Error: option -s requires an argument\n\n"); return EXIT_FAILURE; }
      else {
	svalue = argv[i+1]; err = sscanf(svalue,"%lf",&s); i++;
	if(err != 1) { display_usage(); printf("Error: could not retrieve properly the argument of option -s\n\n"); return EXIT_FAILURE; }
      }
    }
    else if (strcmp(argv[i],"-n") == 0) {
      if(i==argc-1) { display_usage(); printf("Error: option -n requires an argument\n\n"); return EXIT_FAILURE; }
      else {
	nvalue = argv[i+1]; err = sscanf(nvalue,"%d",&n); i++;
	if(err != 1) { display_usage(); printf("Error: could not retrieve properly the argument of option -n\n\n"); return EXIT_FAILURE; }
      }
    }
    else if (strcmp(argv[i],"-s0") == 0) {
      if(i==argc-1) { display_usage(); printf("Error: option -s0 requires an argument\n\n"); return EXIT_FAILURE; }
      else { filename_s0 = argv[i+1]; i++; }
    }
    else if (strcmp(argv[i],"-d0") == 0) {
      if(i==argc-1) { display_usage(); printf("Error: option -d0 requires an argument\n\n"); return EXIT_FAILURE; }
      else { filename_d0 = argv[i+1]; i++; }
    }
    else if (NULL == filename_u0) filename_u0 = argv[i];
    else if (NULL == filename_w0) filename_w0 = argv[i];
    else if (NULL == filename_targets) filename_targets = argv[i];
    else { display_usage(); printf("Error: could not retrieve properly input arguments.\n\n"); return EXIT_FAILURE; }
  }

  /* check consistency */
  if ((NULL == filename_u0) || (NULL == filename_w0) || (NULL == filename_targets)) {
    display_usage();
    if (NULL == filename_u0) printf("Error: input 'u0' is missing\n\n");
    else if(NULL == filename_w0) printf("Error: output 'w0' is missing\n\n");
    else printf("Error: output 'targets' is missing\n\n");
    return EXIT_FAILURE;
  }
  if(evalue && epsilon<=0.){
    display_usage();
    printf("Error: input 'e' must be positive\n\n");
    return EXIT_FAILURE;
  }
  if(s<=0.){
    display_usage();
    printf("Error: input 's' must be positive\n\n");
    return EXIT_FAILURE;
  }

  /* set default values */
  if(NULL == Nvalue) Ntransl = 20;
  if(NULL == Kvalue) K = 25;
  if(NULL == evalue) epsilon = 1.;
  if(NULL == nvalue) n = -1; // we will not limit the number of target
			     // extractions (the process will stop only when no
			     // target with number of false alarms (NFA) less
			     // than epsilon can be found in the data).
  if(NULL == svalue) s = -1; // s will be learnt

  /* deal with verbose mode */
  verbose = vflag || vvflag;
  if(verbose) {
    printf("\nInput parameters\n");
    printf("================\n");
    printf("number of tested translations: N = %d\n",Ntransl);
    printf("locality parameter: K = %d\n",K);
    if(nvalue) printf("maximal number of extractions: n = %d\n",n);
    if(svalue) printf("parameter sigma of the NFA model: sigma = %g\n",s);
    else printf("parameter sigma of the NFA model will be learnt\n");
    printf("maximal number of false detections: e = %g\n\n",epsilon);
  }

  /************************************************************/
  /* Open TIFF Image & get its graylevels in double precision */
  /************************************************************/
  if(verbose) {
    printf("Load input image and allocate memory\n");
    printf("====================================\n");
  }

  // load input image u0 into w0 (avoids initialization process u0 = w0)
  if(EXIT_FAILURE == load_2pages_tiff_into_fftw_complex(&w0,&nx,&ny,filename_u0,vvflag)) {
    if(!vvflag) printf("An error occured (try option -vv to get more information)\n");
    if(verbose) printf("Free memory and exit the program.\n\n");
   return EXIT_FAILURE;
  }

  if(vflag && !vvflag) printf("done\n\n");
  if(!vflag && vvflag) printf("\n");


  /*********************************/
  /* Memory allocation for outputs */
  /*********************************/

  if(NULL == (targets = fopen(filename_targets,"w"))) {
    printf("Error: could not open file '%s'\n",filename_targets);
    return EXIT_FAILURE;
  }

  if(filename_s0) ASSERT_ALLOC(s0 = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex)));
  if(filename_d0) ASSERT_ALLOC(d0 = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex)));



  /***********************************/
  /* KERNEL: process the input image */
  /***********************************/
  sardecomp(w0,s0,d0,targets,&s,K,Ntransl,epsilon,n,nx,ny,nfaorder,verbose);

  /*******************************************/
  /* Write output images (TIFF-float Stacks) */
  /*******************************************/
  if(verbose) {
    printf("Save output(s) as Stack(s) of two TIFF images\n");
    printf("=============================================\n");
  }
  if(EXIT_FAILURE == save_2pages_striped_tiff_float_from_fftw_complex(filename_w0,w0,nx,ny,vvflag)) {
    if(!vvflag) printf("An error occured (try option -vv to get more information)\n");
    if(vflag || vvflag) printf("Free memory and exit the program.\n\n");
    free(w0);
    if(s0) free(s0);
    if(d0) free(d0);
    return EXIT_FAILURE;
  }
  if(s0 && EXIT_FAILURE == save_2pages_striped_tiff_float_from_fftw_complex(filename_s0,s0,nx,ny,vvflag)) {
    if(!vvflag) printf("An error occured (try option -vv to get more information)\n");
    if(vflag || vvflag) printf("Free memory and exit the program.\n\n");
    free(w0);
    free(s0);
    if(d0) free(d0);
    return EXIT_FAILURE;
  }
  if(d0 && EXIT_FAILURE == save_2pages_striped_tiff_float_from_fftw_complex(filename_d0,d0,nx,ny,vvflag)) {
    if(!vvflag) printf("An error occured (try option -vv to get more information)\n");
    if(vflag || vvflag) printf("Free memory and exit the program.\n\n");
    free(w0);
    if(s0) free(s0);
    free(d0);
    return EXIT_FAILURE;
  }

  /***************/
  /* free memory */
  /***************/
  fftw_free(w0);
  if(s0) fftw_free(s0);
  if(d0) fftw_free(d0);
  fclose(targets);

  if(vflag && !vvflag) printf("done\n\n");
  if(!vflag && vvflag) printf("\n");

  return EXIT_SUCCESS;
}
