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

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

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

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

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

   To simplify interfacing (for instance with Matlab) the source file
   'sardecomp_kernel.c' only contains the core functions of the SPECKLE PLUS
   TARGETS DECOMPOSITION algorithm.

   A 'main' function (including argument parser and TIFF input/output support)
   corresponding to this module can be found into the source file
   'sardecomp.c', it can be compiled using the following command line.

   gcc -O3 sardecomp.c sardecomp_kernel.c tools.c tiffread.c tiffwrite.c -lfftw3 -lm -ltiff -o sardecomp


   Short description of this module:
   =================================

   The module 'sardecomp' is an implementation of 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).

   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 details in [1].


   Target extraction process:
   ==========================

   The targets are extracted according an "a contrario" criterion, which boils
   done to the design of a Number of False Alarms (NFA) function which can be
   used as a target detector that provides a precise control of the false
   detections.

   More precisely, a NFA-score NFA(k,l) is associated to each pixel (k,l) of the
   image, and we can mathematically prove that, in pure speckle noise H0,

   for any epsilon > 0, we have

   E[#{(k,l) in {0,...,M-1}x{0,...,N-1}, such as NFA(k,l) <= epsilon}] <= epsilon

   where E denotes the mathematical esperance, and # the cardinality. This means
   that when thresholding the NFA function using a threshold epsilon, the
   average number of false detections (i.e., the number of detection made in
   pure noise H0) is less than epsilon.

   The parameter epsilon can be set by the user (the usual setting is epsilon =
   1) and it must be understood that this number represents an upper bound of
   the (average) number of allowed false detections.

   The targets are extracted iteratively from the input image u0, the extraction
   process stops when no target with NFA smaller than epsilon can be found in
   the image.

   See more details and mathematical proofs in [1].


   Inputs/outputs format:
   ======================

   Complex-valued images (u0, w0, s0 and d0) are stored as a stack of 2
   TIFF-float images, the first (respectively the second) image of the stack
   corresponds to the real (respectively imaginary) part of the complex signal.

   Float precision for the pixel values is non-mandatory for input u0 (8-bit, 16
   bits formats are supported), however, using float precision is strongly
   recommanded since saturation strongly affects the Shannon interpolate of the
   image.

   The target features (position, amplitude, NFA-score) are stored in plain text
   into the output 'targets' which as type FILE*. Each line of the the file
   'targets' corresponds to one extracted target, it contains the following
   numbers, separated by a blank space:

                       x  y  tx  ty  Ar  Ai  log10nfa

   those seven real numbers are such as:

   + (x,y) represents the integer coordinates of the pixel that contains the
     target center

   + (tx,ty) are two real numbers in [-0.5,0.5), they represents the relative
     position of the target within pixel (x,y). More precisely, they are such as
     the sub-pixellic position of the target center is (x-tx,y-ty).

   + Ar and Ai represent the real and imaginary parts of the target amplitude.

   + log10nfa represents the log10 of the NFA-score associated to the target.


   Explanations about the main input parameters:
   =============================================

   The input 'epsilon' represents the NFA threshold parameter (= maximal allowed
   number of false detections) discussed above. The algorithm stops
   automatically when no target with NFA less than epsilon can be found in the
   image.

   The input 'nmax' can be used to limit the number of extraction. If nmax is
   nonpositive, the number of extraction is not limited, the extraction process
   only stops when no target with NFA less than epsilon can be found in the
   image.

   The input 'nfaorder' can be used to decide how the extraction is done. If
   nfaorder != 0, the pixels of the image are sorted in decreasing intensity
   order, a target is extracted when its NFA is less than epsilon. Otherwise,
   when nfaorder == 0, the algorithm

     (i) looks for the target having highest NFA value
     (ii) if this NFA value is higher than epsilon, the target is extracted, and
     we go back to (i), otherwise, the algorithm stops.

   'K' and 'Ntransl' are integer inputs of the algorithm. 'K' is a locality
   parameter and 'Ntransl' denotes the number of tested displacements, such as
   the set of tested displacements (= candidate values of Tx(k,l) and Ty(k,l)
   for all (k,l)) is

             {-0.5 + id/Ntransl, for id in {0,...,Ntransl-1}}.

   See [1] for more details about those two parameters.

 */

#include <math.h>
#include <fftw3.h>
#include <time.h>

#ifndef MATLAB_MEX_FILE // STANDARD COMPILATION

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MATLAB_DISPLAY 0

#else // MEX COMPILATION (use several tricks for Matlab printf)

#include <mex.h>
#include <fftw3.h>
#define fflush(stdout) mexEvalString("pause(0);") // force Matlab to display message in real time (otherwise printf messages are displayed at the end of the mex execution)
#define MATLAB_DISPLAY 1

#endif

#define ASSERT_ALLOC(cmd) if(NULL == (cmd)) { printf("Not enough memory\n"); exit(EXIT_FAILURE); }
#define MAX(a,b)        (((a)>(b)) ? (a) : (b))

/* internal structure */
typedef struct sarpixel {
  int x,y; // integer coordinates of the pixel
  double real,imag; // real and imaginary parts of the complex pixel
} *Sarpixel;


/**************/
/* prototypes */
/**************/

// internal modules
void compute_square_ratios(double*,double*,double*,double*,int*,int*,fftw_complex**,fftw_complex**,int,int,int,int,int,int);
void sardecomp(fftw_complex*,fftw_complex*,fftw_complex*,FILE*,double*,int,int,double,int,int,int,char,char);
static void extract_target(fftw_complex*,struct sarpixel*,double*,double*,fftw_complex*,fftw_complex*,fftw_complex*,fftw_complex*,int,int,double,double,int,int,char);
static int compare_intensities(const void*,const void*);
static double learn_distribution_parameter(int,int,int,int);

// external modules (see source file 'tools.c')
extern void horizontal_translation(fftw_complex*,fftw_complex*,int,int,double,fftw_plan);
extern void vertical_translation(fftw_complex*,fftw_complex*,int,int,double,fftw_plan);
extern double sinc(double);


/********************/
/* kernel functions */
/********************/

/* compare_intensities: comparison function for qsort, to sort intensities (or amplitudes) in **decreasing** order. */
static int compare_intensities(const void *a, const void *b)
{
  struct sarpixel *sarpix1,*sarpix2;
  double intensity1,intensity2;
  int out;

  sarpix1 = (struct sarpixel *)a;
  sarpix2 = (struct sarpixel *)b;
  intensity1 = (sarpix1->real)*(sarpix1->real)+(sarpix1->imag)*(sarpix1->imag);
  intensity2 = (sarpix2->real)*(sarpix2->real)+(sarpix2->imag)*(sarpix2->imag);

  if(intensity1 < intensity2) out = 1;
  else if (intensity1 > intensity2) out = -1;
  else out = 0;

  return out;
}


/* learn_distribution_parameter: learn the sigma parameter of the NFA formula.

   This module can be used to learn the sigma parameter of the NFA formula. It
   consists to compute the (square) ratios Rx^2 and Ry^2 in a pure speckle
   noise, i.e. using an image 'u0sample' (with width 'nx' and height 'ny')
   following the H0 model described in the research paper [1].

   The ratios 'Rx' and 'Ry' computed from 'u0sample' both follow a Rayleigh
   distribution with parameter sigma (which depends on 'K' and 'Ntransl'), that
   we estimate using a maximum likelihood estimator, leading to

   sigma^2 =  (||Rx||^2 + ||Ry||^2) / (4.*nx*ny)

   Inputs description:
   ===================

   + 'nx' and 'ny': horizontal and vertical dimensions of the sample image
     u0sample

   + 'K 'and 'Ntransl': the locality parameter and the number of tested
     translations (see the main documentation, above in this module).

   Output description:
   ===================

   This module returns the estimated value of sigma^2.

 */
static double learn_distribution_parameter(nx,ny,K,Ntransl)
     int nx,ny,K,Ntransl;
{
  fftw_complex *u0sample,*dft_u0sample,*dft_vt,**u0sample_tx,**u0sample_ty;
  fftw_plan plan_u0sample,*iplan_u0sample_tx,*iplan_u0sample_ty;
  double sig_square,unused_dble,rx_square,ry_square,delta,a,b;
  int x,y,adr,id_t,unused_int;

  // memory allocation //
  ASSERT_ALLOC(dft_u0sample = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex)));
  ASSERT_ALLOC(dft_vt = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex)));
  ASSERT_ALLOC(u0sample = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex)));
  ASSERT_ALLOC(u0sample_tx = (fftw_complex**) fftw_malloc (Ntransl*sizeof(fftw_complex*)));
  ASSERT_ALLOC(u0sample_ty = (fftw_complex**) fftw_malloc (Ntransl*sizeof(fftw_complex*)));
  ASSERT_ALLOC(plan_u0sample = fftw_plan_dft_2d(ny,nx,u0sample,dft_u0sample,FFTW_FORWARD,FFTW_ESTIMATE));
  ASSERT_ALLOC(iplan_u0sample_tx = (fftw_plan*) malloc (Ntransl*sizeof(fftw_plan)));
  ASSERT_ALLOC(iplan_u0sample_ty = (fftw_plan*) malloc (Ntransl*sizeof(fftw_plan)));
  for(id_t=0;id_t<Ntransl;id_t++) {
    ASSERT_ALLOC(u0sample_tx[id_t]  = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex)));
    ASSERT_ALLOC(u0sample_ty[id_t]  = (fftw_complex*) fftw_malloc (nx*ny*sizeof(fftw_complex)));
    ASSERT_ALLOC(iplan_u0sample_tx[id_t] = fftw_plan_dft_2d(ny,nx,dft_vt,u0sample_tx[id_t],FFTW_BACKWARD,FFTW_ESTIMATE));
    ASSERT_ALLOC(iplan_u0sample_ty[id_t] = fftw_plan_dft_2d(ny,nx,dft_vt,u0sample_ty[id_t],FFTW_BACKWARD,FFTW_ESTIMATE));
  }

  // fill u0sample with random values (normal distribution for real and imaginary parts) //
  srand48((long int) time (NULL));
  for(adr=0;adr<nx*ny;adr++) {
    a = drand48(); b = drand48();
    u0sample[adr][0] = sqrt(-2.*log(a))*cos(2.*M_PI*b);
    a = drand48(); b = drand48();
    u0sample[adr][1] = sqrt(-2.*log(a))*cos(2.*M_PI*b);
  }

  // compute horizontal/vertical translations of u0sample //
  fftw_execute(plan_u0sample);
  for(id_t=0;id_t<Ntransl;id_t++) {
    delta = -0.5 + (double)id_t/(double)Ntransl;
    horizontal_translation(dft_vt,dft_u0sample,nx,ny,delta,iplan_u0sample_tx[id_t]);
    vertical_translation(dft_vt,dft_u0sample,nx,ny,delta,iplan_u0sample_ty[id_t]);
  }

  // compute ratiomaps in u0sample and use maximum likelihood estimator of the Rayleigh distribution //
  sig_square = 0;
  for(adr=y=0;y<ny;y++)
    for(x=0;x<nx;x++,adr++) {
      compute_square_ratios(&rx_square,&ry_square,&unused_dble,&unused_dble,&unused_int,&unused_int,u0sample_tx,u0sample_ty,x,y,K,Ntransl,nx,ny);
      sig_square += 0.5*(rx_square + ry_square);
    }
  sig_square /= (2.*(double)nx*(double)ny);

  // free memory //
  fftw_free(u0sample);
  fftw_free(dft_u0sample);
  fftw_destroy_plan(plan_u0sample);
  fftw_free(dft_vt);
  for(id_t=0;id_t<Ntransl;id_t++) {
    fftw_free(u0sample_tx[id_t]);
    fftw_free(u0sample_ty[id_t]);
    fftw_destroy_plan(iplan_u0sample_tx[id_t]);
    fftw_destroy_plan(iplan_u0sample_ty[id_t]);
  }
  fftw_free(u0sample_tx);
  fftw_free(u0sample_ty);
  free(iplan_u0sample_tx);
  free(iplan_u0sample_ty);

  return sig_square;
}


/* extract_target: remove a target from u0.

   Given the sub-pixellic position of a target present in u0, we estimate the
   target complex amplitude using a maximum likelihood estimator, and we remove
   the target from 'u0'.

   Optionally, the target can be added to 's0' (cardinal sine) and/or to 'd0'
   (discrete Dirac).

   Inputs/outputs description:
   ===========================

   + 'x0' and 'y0': horizontal and vertical integer coordinates of the pixel
     containing the target center.

   + 'tx' and 'ty': the horizontal and vertical relative position of the target
     center within pixel (x0,y0), this means that the sub-pixellic position of
     the center of the target is (x0-tx,y0-ty).

   + 'u0': image where the target will be removed

   + 'nx' and 'ny': width and height of 'u0'

   + 'sarpix': (optional) copy of 'u0' with Sarpixel structure (this structure
     allow to sort the pixels according to their intensity), if NULL != sarpix,
     the target will be also removed from sarpix.

   + 'A_re' and 'A_im': two double pointers, *A_re and *A_im are the real and
     imaginary parts of the estimation of the complex amplitude of the target.

   + 's0': if NULL != s0, all signal removed from 'u0' is added to 's0'.

   + 'd0': if NULL != d0, *A_re is added to the real part of 'd0' at the integer
     position (x0,y0), and *A_im is added to the imaginary part of 'd0' at the
     integer location (x0,y0).

   + 'u0_tx' and 'u0_ty' denote the horizontal and vertical translations of u0
     using the sub-pixellic translation steps 'tx' and 'ty'. More precisely,
     noting U0 the Shannon interpolate of u0, we have

     for all (x,y) in {0,...,nx-1} x {0,...,ny-1},

                     u0_tx[y*nx+x] = U0(x-tx,y) and u0_ty[y*nx+x] = U0(x,y-ty)

   + 'dir' should be set to 0 if nx >= ny, and any non-zero value otherwise.

 */
static void extract_target(u0,sarpix,A_re,A_im,s0,d0,u0_tx,u0_ty,x0,y0,tx,ty,nx,ny,dir)
     fftw_complex *u0,*u0_tx,*u0_ty,*s0,*d0;
     struct sarpixel *sarpix;
     double *A_re,*A_im,tx,ty;
     int x0,y0,nx,ny;
     char dir;
{
  double val_re,val_im,val;
  int x,y,adr;

  // first compute A_re and A_im, the real and imaginary parts of the target
  // amplitude using a maximum likelihood estimator
  if(dir) { // let us loop on the horizontal dimension which is the smallest
    val_re = val_im = 0.;
    for(x=0;x<nx;x++) {
      val = sinc(x-x0+tx);
      val_re += u0_ty[y0*nx+x][0] * val;
      val_im += u0_ty[y0*nx+x][1] * val;
    }
  }
  else { // let us loop on the vertical dimension which is the smallest
    val_re = val_im = 0.;
    for(y=0;y<ny;y++) {
      val = sinc(y-y0+ty);
      val_re += u0_tx[y*nx+x0][0] * val;
      val_im += u0_tx[y*nx+x0][1] * val;
    }
  }
  *A_re = val_re;
  *A_im = val_im;

  // now let us remove from u0 (and sarpix if NULL != sarpix) the target with
  // amplitude A_re + i*A_im located at the sub-pixellic position
  // (x0-tx,y0-ty). Besides if NULL!=s0, then the removed target is added to s0.
  for(adr=y=0;y<ny;y++)
    for(x=0;x<nx;x++,adr++) {
      val = sinc(x-x0+tx)*sinc(y-y0+ty);
      u0[adr][0] -= val_re*val;
      u0[adr][1] -= val_im*val;
      if(s0) {
	s0[adr][0] += val_re*val;
	s0[adr][1] += val_im*val;
      }
      if(sarpix) {
	sarpix[adr].x = x;
	sarpix[adr].y = y;
	sarpix[adr].real = u0[adr][0];
	sarpix[adr].imag = u0[adr][1];
      }
    }

  // if NULL != d0, add to d0 a discrete Dirac same amplitude as the extracted
  // target at the (discrete) location (x0,y0)
  if(d0) {
    adr = y0*nx+x0;
    d0[adr][0] += val_re;
    d0[adr][1] += val_im;
  }

}


/* compute_square_ratios: returns the value of Rx^2 and Ry^2 at the specified
   location.

   Inputs description:
   ===================

   + 'nx' and 'ny' denotes the horizontal (width) and vertical (heigh) dimensions of
   the image.

   + 'K 'and 'Ntransl': the locality parameter and the number of tested
     translations (see the main documentation, above in this module).

   + 'x0' and 'y0': the pixel coordinate position where the square ratios Rx^2
     and Ry^2 must be computed.

   + 'u0_tx' (respectively 'u0_ty'): a stack of images, such as, for any integer
     index integer index 'id' in {0,..,Ntransl-1}, u0_tx[id] (respectively
     u0_ty[id]) contains the horizontal (respectively vertical) translation of
     the pseudo-raw image u0 with subpixellic translation step

     t = -0.5+(double)id/(double)Ntransl.

     More precisely, for all integer pixel positions (x,y) in
     {0,...,nx-1}x{0,...,ny-1}, we have

         u0_tx[id][y*nx+x] = U0(x-t,y), and u0_tx[id][y*nx+x] = U0(x,y-t),

     where U0 denotes the Shannon interpolate of the image u0.

     Outputs description:
     ====================

     The outputs of this modules are pointers to scalar values (double or
     int). The explanations below refer to the values toward which points those
     pointers.

     + 'tx_opt' and 'ty_opt': the best candidate translations found at position (x0,y0), i.e.,

       tx_opt = argmin_{tx in T} J(U0(x0-tx+[-K,K],y0),
       ty_opt = argmin_{ty in T} J(U0(x0,y0-tx+[-K,K]),

       where T = {-0.5+(double)id/(double)Ntransl, for id in {0,...,Ntransl-1}}
       and for any discrete complex valued signal v = v^r + i*v^i with size
       2K+1,

          J(v) = TV(v^r) + TV(v^i),
	  noting TV(s) = sum_{-K <= k <= K, k not in {-1,0}} |s(k+1)-s(k)|

     + 'id_tx_opt' and 'id_ty_opt': translation identifiers associated to
       'tx_opt' and 'ty_opt', such as

       tx_opt = -0.5 + (double)id_tx_opt/(double)Ntransl
       ty_opt = -0.5 + (double)id_ty_opt/(double)Ntransl

     + 'Rx_square' and 'Ry_square': the computed square ratios.

 */
void compute_square_ratios(Rx_square,Ry_square,tx_opt,ty_opt,id_tx_opt,id_ty_opt,u0_tx,u0_ty,x0,y0,K,Ntransl,nx,ny)
  double *Rx_square,*Ry_square,*tx_opt,*ty_opt;
  fftw_complex **u0_tx,**u0_ty;
  int *id_tx_opt,*id_ty_opt,x0,y0,K,Ntransl,nx,ny;
{
  int x,y,adr1,adr2,id_t,delta,id_tx,id_ty;
  double tvx,tvy,tvx_re,tvx_im,tvy_re,tvy_im,vx_re,vx_im,vy_re,vy_im,tvx_opt,tvy_opt,val;

  // loop on the translations (compute masked TV for each translation) //
  tvx_opt = tvy_opt = INFINITY;
  id_tx = id_ty = 0; // to avoid gcc warnings.
  for(id_t=0;id_t<Ntransl;id_t++) {
    tvx_re = tvx_im = tvy_re = tvy_im = 0.;
    for(delta=-K;delta<=K;delta++) {
      if((delta != -1) && (delta != 0)) {
	x = x0+delta; y = y0+delta;
	tvx_re += fabs(u0_tx[id_t][y0*nx+(nx+x+1)%nx][0] - u0_tx[id_t][y0*nx+(nx+x)%nx][0]);
	tvx_im += fabs(u0_tx[id_t][y0*nx+(nx+x+1)%nx][1] - u0_tx[id_t][y0*nx+(nx+x)%nx][1]);
	tvy_re += fabs(u0_ty[id_t][((ny+y+1)%ny)*nx+x0][0] - u0_ty[id_t][((ny+y)%ny)*nx+x0][0]);
	tvy_im += fabs(u0_ty[id_t][((ny+y+1)%ny)*nx+x0][1] - u0_ty[id_t][((ny+y)%ny)*nx+x0][1]);
      }
    }
    tvx = tvx_re + tvx_im;
    tvy = tvy_re + tvy_im;
    if(tvx_opt > tvx) { tvx_opt = tvx; id_tx = id_t; }
    if(tvy_opt > tvy) { tvy_opt = tvy; id_ty = id_t; }
  }

  // compute local masked variances (in both directions) of the optimally translated signal //
  vx_re = vx_im = vy_re = vy_im = 0.;
  for(delta=-K;delta<=K;delta++) {
    x = x0+delta; y = y0+delta;
    adr1 = y0*nx + (nx+x)%nx;
    adr2 = ((ny+y)%ny)*nx + x0;
    if(delta !=0) {
      // horizontal axis
      val = u0_tx[id_tx][adr1][0];
      vx_re += val*val;
      val = u0_tx[id_tx][adr1][1];
      vx_im += val*val;
      // vertical axis
      val = u0_ty[id_ty][adr2][0];
      vy_re += val*val;
      val = u0_ty[id_ty][adr2][1];
      vy_im += val*val;
    }
  }

  // output square ratios (Rx^2,Ry^2) and the corresponding optimal translations (tx_opt,ty_opt) //
  adr1 = y0*nx+x0;
  *Rx_square = 2.* (double)K * (pow(u0_tx[id_tx][adr1][0],2.)/vx_re + pow(u0_tx[id_tx][adr1][1],2.)/vx_im);
  *Ry_square = 2.* (double)K * (pow(u0_ty[id_ty][adr1][0],2.)/vy_re + pow(u0_ty[id_ty][adr1][1],2.)/vy_im);
  *id_tx_opt = id_tx;
  *id_ty_opt = id_ty;
  *tx_opt = -0.5 + (double)(id_tx)/(double)(Ntransl);
  *ty_opt = -0.5 + (double)(id_ty)/(double)(Ntransl);
}


/* sardecomp: speckle plus targets decomposition.

   This module correspond to Algorithm 4 of the research paper [1], it remove
   from the complex pseudo raw image 'u0' all the targets with number of false
   alarms smaller than 'epsilon' (or after 'nmax' extractions if nmax > 0).

   At the end of the process, the image 'u0' corresponds to the speckle
   component 'w0' of the decomposition. The file 'targets' contains the features
   of each extracted target.

   Optionally, the information removed from 'u0' can be saved in 's0' (linear
   combination of sampled cardinal sines), and the linear recombination of the
   targets into discrete Diracs can be saved in 'd0'.

   Inputs description:
   ===================

   + 'epsilon': NFA threshold, this number represents the maximal allowed number
     of false detections.

   + 'verbose': set verbose = 0 to avoid printf, and verbose != 0 otherwise.

   + 'nfaorder': set nfaorder = 0 to extract the targets in decreasing intensity
     order, or set nfaorder != 0 to extract the targets in decreasing NFA
     order. In both cases, only the targets with NFA smaller than 'epsilon' are
     extracted. See more details above in the main documentation of this module.

   + 'nmax': set nmax <= 0 to NOT impose a maximal number of extraction (the
     process stops when no target with NFA smaller than epsilon can be found in
     'u0'), otherwise, at most 'nmax' extractions can be done in 'u0' (the NFA
     criterion is still used to stop the process).

   + 'nx' and 'ny': width and height of 'u0'

   + 'K 'and 'Ntransl': the locality parameter and the number of tested
     translations (see the main documentation, above in this module).

   + 's' s > 0, then set the sigma parameter of the NFA model equal to
     *s, otherwise sigma is learnt in pure speckle noise using a maximum
     likelihood approach.

   Outputs description:
   ====================

   + 'u0': since the detected targets are removed from 'u0', at the end of the
     extraction process, this image represents the speckle component of the
     decomposition 'w0'.

   + 'targets': a FILE in which are saved in plain text the features of each
     extracted target (see above the main documentation of the module for the
     description of the features).

   + 's0': if allocated (s0 != NULL), all information removed from 'u0' is added
     to 's0', so that 's0' is a linear combination of sampled cardinal sine
     functions (one cardinal sine per target extraction), and satisfies

                     input 'u0' = output 'w0' + output 's0'

   + 'd0': if allocated (d0 != NULL), after each extraction of a target with
     sub-pixellic position (x,y) and complex amplitude A, we add the complex
     value 'A' to 'd0' at the integer position (round(x),round(y)). Therefore,
     'd0' represents the linear recombination of the targets into discrete
     Diracs, it is such as the imare R(u0) defined by

                              R(u0) = w0 + d0

     is free of sidelobes effects.

   + 's': a pointer that points to (*s) = the sigma value used in the NFA model.

   Of course, both images 's0' and 'd0' can be computed only given the target
   features (the reciprocal is false), this is why those two outputs are
   optional.

*/
void sardecomp(u0,s0,d0,targets,s,K,Ntransl,epsilon,nmax,nx,ny,nfaorder,verbose)
     fftw_complex *u0,*s0,*d0;
     int nx,ny,Ntransl,K,nmax;
     double *s,epsilon;
     char nfaorder,verbose;
     FILE *targets;
{
  double sig_square,Rmax_square,R_square,Rx_square,Ry_square,R_square_opt,delta,tx_opt,ty_opt,tx_opt2,ty_opt2,A_re,A_im,log10nfa;
  fftw_complex *dft_u0,*dft_vt,**u0_tx,**u0_ty;
  fftw_plan plan_u0,*iplan_u0_tx,*iplan_u0_ty;
  int x0,y0,x0_opt,y0_opt,adr,id_t,id_tx_opt,id_ty_opt,id_tx_opt2,id_ty_opt2,n,Omega=nx*ny;
  const char no_nflag=(nmax <= 0);
  const double coef = log10(2.*(double)Omega);
  const double l10 = log(10.);
  struct sarpixel *sarpix = NULL;
  char dir,found;

  // memory allocation //
  ASSERT_ALLOC(dft_u0 = (fftw_complex*) fftw_malloc (Omega*sizeof(fftw_complex)));
  ASSERT_ALLOC(dft_vt = (fftw_complex*) fftw_malloc (Omega*sizeof(fftw_complex)));
  ASSERT_ALLOC(u0_tx = (fftw_complex**) fftw_malloc (Ntransl*sizeof(fftw_complex*)));
  ASSERT_ALLOC(u0_ty = (fftw_complex**) fftw_malloc (Ntransl*sizeof(fftw_complex*)));
  ASSERT_ALLOC(plan_u0 = fftw_plan_dft_2d(ny,nx,u0,dft_u0,FFTW_FORWARD,FFTW_ESTIMATE));
  ASSERT_ALLOC(iplan_u0_tx = (fftw_plan*) malloc (Ntransl*sizeof(fftw_plan)));
  ASSERT_ALLOC(iplan_u0_ty = (fftw_plan*) malloc (Ntransl*sizeof(fftw_plan)));
  for(id_t=0;id_t<Ntransl;id_t++) {
    ASSERT_ALLOC(u0_tx[id_t]  = (fftw_complex*) fftw_malloc (Omega*sizeof(fftw_complex)));
    ASSERT_ALLOC(u0_ty[id_t]  = (fftw_complex*) fftw_malloc (Omega*sizeof(fftw_complex)));
    ASSERT_ALLOC(iplan_u0_tx[id_t] = fftw_plan_dft_2d(ny,nx,dft_vt,u0_tx[id_t],FFTW_BACKWARD,FFTW_ESTIMATE));
    ASSERT_ALLOC(iplan_u0_ty[id_t] = fftw_plan_dft_2d(ny,nx,dft_vt,u0_ty[id_t],FFTW_BACKWARD,FFTW_ESTIMATE));
  }
  if(!nfaorder) ASSERT_ALLOC(sarpix = (struct sarpixel *) malloc (Omega*sizeof(struct sarpixel)));

  // initialization //
  for(adr=y0=0;y0<ny;y0++)
    for(x0=0;x0<nx;x0++,adr++) {
      if(sarpix) {
	sarpix[adr].x = x0;
	sarpix[adr].y = y0;
	sarpix[adr].real = u0[adr][0];
	sarpix[adr].imag = u0[adr][1];
      }
      if(s0) s0[adr][0] = s0[adr][1] = 0.;
      if(d0) d0[adr][0] = d0[adr][1] = 0.;
    }
  dir = (nx < ny);
  if(verbose && (*s<=0)) {
    printf("Learn parameter 'sigma' of the NFA model\n");
    printf("========================================\n");
    fflush(stdout);
  }
  sig_square = ((*s<=0) ? learn_distribution_parameter(512,512,K,Ntransl) : (*s)*(*s));
  if(*s <= 0) {
    *s = sqrt(sig_square);
    if(verbose) {
      printf("computed sigma = %.17g\n\n",*s);
      fflush(stdout);
    }
  }

  Rmax_square = 2.*sig_square*log(2.*(double)Omega/epsilon);
  A_re = A_im = tx_opt = ty_opt = tx_opt2 = ty_opt2 = Rx_square = Ry_square = 0.; // to avoid gcc warnings
  x0_opt = y0_opt = x0 = y0 = id_tx_opt2 = id_ty_opt2 = 0; // to avoid gcc warnings
  n=0;

  //---------------- MAIN LOOP ----------------//
  if(verbose) {
    printf("Start the speckle plus targets decomposition\n");
    printf("============================================\n");
    fflush(stdout);
  }
  if(nfaorder) { // extract targets in decreasing nfa order (heavy computation time)
    do {

      // precompute horizontal & vertical translations //
      fftw_execute(plan_u0);
      for(id_t=0;id_t<Ntransl;id_t++) {
	delta = -0.5 + (double)id_t/(double)Ntransl;
	horizontal_translation(dft_vt,dft_u0,nx,ny,delta,iplan_u0_tx[id_t]);
	vertical_translation(dft_vt,dft_u0,nx,ny,delta,iplan_u0_ty[id_t]);
      }

      // find position (discrete + sub-pixellic) of the target having highest square ratio R^2 = max(Rx,Ry)^2 //
      R_square_opt = 0;
      for(adr=y0=0;y0<ny;y0++)
	for(x0=0;x0<nx;x0++,adr++) {
	  compute_square_ratios(&Rx_square,&Ry_square,&tx_opt,&ty_opt,&id_tx_opt,&id_ty_opt,u0_tx,u0_ty,x0,y0,K,Ntransl,nx,ny);
	  R_square = MAX(Rx_square,Ry_square); // R-score for pixel (x0,y0)
	  if(R_square > R_square_opt) {
	    R_square_opt = R_square;
	    if(Rx_square < Rmax_square) { id_tx_opt2 = Ntransl/2; tx_opt2 = -0.5+(double)id_tx_opt2/(double)Ntransl; }
	    else { id_tx_opt2 = id_tx_opt; tx_opt2 = tx_opt; }
	    if(Ry_square < Rmax_square) { id_ty_opt2 = Ntransl/2; ty_opt2 = -0.5+(double)id_ty_opt2/(double)Ntransl; }
	    else { id_ty_opt2 = id_ty_opt; ty_opt2 = ty_opt; }
	    x0_opt = x0;
	    y0_opt = y0;
	  }
	}

      // check wether the target is meaningful or not (NFA criterion) //
      found = (R_square_opt >= Rmax_square);

      // if the target candidate is meaningful, extract it from u0 //
      if(found) {
	extract_target(u0,NULL,&A_re,&A_im,s0,d0,u0_tx[id_tx_opt2],u0_ty[id_ty_opt2],x0_opt,y0_opt,tx_opt2,ty_opt2,nx,ny,dir);
	log10nfa = coef - R_square_opt/(2.*sig_square*l10);
	if(targets) fprintf(targets,"%d %d %.10g %.10g %.10g %.10g %.10e\n",x0_opt,y0_opt,tx_opt2,ty_opt2,A_re,A_im,log10nfa);
	if(verbose) { printf("target #%d has been extracted: position=(x=%g,y=%g), complex amplitude A=%g%s%g*i, log10(NFA-score)=%g\n",n,x0_opt+tx_opt,y0_opt+ty_opt,A_re,A_im>0?"+":"-",fabs(A_im),log10nfa); fflush(stdout); }
      }
      n++;

    } while ((found) && (no_nflag || (n < nmax)));
  }
  else { // extract meaningful targets in decreasing intensity order (save computation time)
    do {

      // precompute horizontal & vertical translations //
      fftw_execute(plan_u0);
      for(id_t=0;id_t<Ntransl;id_t++) {
	delta = -0.5 + (double)id_t/(double)Ntransl;
	horizontal_translation(dft_vt,dft_u0,nx,ny,delta,iplan_u0_tx[id_t]);
	vertical_translation(dft_vt,dft_u0,nx,ny,delta,iplan_u0_ty[id_t]);
      }

      // sort sarpixels in decreasing intensity order //
      qsort((void*)sarpix,Omega,sizeof(struct sarpixel),compare_intensities);

      // find a target candidate (= that having the highest intensity among meaningful targets) //
      found = 0; adr = 0;
      while((!found) && adr < Omega) {
	x0 = sarpix[adr].x;
	y0 = sarpix[adr].y;
	compute_square_ratios(&Rx_square,&Ry_square,&tx_opt,&ty_opt,&id_tx_opt,&id_ty_opt,u0_tx,u0_ty,x0,y0,K,Ntransl,nx,ny);
	R_square = MAX(Rx_square,Ry_square); // R-score for pixel (x0,y0)
	found = (R_square >= Rmax_square);
	adr++;
      }


      // if the target candidate is meaningful, extract it from u0 (and sarpix) //
      if(found) {
	if(Rx_square < Rmax_square) { id_tx_opt = Ntransl/2; tx_opt = -0.5+(double)id_tx_opt/(double)Ntransl; } // the target is not meaningfull along the X-axis (prefer zero translation)
	if(Ry_square < Rmax_square) { id_ty_opt = Ntransl/2; ty_opt = -0.5+(double)id_ty_opt/(double)Ntransl; } // the target is not meaningfull along the Y-axis (prefer zero translation)
	extract_target(u0,sarpix,&A_re,&A_im,s0,d0,u0_tx[id_tx_opt],u0_ty[id_ty_opt],x0,y0,tx_opt,ty_opt,nx,ny,dir);
	log10nfa = coef - R_square/(2.*sig_square*l10);
	if(targets) fprintf(targets,"%d %d %.10g %.10g %.10g %.10g %.10e\n",x0,y0,tx_opt,ty_opt,A_re,A_im,log10nfa);
	if(verbose) { printf("target #%d has been extracted: position=(x=%g,y=%g), complex amplitude A=%g%s%g*i, log10(NFA-score)=%g\n",n,x0+tx_opt,y0+ty_opt,A_re,A_im>0?"+":"-",fabs(A_im),log10nfa); fflush(stdout); }
      }
      n++;

    } while ((found) && (no_nflag || (n < nmax)));
  }
  if(verbose) {
    if(!found) printf("All targets with number of false alarms smaller than %g have been extracted.\n\n",epsilon);
    else printf("The process stopped after %d extractions, notice that it may remain some\ntargets with number of false alarms smaller than %g in the processed image.\n\n",nmax,epsilon);
    fflush(stdout);
  }

  // free memory //
  fftw_free(dft_u0);
  fftw_destroy_plan(plan_u0);
  fftw_free(dft_vt);
  for(id_t=0;id_t<Ntransl;id_t++) {
    fftw_free(u0_tx[id_t]);
    fftw_free(u0_ty[id_t]);
    fftw_destroy_plan(iplan_u0_tx[id_t]);
    fftw_destroy_plan(iplan_u0_ty[id_t]);
  }
  fftw_free(u0_tx);
  fftw_free(u0_ty);
  free(iplan_u0_tx);
  free(iplan_u0_ty);
  if(!nfaorder) free(sarpix);

}
