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

   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/>.

   This is a MEX file for MATLAB.

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

/*                    ************************************
                      * MAIN DOCUMENTATION OF THE MODULE *
                      ************************************

   Compilation (using mex from Matlab terminal):
   =============================================

   mex -lm -lfftw3 CFLAGS="\$CFLAGS -std=c99" tools.c sardecomp_kernel.c sardecomp_mexinterface.c -output sardecomp_mexinterface


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

   This module is a mex interface for MATLAB 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 <mex.h>
#include <fftw3.h>

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'

/* usage displayer */
static void display_usage()
{
  mexPrintf("\nSpeckle plus targets decomposition (MATLAB MEX INTERFACE).\n\n");
  mexPrintf("Usage: [w0,s0,d0,sigma] = sardecomp_mexinterface(u0,epsilon,K,Ntransl,nmax,nfaorder,verbose,filename,s)\n\n");
  mexPrintf("   u0       : (matrix of complex double) input complex pseudo-raw image\n");
  mexPrintf("   epsilon  : (scalar double) maximal allowed number of false detections\n");
  mexPrintf("   K        : (scalar int32) locality parameter (must be positive)\n");
  mexPrintf("   Ntransl  : (scalar int32) number of tested translations (must be positive)\n");
  mexPrintf("   nmax     : (scalar int32) maximal number of extraction (set nmax=-1 for no limit)\n");
  mexPrintf("   nfaorder : (scalar int32) set to nonzero to extract target in decreasing nfa order\n");
  mexPrintf("   verbose  : (scalar int32) enter any non-zero value to enable verbose mode\n");
  mexPrintf("   filename : (string) path to output file containing the target features\n\n");
  mexPrintf("   s        : sigma parameter of the NFA model (set s<=0 to learn this parameter)\n");
  mexPrintf("   w0       : (matrix of complex double) output speckle component\n");
  mexPrintf("   s0       : (matrix of complex double) cardinal sines recombination s0=u0-w0\n");
  mexPrintf("   d0       : (matrix of double) recombination of the targets into discrete Diracs\n");
  mexPrintf("   sigma    : sigma parameter used in the NFA model (sigma = s if s > 0)\n\n");
  mexPrintf("Inputs/Outputs description:\n\n");
  mexPrintf("   This function is an implementation of the ''SPECKLE PLUS TARGETS\n");
  mexPrintf("   DECOMPOSITION'' algorithm coming from the research paper [1] (see Algorithm 4\n");
  mexPrintf("   therein for a pseudo-code description)\n\n");
  mexPrintf("   [1] R. Abergel, L. Denis, S. Ladjal and F. Tupin, ``Sidelobes Suppression and\n");
  mexPrintf("       Strong Targets Extraction in Single Look Complex SAR Images with Several\n");
  mexPrintf("       Applications'', preprint HAL-IMT, submitted in July 2017.\n\n");
  mexPrintf("   This algorithm can be used to decompose a single look complex pseudo raw\n");
  mexPrintf("   image (u0) into two components: a speckle-dominated image (w0) and the\n");
  mexPrintf("   point-like targets (a set of 3-uple of the type (x,y,A), such as (x,y)\n");
  mexPrintf("   denotes the sub-pixellic coordinates of the target center, and A its complex\n");
  mexPrintf("   amplitude).\n\n");
  mexPrintf("   The algorithm can also produce the images s0 and d0. Image s0 is such as\n\n");
  mexPrintf("                                u0 = w0 + s0,\n\n");
  mexPrintf("   or equivalently, s0 = u0 -w0. Image s0 corresponds to a linear combination of\n");
  mexPrintf("   sampled cardinal sine functions (each cardinal sine represents one extracted\n");
  mexPrintf("   target). The image d0 corresponds to the linear combination of discrete\n");
  mexPrintf("   Diracs (each discrete Dirac correspond to one extracted target whose\n");
  mexPrintf("   sub-pixellic position was projected on the initial sampling grid), so that\n");
  mexPrintf("   the image R(u0), defined by\n\n");
  mexPrintf("                              R(u0) = w0 + d0,\n\n");
  mexPrintf("   is free of sidelobe effects. See more details in [1].\n\n");
}

/* The gateway function */
void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{

  double *u0_re, *u0_im, *w0_re, *w0_im, *s0_re, *s0_im, *d0_re, *d0_im, epsilon, s;
  fftw_complex *w0,*s0,*d0,*T;
  size_t nx,ny;
  mwSize adr,x,y,K,Ntransl,nfaorder,verbose,nmax;
  char *filename = NULL;
  FILE *targets = NULL;

  /*****************************************/
  /* retrieve inputs and check consistency */
  /*****************************************/

  // number of inputs //
  if(nrhs != 9) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputnumber","Incorrect number of inputs");
  }

  // number of outputs //
  if(nlhs > 4) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badoutputnumber","Incorrect number of outputs");
  }

  // first input (u0 = matrix of complex double) //
  if( !mxIsDouble(prhs[0]) || !mxIsComplex(prhs[0]) ) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the first input (u0) must be of type 'matrix of complex double'.");
  }
  nx = mxGetN(prhs[0]);
  ny = mxGetM(prhs[0]);
  if(!(nx > 1 && ny > 1)) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the first input (u0) must have two dimensions (row and column vectors are not allowed).");
  }
  u0_re = mxGetPr(prhs[0]);
  u0_im = mxGetPi(prhs[0]);

  // second input (epsilon = scalar double) //
  if( !mxIsScalar(prhs[1]) || !mxIsDouble(prhs[1]) ) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the second input (epsilon) must be of type 'scalar double'.");
  }
  epsilon = (double) mxGetScalar(prhs[1]);
  if(K < 0) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the second input (espilon) must be positive.");
  }

  // third input (K = scalar 32 bits integer) //
  if( !mxIsScalar(prhs[2]) || !mxIsInt32(prhs[2]) ) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the third input (K) must be of type 'scalar int32'.");
  }
  K = (mwSize) mxGetScalar(prhs[2]);
  if(K < 1) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the third input (K) must be greater or equal to 1.");
  }

  // fourth input (Ntransl = scalar 32 bits integer) //
  if( !mxIsScalar(prhs[3]) || !mxIsInt32(prhs[3]) ) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the fourth input (Ntransl) must be of type 'scalar int32'.");
  }
  Ntransl = (mwSize) mxGetScalar(prhs[3]);
  if(Ntransl < 1) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the fourth input (Ntransl) must be greater or equal to 1.");
  }

  // fifth  input (nmax = scalar 32 bits integer) //
  if( !mxIsScalar(prhs[4]) || !mxIsInt32(prhs[4]) ) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the fifth input (nmax) must be of type 'scalar int32'.");
  }
  nmax = (mwSize) mxGetScalar(prhs[4]);
  if(nmax < 0 && nmax != -1) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the fifth input (nmax) must be positive or equal to -1.");
  }

  // sixth input (nfaorder = scalar 32 bits integer) //
  if( !mxIsScalar(prhs[5]) || !mxIsInt32(prhs[5]) ) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the sixth input (nfaorder) must be of type 'scalar int32'.");
  }
  nfaorder = (mwSize) mxGetScalar(prhs[5]);

  // seventh input (verbose = scalar 32 bits integer) //
  if( !mxIsScalar(prhs[6]) || !mxIsInt32(prhs[6]) ) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the seventh input (verbose) must be of type 'scalar int32'.");
  }
  verbose = (mwSize) mxGetScalar(prhs[6]);

  // eighth input (filename = string) //
  if( !mxIsChar(prhs[7]) ) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the eigth input (filename) must be of type 'string'.");
  }
  filename = mxArrayToString(prhs[7]);
  if(NULL == filename || NULL == (targets = fopen(filename,"w"))) {
    printf("WARNING: could not open file '%s', the target features will not be saved\n",filename);
  }
  mxFree(filename);

  // nineth input (s = scalar double) //
  if( !mxIsScalar(prhs[8]) || !mxIsDouble(prhs[8]) ) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:sardecomp_mexinterface:badinputdatatype","the nineth input (s) must be of type 'scalar double'.");
  }
  s = (double) mxGetScalar(prhs[8]);

  /**********************/
  /* memory allocation */
  /**********************/
  // output w0
  w0 = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex));
  plhs[0] = mxCreateDoubleMatrix((mwSize)ny,(mwSize)nx,mxCOMPLEX);
  // output s0
  if(nlhs >= 2) {
    s0 = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex));
    plhs[1] = mxCreateDoubleMatrix((mwSize)ny,(mwSize)nx,mxCOMPLEX);
  }
  else s0 = NULL;
  // output d0
  if(nlhs >= 3) {
    d0 = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex));
    plhs[2] = mxCreateDoubleMatrix((mwSize)ny,(mwSize)nx,mxCOMPLEX);
  }
  else d0 = NULL;


  /**********************/
  /*       KERNEL       */
  /**********************/

  // intitialize output w0 = input u0 //
  for(adr=y=0;y<ny;y++)
    for(x=0;x<nx;x++,adr++) {
      // recall that Matlab stores matrices in column-major order (as in
      // Fortran) whereas row-major order is rather used in C language.
      w0[adr][0] = u0_re[x*ny+y];
      w0[adr][1] = u0_im[x*ny+y];
    }

  // run the irregularsampling algorithm //
  if(verbose) mexPrintf("\n");
  sardecomp(w0,s0,d0,targets,&s,K,Ntransl,epsilon,nmax,nx,ny,nfaorder!=0,verbose!=0);

  // recopy output w0 into Matlab memory (I hate myself!) //
  w0_re = mxGetPr(plhs[0]);
  w0_im = mxGetPi(plhs[0]);
  for(adr=y=0;y<ny;y++)
    for(x=0;x<nx;x++,adr++) {
      w0_re[x*ny+y] = w0[adr][0];
      w0_im[x*ny+y] = w0[adr][1];
    }

  // deal with optional outputs s0 & d0 //
  if(nlhs >= 3) {
    s0_re = mxGetPr(plhs[1]);
    s0_im = mxGetPi(plhs[1]);
    d0_re = mxGetPr(plhs[2]);
    d0_im = mxGetPi(plhs[2]);
    for(adr=y=0;y<ny;y++)
      for(x=0;x<nx;x++,adr++) {
	s0_re[x*ny+y] = s0[adr][0];
	s0_im[x*ny+y] = s0[adr][1];
	d0_re[x*ny+y] = d0[adr][0];
	d0_im[x*ny+y] = d0[adr][1];
      }
  }
  else if(nlhs == 2) {
    s0_re = mxGetPr(plhs[1]);
    s0_im = mxGetPi(plhs[1]);
    for(adr=y=0;y<ny;y++)
      for(x=0;x<nx;x++,adr++) {
	s0_re[x*ny+y] = s0[adr][0];
	s0_im[x*ny+y] = s0[adr][1];
      }
  }
  if(nlhs == 4) plhs[3] = mxCreateDoubleScalar(s);

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

}
