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

   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 irregularsampling_kernel.c irregularsampling_mexinterface.c -output irregularsampling_mexinterface

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

   This module is a mex interface for MATLAB associated to the irregular
   resampling scheme coming from the research papers

   [1] R. Abergel, S. Ladjal, F. Tupin and J.-M. Nicolas, ``A Complex Spectrum
       Based SAR Image Resampling Method with Restricted Target Sidelobes and
       Statistics Preservation'', in proceedings of the 2017 IEEE International
       Geoscience and Remote Sensing Symposium (IGARSS 2017).

   [2] 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 [2] 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
   'irregularsampling_kernel.c'.


   Documentation of the irregularsampling module:
   ==============================================

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

   Please refer to [1,2] for more details and explanations about the irregular
   resampling algorithm.

*/

#include <mex.h>
#include <fftw3.h>

extern void irregularsampling(fftw_complex*,fftw_complex*,fftw_complex*,int,int,int,int,char); // see source file 'irregularsampling_kernel.c'

/* usage displayer */
static void display_usage()
{
  mexPrintf("\nIrregular resampling scheme for pseudo-raw SAR images (MATLAB MEX INTERFACE).\n\n");
  mexPrintf("Usage: [v0,Tx,Ty] = irregularsampling_mexinterface(u0,K,Ntransl,verbose)\n\n");
  mexPrintf("   u0       : (matrix of complex double) input image\n");
  mexPrintf("   K        : (scalar int32) locality parameter (must be positive)\n");
  mexPrintf("   Ntransl  : (scalar int32) number of tested translations (must be positive)\n");
  mexPrintf("   verbose  : (scalar int32) enter any non-zero value to enable verbose mode\n\n");
  mexPrintf("   v0       : (matrix of complex double) output irregularly resampled image\n");
  mexPrintf("   Tx       : (matrix of double) output horizontal local displacement field\n");
  mexPrintf("   Ty       : (matrix of double) output vertical local displacement field\n\n");
  mexPrintf("Description:\n\n");
  mexPrintf("   This function is an implementation of the irregular resampling scheme\n");
  mexPrintf("   proposed in the research paper [1] and also described in [2]\n\n");
  mexPrintf("    [1] R. Abergel, S. Ladjal, F. Tupin and J.-M. Nicolas, ``A Complex Spectrum\n");
  mexPrintf("        Based SAR Image Resampling Method with Restricted Target Sidelobes and\n");
  mexPrintf("        Statistics Preservation'', in proceedings of the 2017 IEEE International\n");
  mexPrintf("        Geoscience and Remote Sensing Symposium (IGARSS 2017).\n\n");
  mexPrintf("    [2] R. Abergel, L. Denis, S. Ladjal and F. Tupin, ``Sidelobes Suppression\n");
  mexPrintf("        and Strong Targets Extraction in Single Look Complex SAR Images with \n");
  mexPrintf("        Several Applications'', preprint HAL-IMT, submitted in July 2017.\n\n");
  mexPrintf("    This irregular resampling algorithm can be used to process pseudo-raw SAR\n");
  mexPrintf("    image acquired in stripmap mode (pseudo-raw means non-apodized and\n");
  mexPrintf("    non-oversampled SAR images, see [1,2] for more details). It consists in\n");
  mexPrintf("    resampling the input (complex-valued) pseudo-raw image (u0) over a \n");
  mexPrintf("    non-regular grid such as the resampled image (v0) does not exhibits \n");
  mexPrintf("    sidelobes in the vicinity of strong target responses.\n\n");
}

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

  double *u0_re, *u0_im, *v0_re, *v0_im, *Tx, *Ty;
  fftw_complex *u0,*v0,*T;
  size_t nx,ny;
  mwSize adr,x,y,K,Ntransl,verbose;

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

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

  // number of outputs //
  if(nlhs > 3) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:irregularsampling_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:irregularsampling_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:irregularsampling_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 (K = scalar 32 bits integer) //
  if( !mxIsScalar(prhs[1]) || !mxIsInt32(prhs[1])) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:irregularsampling_mexinterface:badinputdatatype","the second input (K) must be of type 'scalar int32'.");
  }
  K = (mwSize) mxGetScalar(prhs[1]);
  if(K < 1) {
    display_usage();
    mexErrMsgIdAndTxt("CPSAR_TOOLBOX:irregularsampling_mexinterface:badinputdatatype","the second input (K) must be greater or equal to 1.");
  }

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

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


  /**********************/
  /* memory allocation */
  /**********************/
  u0 = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex));
  v0 = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex));
  plhs[0] = mxCreateDoubleMatrix((mwSize)ny,(mwSize)nx,mxCOMPLEX);
  if(nlhs > 1) T = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex)); else T = NULL;
  if(nlhs >= 2) plhs[1] = mxCreateDoubleMatrix((mwSize)ny,(mwSize)nx,mxREAL);
  if(nlhs == 3) plhs[2] = mxCreateDoubleMatrix((mwSize)ny,(mwSize)nx,mxREAL);



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

  // recopy u0_re and u0_im into fftw_complex* datatype u0 (I hate myself!) //
  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.
      u0[adr][0] = u0_re[x*ny+y];
      u0[adr][1] = u0_im[x*ny+y];
    }

  // run the irregularsampling algorithm //
  if(verbose) mexPrintf("\n");
  irregularsampling(v0,T,u0,nx,ny,Ntransl,K,verbose!=0);

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

  // deal with optional outputs (Tx,Ty) //
  if(nlhs == 3) {
    Tx = mxGetData(plhs[1]);
    Ty = mxGetData(plhs[2]);
    for(adr=y=0;y<ny;y++)
      for(x=0;x<nx;x++,adr++) {
	Tx[x*ny+y] = T[adr][0];
	Ty[x*ny+y] = T[adr][1];
      }
  }
  else if (nlhs == 2) {
    Tx = mxGetData(plhs[1]);
    for(adr=y=0;y<ny;y++)
      for(x=0;x<nx;x++,adr++)
	Tx[x*ny+y] = T[adr][0];

  }

  /***************/
  /* free memory */
  /***************/
  fftw_free(u0);
  fftw_free(v0);
  if(T) fftw_free(T);

}
