#+TITLE: CPSAR-TOOLS: Complex Pseudo-RAW SAR Image Processing Toolbox
#+author: Rémy Abergel, Télécom ParisTech, Université Paris Saclay, Paris, France
#+STARTUP: overview hideblocks

* Note

  We kindly advise you to open this .org file using Emacs and benefit of the
  power of *Orgmode*, making for instance possible to *show/hide* all sections
  by placing the cursor on the corresponding line (begins with character `*`)
  and pressing the <TAB> key on your keyboard. The same functionality holds for
  list items (starting with `+`) or code blocks (starting with
  "#+begin_src"). You will also enjoy other surprises (such as code coloration,
  bold and italic fonts, etc.), and even much more (see [[http://orgmode.org/]]).

* License

  CPSAR-TOOLS : 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 software provides a C implementation (+ Matlab interface) of the
  algorithms described in 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).

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

* Installation
** Content and dependencies

   This software was designed with *as less dependencies as possible*. The
   following libraries are required :

   + fftw :: software + installation instructions available at http://www.fftw.org/
   + libtiff :: For Debian GNU/Linux Operating Systems, libtiff can be installed
                by running the following command line (displayed between
                the #+begin_src and #+end_src tags) from the terminal:

		#+begin_src sh
                  sudo apt-get install libtiff5-dev
		#+end_src

   The remaining dependencies are the following standard libraries :

   + stdio
   + stdlib
   + math
   + string
   + time

   The Input/Output image type is TIFF, you may be able to visualize those
   images using your favourite image displayer. For instance, images can be
   displayed using the software ImageJ which can be installed under GNU/Linux
   distributions using =apt-get=, with the following command :

   #+begin_src sh
    sudo apt-get install imagej
   #+end_src

   Of course the installation of ImageJ is NOT MANDATORY to use this software,
   it is only suggested as a tool for displaying TIFF images.

   Beware that all input/output complex images are represented as a Stack of two
   TIFF images in float precision (first image of the stack = real part, second
   image of the stack = imaginary part). See the Input/Output section of this
   documentation for more details and simple tips to manipulate those images.

** Compilation of the C modules with gcc

   All the C modules of this software can be installed at once by executing the
   'CPSAR_SETUP_C_MODULES' executable file (bash) contained in the 'setup'
   directory of this software. For instance, from the root directory of this
   software, run the following command :

   #+begin_src sh
     ./setup/CPSAR_SETUP_C_MODULES
   #+end_src

   or equivalently, compile each module using gcc (the compilation commands can
   be adapted to other C-compilers). To that aim, place yourself into the 'src'
   directory of this software, and run the following commands:

   #+begin_src sh
     gcc -w -O3 irregularsampling.c irregularsampling_kernel.c tools.c tiffread.c tiffwrite.c -lfftw3 -lm -ltiff -o irregularsampling
     gcc -w -O3 sardecomp.c sardecomp_kernel.c tools.c tiffread.c tiffwrite.c -lfftw3 -lm -ltiff -o sardecomp
     gcc -w -O3 tiffprintasc.c tiffread.c -ltiff -o tiffprintasc
     gcc -w -O3 tiffreadasc.c tiffwrite.c -ltiff -o tiffreadasc
     gcc -w -O3 tiffmerge.c tiffread.c tiffwrite.c -ltiff -o tiffmerge
     gcc -w -O3 tiffsplit.c tiffread.c tiffwrite.c -ltiff -o tiffsplit
     gcc -w -O3 fftspectrum.c tiffread.c tiffwrite.c -lfftw3 -lm -ltiff -o fftspectrum
     gcc -w -O3 modulus.c tiffread.c tiffwrite.c -lfftw3 -lm -ltiff -o modulus
     gcc -w -O3 tiffadd.c tiffread.c tiffwrite.c -lfftw3 -lm -ltiff -o tiffadd
   #+end_src

   When the installation is done, you can find the executable files listed below
   in the 'src' directory of the software:

   + irregularsampling
   + sardecomp
   + tiffprintasc
   + tiffreadasc
   + tiffmerge
   + tiffsplit
   + fftspectrum
   + modulus
   + tiffadd

   You can check the usage and see a (short) documentation of each module by
   running the corresponding executable file without argument. For instance, if
   you place yourself in the 'src' directory, the following instruction

   #+begin_src
     ./irregularsampling
   #+end_src

   should result in the following message in your terminal:

   #+begin_output

     Irregular resampling scheme for pseudo-raw SAR images.

     Usage: irregularsampling [-N N] [-K K] [-T T] [-v] [-vv] u0 v0

        -N N    : (default = 20) number of tested translations
        -K K    : (default = 25) locality parameter
        -T T    : output local displacement field (stack of 2 TIFF Float images)
        -v      : light verbose mode (display algorithm progression)
        -vv     : heavy verbose mode (also display information about TIFF I/O)
        u0      : input image (a stack of 2 TIFF images, see description below)
        v0      : output image (a stack of 2 TIFF images, see description below)

     Description of the 2-pages TIFF I/O format :

        Input (u0) is a complex-valued pseudo-raw images, it must be stored as a
        stack of two images into a single TIFF file, such as the first (respectively
        second) image of the stack represents the real (resp. imaginary) part of u0.

        Output (v0) is the irregularly resampled image produced by this algorithm,
        it is stored as a stack of two TIFF-float images using the same convention.

        Optional output (T) is a stack of two TIFF-Float images, the first (resp.
        second) image of the stack contains the horizontal (resp. vertical)
        component of the displacement field.

     Error: input 'u0' is missing

   #+end_output

** Compilation of the MATLAB interfaces using mex

   To compile the Matlab interfaces provided in CPSAR-TOOLS, open a Matlab, add
   to your Matlab search path the root directory of this software and run the
   'CPSAR_SETUP_MATLAB_MODULES.m' (the file is located in the 'setup' directory
   of the toolbox).

   #+begin_src matlab
     cd PATH_TO_THE_ROOT_DIRECTORY_OF_THIS_TOOLBOX % change this line
     addpath(genpath('.')) % add CPSAR-TOOLS (and subdirectory) to your current search path
     run CPSAR_SETUP_MATLAB_MODULES.m % compilation of the mex interfaces
   #+end_src

   or equivalently, compile each module using mex. To that aim, place yourself
   into the 'src' directory of this software, and run the following commands in
   your Matlab terminal:

   #+begin_src matlab
     mex -silent -lm -lfftw3 CFLAGS="\$CFLAGS -std=c99" tools.c irregularsampling_kernel.c irregularsampling_mexinterface.c -output irregularsampling_mexinterface
     mex -silent -lm -lfftw3 CFLAGS="\$CFLAGS -std=c99" tools.c sardecomp_kernel.c sardecomp_mexinterface.c -output sardecomp_mexinterface
   #+end_src

   When the installation is done, the functions 'irregularsampling_mexinterface'
   and 'sardecomp_mexinterface' are available from your Matlab terminal. To
   check their corresponding usage and see a (short) documentation, run the
   function without argument. More precisely, the commands

   #+begin_src matlab
     irregularsampling_mexinterface;
     sardecomp_mexinterface;
   #+end_src

   should respectively result in the following messages in your Matlab terminal :

   #+begin_output

     Irregular resampling scheme for pseudo-raw SAR images (MATLAB MEX INTERFACE).

     Usage: [v0,Tx,Ty] = irregularsampling_mexinterface(u0,K,Ntransl,verbose)

        u0       : (matrix of complex double) input image
        K        : (scalar int32) locality parameter (must be positive)
        Ntransl  : (scalar int32) number of tested translations (must be positive)
        verbose  : (scalar int32) enter any non-zero value to enable verbose mode

        v0       : (matrix of complex double) output irregularly resampled image
        Tx       : (matrix of double) output horizontal local displacement field
        Ty       : (matrix of double) output vertical local displacement field

     Description:

        This function is an implementation of the irregular resampling scheme
        proposed in the research paper [1] and also described in [2]

         [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, ``Sidelobes Suppression
             and Strong Targets Extraction in Single Look Complex SAR Images with
             Several Applications'', preprint HAL-IMT, submitted in July 2017.

         This irregular resampling algorithm can be used to process pseudo-raw SAR
         image acquired in stripmap mode (pseudo-raw means non-apodized and
         non-oversampled SAR images, see [1,2] for more details). It consists in
         resampling the input (complex-valued) pseudo-raw image (u0) over a
         non-regular grid such as the resampled image (v0) does not exhibits
         sidelobes in the vicinity of strong target responses.

     Error using irregularsampling_mexinterface
     Incorrect number of inputs

   #+end_output

   and

   #+begin_output

     Speckle plus targets decomposition (MATLAB MEX INTERFACE).

     Usage: [w0,s0,d0,sigma] = sardecomp_mexinterface(u0,epsilon,K,Ntransl,nmax,nfaorder,verbose,filename,s)

        u0       : (matrix of complex double) input complex pseudo-raw image
        epsilon  : (scalar double) maximal allowed number of false detections
        K        : (scalar int32) locality parameter (must be positive)
        Ntransl  : (scalar int32) number of tested translations (must be positive)
        nmax     : (scalar int32) maximal number of extraction (set nmax=-1 for no limit)
        nfaorder : (scalar int32) set to nonzero to extract target in decreasing nfa order
        verbose  : (scalar int32) enter any non-zero value to enable verbose mode
        filename : (string) path to output file containing the target features

        s        : sigma parameter of the NFA model (set s<=0 to learn this parameter)
        w0       : (matrix of complex double) output speckle component
        s0       : (matrix of complex double) cardinal sines recombination s0=u0-w0
        d0       : (matrix of double) recombination of the targets into discrete Diracs
        sigma    : sigma parameter used in the NFA model (sigma = s if s > 0)

     Inputs/Outputs description:

        This function 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, ``Sidelobes Suppression and
            Strong Targets Extraction in Single Look Complex SAR Images with Several
            Applications'', preprint HAL-IMT, submitted in July 2017.

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

     Error using sardecomp_mexinterface
     Incorrect number of inputs

   #+end_output


   _IMPORTANT NOTE_

   once the mex compilation is done, you can use the MATLAB functions
   'irregularsampling' and 'sardecomp' (see both .m files of the 'src' directory
   of this software), which is better fitted to MATLAB standards than the
   *_mexinterface described above.

   In particular, the 'irregularsampling' and 'sardecomp' Matlab functions
   provide the handling of optional input arguments (traditional 'Name-Value'
   pair of arguments) and proper datatype conversion before calling the mex
   interface.

   You can check their corresponding usage and read a (short) documentation by
   using the Matlab keyword 'help' followed by the name of the function, or
   equivalently by running the function without argument. More precisely, the
   commands

   #+begin_src matlab
     irregularsampling; % or equivalently, type: help irregularsampling
     sardecomp; % or equivalently, type: help sardecomp
   #+end_src

   yields the following messages in your Matlab terminal:

   #+begin_output

     Usage: [v0,Tx,Ty] = irregularsampling(u0,Name,Value)

     Input(s)/Output(s):

       u0       : (input complex matrix) complex pseudo-raw image
       v0       : (output complex matrix) irregularly resampled image
       Tx       : (output matrix of double) horizontal component of the local
                  displacement field
       Ty       : (output matrix of double) vertical component of the local
                  displacement field

     Optional Name-Value pair arguments:

       'K' K        : (scalar, default K=25) locality parameter
       'Ntransl' N  : (scalar, default N=20) number of tested translations
       'verbose' v  : (scalar logical, default v=false) set v=true to display
                       progression in the Matlab console

     Description: irregular resampling scheme for pseudo-raw SAR images.

       This function is an implementation of the irregular resampling scheme
       proposed in the research paper [1] and also described in [2]

        [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, ``Sidelobes Suppression
            and Strong Targets Extraction in Single Look Complex SAR Images with
            Several Applications', preprint HAL-IMT, submitted in July 2017.

        This irregular resampling algorithm can be used to process pseudo-raw SAR
        image acquired in stripmap mode (pseudo-raw means non-apodized and
        non-oversampled SAR images, see [1,2] for more details). It consists in
        resampling the input (complex-valued) pseudo-raw image (u0) over a
        non-regular grid such as the resampled image (v0) does not exhibits
        sidelobes in the vicinity of strong target responses.

     Error using irregularsampling (line 86)
     Incorrect number of input(s)
   #+end_output

   and

   #+begin_output

     Usage: [w0,s0,d0,targets,sigma] = sardecomp(u0,Name,Value)

     Input(s)/Output(s):

       u0      : (input complex matrix) complex pseudo-raw image
       w0      : (output complex matrix) speckle component of u0
       s0      : (output complex matrix) recombination of the extracted targets
                 s0=u0-w0 (linear combination of cardinal sines)
       d0      : (output complex matrix) recombination of the extracted targets
                 into discrete Diracs
       targets : (output matrix of double) target(i,:) contains the 7 features
                 of the target #(i-1), we have

                            target(i,:) = [x,y,tx,ty,Ar,Ai,log10nfa],
                 where,

                 + (x,y)   : integer position of the pixel that contains the
                             target center (top left corner of the image
                             correspond to coordinate (0,0)

                 + (tx,ty) : relative position of the target center in the pixel
                             (x,y), the pixel that contains the target center,
                             which means that the sub-pixellic position of the
                             target is (x-tx,y-ty)

                 + (Ar,Ai) : real and imaginary parts of the complex amplitude
                             A=Ar+1i*Ai of the target

                 + log10nfa: log10 of the NFA-score of the target.
       sigma : (scalar double) return the value of the sigma parameter used in
               the NFA model (sigma = s when optional input s is set to a
               positive value)

     Optional Name-Value pair arguments:

       'epsilon' eps : (scalar, eps=1.) maximal allowed number of false
                     detection (or number of false alarms, NFA), the target
                     extraction process stops when no target with NFA
                     smaller than eps can be found in the image
       'K' K        : (scalar, default K=25) locality parameter
       'Ntransl' N  : (scalar, default N=20) number of tested translations
       'sigma' s    : (scalar, default s=-1) set a positive value for s
                       to set sigma = s in the NFA model (otherwise the sigma
                       parameter of the NFA model is learnt)
       'nmax' n     : (scalar) if used, the extraction process will stop
                      after at most n extractions
       'nfaorder' v : (scalar logical, default v=false) extract targets
                      in decreasing NFA order (this operation is time
                      consuming)
       'verbose' v  : (scalar logical, default v=false) set v=true to display
                      progression in the Matlab console

     Description: Speckle plus targets decomposition algorithm for pseudo-raw
     SAR images.

       This function 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, ``Sidelobes Suppression and
           Strong Targets Extraction in Single Look Complex SAR Images with Several
           Applications', preprint HAL-IMT, submitted in July 2017.

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

   #+end_output

* Software description

  This software is composed of several source files dedicated to

  + Irregular resampling scheme for complex pseudo-raw SAR images ::

    This is an implementation of the irregular resampling scheme proposed in the
    research paper [1] (Algorithm 1), and also described in [2] (see Algorithm
    1).

    This irregular resampling algorithm can be used to process pseudo-raw SAR
    images acquired in stripmap mode (pseudo-raw means non-apodized and
    non-oversampled SAR images, see [1,2] for more details). It consists in
    resampling the input (complex-valued) pseudo-raw image over a non-regular
    grid such as the resampled image does not exhibit sidelobes in the vicinity
    of strong target responses. It was observed in [1] and [2] that this scheme
    preserves the spatial resolution of the input pseudo-raw image, as well as
    the signal statistics in fully developped speckle area.

  + Speckle + targets decomposition of complex pseudo-raw SAR images ::

    This is is an implementation of the 'SPECKLE PLUS TARGETS DECOMPOSITION'
    algorithm coming from the research paper [2] (see Algorithm 4 therein for a
    pseudo-code description)

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

  + Simple tools for the manipulation of complex-valued image in TIFF format ::

    including pointwise addition between two TIFF images, evaluation of the
    modulus of a complex-valued image, image saturation (to visualize images
    with large dynamic), discrete Fourier transform evaluation, etc. Some of
    these modules are inspired from the MegaWave image processing library [3].

  We give below a short description of those modules.

** Irregular resampling algorithm

   + irregularsampling_kernel.c ::

     This source file contains all the core functions of the irregular
     resampling algorithm (does not contain any main function). This module uses
     some routines available in the 'tools.c' source file (those routines are
     shared by both the 'irregular resampling' algorithm and the 'speckle +
     targets decomposition' algorithm)

   + irregularsampling.c ::

     This source files provides a main function (including argument parser) to
     use the 'irregular resampling' algorithm coded in the source file
     'irregularsampling_kernel.c'

   + irregularsampling_mexinterface.c ::

     This source file provides a MEX interface to use in MATLAB the irregular
     resampling algorithm coded in the source file 'irregularsampling_kernel.c'

   + irregularsampling.m ::

     This source file is a MATLAB function that calls the
     'irregularsampling_mexinterface' but provides the handling of optional
     input arguments (traditional 'Name-Value' pair of arguments) and proper
     datatype conversion before calling the mex interface.

** Speckle + targets decomposition algorithm

   + sardecomp_kernel.c ::

     This source file contains all the core functions of the speckle + targets
     decomposition algorithm (does not contain any main function). This module
     uses some routines available in the 'tools.c' source file (those routines
     are shared by both the 'irregular resampling' algorithm and the 'speckle +
     targets decomposition' algorithm)

   + sardecomp.c ::

     This source files provides a main function (including argument parser) to
     use the 'speckle + targets decomposition' algorithm coded in the source file
     'sardecomp_kernel.c'

   + sardecomp_mexinterface.c ::

     This source file provides a MEX interface to use in MATLAB the 'speckle +
     targets decomposition' algorithm coded in the source file
     'sardecomp_kernel.c'

   + sardecomp.m ::

     This source file is a MATLAB function that calls the
     'sardecomp_mexinterface' but provides the handling of optional input
     arguments (traditional 'Name-Value' pair of arguments) and proper datatype
     conversion before calling the mex interface.

** Other modules

  + tools.c :: A source file shared by both algorithms ('irregular resampling'
               and 'speckle + targets decomposition'), it contains the
               sub-pixellic translation routines and the cardinal sine function.
  + fftspectrum.c :: compute the DFT spectrum of a complex-valued TIFF image
  + modulus.c :: Compute the modulus of a complex-valued TIFF image (including
                 saturation options)
  + tiffadd.c :: Compute the pointwise addition of two complex-valued TIFF
                 images
  + tiffmerge.c :: Merge two (monopage) TIFF images into a stack of two TIFF
                   images (float precision).
  + tiffsplit.c :: Split a stack of two TIFF images into two monopage TIFF
                   images (float precision).
  + tiffprintasc.c :: print the graylevels of a TIFF image (row major order) in
                      ascii format. This module is useful to convert TIFF image
                      in another format, for instance into a MegaWave Fimage (by
                      a joint use of this module and the MegaWave function
                      freadasc).
  + tiffreadasc.c :: read an image in ascii format (in row-major order). This
                     module is useful to convert in TIFF an image given with
                     another image format, for instance to convert a MegaWave
                     Fimage into TIFF or TIFF-Float (by the joint use of this
                     module and the MegaWave function fprintasc).
  + tiffread.c :: routines dedicated to TIFF image reading. This module does not
                  contain any main function.
  + tiffwrite.c :: routines dedicated to TIFF image writting. This module does
                   not contain any main function.

* Input/output TIFF image format

  Complex-valued input/output images are stored as a stack of two images into a
  single TIFF file, such as the first (respectively the second) image of the
  stack represents the real (respectively the imaginary) part of the complex
  image.

  Output images are saved in float precision (unless explicitly stated, see the
  function usage in you terminal by typing the function name without argument).

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

  Several tools included to this toolbox can be used to efficiently manipulate
  this image datatype:

  you can use 'modulus' to compute the modulus of the complex-valued image and
  save it as a standard real-valued TIFF image (8-bits of float precision,
  saturation options are also available).

  You can also use 'tiffsplit' to split a Stack of two TIFF images into two
  single TIFF-float images.

  And you can use 'tiffmerge' to merge two single TIFF images into a Stack of
  two TIFF-float images.

* Examples
** Important note about the example SAR image embedded to this software

   In order to respect image copyrights, we do not provide real SAR images in
   this software. However, a *synthetic* pseudo-raw SAR image is provided in the
   'data' directory in two different formats, yielding the two files:

   + '/data/u0.tif': a multipage TIFF image (first page = real part, second page
     = imaginary part)

   + '/data/u0.mat': same image saved as a MAT-file (MATLAB formatted binary
     file).

   The image has been synthesized from a stack of 26 registered looks of
   complex-valued TerraSAR-X images (*those original TerraSAR-X images are not
   provided in this software*) according to the following procedure:

   1. We resampled all images of the stack to the Shannon-Nyquist critical
      sampling rate (or equivalently, we removed the zero-padding from the
      Fourier spectrum) and we inverted the Hamming apodization from each image
      to retrieve the corresponding pseudo-raw.

   2. We applied our speckle plus target decomposition algorithm to each image
      of the stack. Thus, each pseudo-raw image (u0) of the stack was decomposed
      into a speckle component (w0) and a target component (s0: linear
      combination of cardinal sine functions), such as u0 = w0 + s0;

   3. We computed a 26-look averaging of all the 26 speckle components and this
      (real-valued) averaged image was denoised using the MuLog+BM3D algorithm
      recently proposed in

      Ch.-A. Deledalle, L. Denis, S. Tabti and F. Tupin, ``MuLoG, or How to
      apply Gaussian denoisers to multi-channel SAR speckle reduction?'', IEEE
      Transactions on Image Processing, 2017,

      yielding a reflectivity map denoted 'R'.

   4. We simulated a single look observation of the reflectivity map 'R' by
      mutliplying its gray levels by a complex-valued synthetic speckle, and we
      added to this noisy image the target component 's0' extracted from the
      first image of the pseudo-raw stack.

   5. The resulting image was saved as 'u0' in TIFF and MAT-File formats.

   This way, we computed a realistic (but still synthetic) pseudo-raw SAR image
   that we use below to illustrate the capabilities of the two algorithms
   provided in this software.

** Visualization of a complex pseudo-raw SAR image

   A sample pseudo-raw SAR image is included to this software in the directory
   'data'. The image is complex but you can visualize its amplitude (and use
   saturation to change the dynamic).

   Given a complex-valued image, let us compute the image v = |u0| and use
   saturation to improve display

   Open a terminal, place yourself in the root directory of this software and
   run the following instructions.

   #+begin_src sh

     # open the complex image using imagej (the image is a Stack of two TIFF floats, 1st image = real part, 2nd image = imaginary part) #
     imagej ./data/u0.tif

     # display usage of the 'modulus' routine #
     ./src/modulus

     # compute the modulus of u0 (see details about several options below) #
     # option '-q q' is used to saturate a ratio of 'q' pixels with maximal values
     # option '-c' is used to to normalize the output pixel values in [0,255]
     # default output datatype is TIFF-8bit (so that it can be read by most of the standard image displayer softwares)
     ./src/modulus -c -q 0.5 ./data/u0.tif /tmp/abs_u0.tif
     imagej /tmp/abs_u0.tif

     # you can also use the option '-f' to save the output image as a TIFF-Float image #
     # here let us saturate 0.5% of pixels with maximal values, but we do not use option -c (so that output values can be outside of [0,255])
     ./src/modulus -f -q 0.5 ./data/u0.tif /tmp/abs_u0.tif
     imagej /tmp/abs_u0.tif

     # here we do not saturate at all (due to its high dynamic, the image is very dark)
     ./src/modulus -f ./data/u0.tif /tmp/abs_u0.tif
     imagej /tmp/abs_u0.tif

   #+end_src

** Irregular resampling

   Open a terminal, place yourself in the root directory of this software and
   run the following instructions.

   #+begin_src sh

     # display usage #
     ./src/irregularsampling

     # run the algorithm with the default setting (K=25, N=20) and verbose mode (option '-v') #
     ./src/irregularsampling -v ./data/u0.tif /tmp/v0.tif

     # compute and display the (saturated) modulus of the input image (u0) and the output image (v0)
     ./src/modulus -c -q 0.5 ./data/u0.tif /tmp/abs_u0.tif
     ./src/modulus -c -q 0.5 /tmp/v0.tif /tmp/abs_v0.tif
     imagej -p 1 /tmp/abs_u0.tif /tmp/abs_v0.tif

     # you can also run the algorithm with other parameters (here, we set K=5) #
     ./src/irregularsampling -K 5 -v ./data/u0.tif /tmp/v0.tif
     ./src/modulus -c -q 0.5 ./data/u0.tif /tmp/abs_u0.tif
     ./src/modulus -c -q 0.5 /tmp/v0.tif /tmp/abs_v0.tif
     imagej -p 1 /tmp/abs_u0.tif /tmp/abs_v0.tif

   #+end_src

   We can perform the same experiment in Matlab (do not forget to add the root
   directory of this software to Matlab search path if not permanently added).

   #+begin_src matlab

     %% Add the toolbox to Matlab search path
     path_cpsar = 'PATH_TO_THE_ROOT_DIRECTORY_OF_THIS_TOOLBOX' % you must change this line (enter absolute or relative path to the root directory of this software)
     addpath(genpath(path_cpsar)); % add CPSAR-TOOLS (and subdirectory) to your current search path

     %% display help
     help irregularsampling; % or simply type irregularsampling

     %% load the complex-valued TIFFimage in the Matlab memory
     u0 = imread(sprintf('%s/data/u0.tif',path_cpsar),1) + 1i*imread(sprintf('%s/data/u0.tif',path_cpsar),2);
     % alternatively, you can load the data from the 'u0.mat' file in the 'data' directory
     % load(sprintf('%s/data/u0.mat',path_cpsar));

     %% run the algorithm with the default setting (K=25, N=20) and verbose mode (option '-v')
     v0 = irregularsampling(u0,'verbose',true);
     figure(); cpsarimshow(abs(u0),'Name','input pseudo-raw image u0','saturation',0.5);
     figure(); cpsarimshow(abs(v0),'Name','output irregularly resampled image v0','saturation',0.5);

     %% you can also run the algorithm with other parameters (here, we set K=5)
     v0 = irregularsampling(u0,'K',5,'verbose',true);
     figure(); cpsarimshow(abs(u0),'Name','input pseudo-raw image u0','saturation',0.5);
     figure(); cpsarimshow(abs(v0),'Name','output irregularly resampled image v0','saturation',0.5);

   #+end_src

   ADDITIONAL EXPERIMENT (in Matlab)

   Let us check that the irregular resampling algorithm *preserves the speckle
   statistics*: in this experiments, we check that the statistics of a fully
   developed speckle are preserved by the irregular resampling procedure.

   #+begin_src matlab

     %% Add the toolbox to Matlab search path
     path_cpsar = 'PATH_TO_THE_ROOT_DIRECTORY_OF_THIS_TOOLBOX' % you must change this line (enter absolute or relative path to the root directory of this software)
     addpath(genpath(path_cpsar)); % add CPSAR-TOOLS (and subdirectory) to your current search path

     %% Apply the irregular resampling algorithm to a synthetic image u0 of fully developed speckle
     u0 = randn(512,512) + 1i*randn(512,512); % synthetic image of fully developed speckle
     v0 = irregularsampling(u0,'verbose',true);

     %% compare the histograms of real(u0) and real(v0): both are Gaussian with zero mean an unitary variance
     figure();
     % plot the empirical histogram of real(u0(:)) and real(v0(:)), same experiment can be done with the imaginary part
     histogram(real(u0(:)),100,'Normalization','pdf'); hold on;
     histogram(real(v0(:)),100,'Normalization','pdf');
     % plot the probability density function of a Gaussian N(0,1)
     a = gca(); x0 = a.XLim(1); x1 = a.XLim(2);
     x = linspace(x0,x1,1e6); pdf = (1/sqrt(2*pi)) * exp(-x.^2/2);
     plot(x,pdf);
     legend('empirical p.d.f. of real(u0)','empirical p.d.f. of real(v0)','theoretical p.d.f: law N(0,1)');

     %% Evaluate the spatial correlation (with distance 1 pixel) in u0 and v0
     corr_u0 = sum(sum(u0(1:end-1,:).*conj(u0(2:end,:))))/norm(u0(:))^2; % correlations in u0
     corr_v0 = sum(sum(v0(1:end-1,:).*conj(v0(2:end,:))))/norm(v0(:))^2; % correlations in v0
     fprintf('correlation in u0 : %g%%\n',100*real(corr_u0)); % should be close to 0
     fprintf('correlation in v0 : %g%%\n',100*real(corr_v0)); % we observe a correlation coefficient close to 4% (which is very few in comparaison with to correlations observed after applying Hamming apodization)

   #+end_src

** Speckle + Targets decomposition

   Open a terminal, place yourself in the root directory of this software and
   run the following instructions.

   #+begin_src sh

     # display usage #
     ./src/sardecomp

     # LONG EXPERIMENT (5404 extractions): extract ALL targets with NFA smaller than epsilon (use the automatic stopping criterion, thus, option -n is not used) #
     # The algorithm is run the algorithm with the default setting (K=25, N=20),
     # verbose mode (use option '-v'), and uses a deterministic value of sigma
     # (option '-s') (otherwise sigma is learnt in pure speckle).  Exactly 5404
     # targets will be extracted from u0, we also retrieve optional output images
     # 's0' and 'd0' using options '-s0' and '-d0' ./src/sardecomp -s 1.116 -s0
     # /tmp/s0.tif -d0 /tmp/d0.tif -v ./data/u0.tif /tmp/w0.tif /tmp/targets.txt

     # ALTERNATIVELY (incomplete extraction) you can limit the process to a given number of extractions using option '-n' #
     # e.g., for a limit of 100 extractions, using
     # ./src/sardecomp -n 100 -s0 /tmp/s0.tif -d0 /tmp/d0.tif -v ./data/u0.tif /tmp/w0.tif /tmp/targets.txt

     # have a look to the file containing the extracted targets (each line contains the features of one extracted target) #
     less /tmp/targets.txt

     # have a look to the speckle component 'w0' #
     # remark that after extracting the most strongest bright targets,
     # the image has a small dynamic and saturation is no longer mandatory
     ./src/modulus -f /tmp/w0.tif /tmp/abs_w0.tif
     imagej /tmp/abs_w0.tif

     # have a look to the target component 's0' (linear combination of sampled cardinal sines)
     # and the recombination of the targets into discrete Diracs 'd0'
     ./src/modulus -q 0.1 -f /tmp/s0.tif /tmp/abs_s0.tif
     ./src/modulus -f /tmp/d0.tif /tmp/abs_d0.tif
     imagej -p 1 /tmp/abs_s0.tif /tmp/abs_d0.tif

     # compute R(u0) = w0 + d0 (replace the cardinal sines by discrete Diracs to remove sidelobes) #
     ./src/tiffadd /tmp/w0.tif /tmp/d0.tif /tmp/r0.tif
     ./src/modulus -q 0.1 -f /tmp/r0.tif /tmp/abs_r0.tif
     imagej /tmp/abs_r0.tif

     # also check the additive decomposition: u0 = w0 + s0
     ./src/tiffadd /tmp/w0.tif /tmp/s0.tif /tmp/sum.tif
     ./src/modulus -c -q 0.5 /tmp/sum.tif /tmp/abs_sum.tif
     ./src/modulus -c -q 0.5 ./data/u0.tif /tmp/abs_u0.tif
     imagej -p 1 /tmp/abs_u0.tif /tmp/abs_sum.tif

   #+end_src

   We can perform the same experiment in Matlab (do not forget to add the root
   directory of this software to Matlab search path if not permanently added).

   #+begin_src matlab

     %% Add the toolbox to Matlab search path
     path_cpsar = 'PATH_TO_THE_ROOT_DIRECTORY_OF_THIS_TOOLBOX' % you must change this line (enter absolute or relative path to the root directory of this software)
     addpath(genpath(path_cpsar)); % add CPSAR-TOOLS (and subdirectory) to your current search path

     %% display help
     help sardecomp; % or simply type sardecomp

     %% load the complex-valued TIFFimage in the Matlab memory
     u0 = imread(sprintf('%s/data/u0.tif',path_cpsar),1) + 1i*imread(sprintf('%s/data/u0.tif',path_cpsar),2);

     %% LONG EXPERIMENT (5404 extractions): extract ALL targets with NFA smaller than epsilon #
     % The algorithm is run the algorithm with the default setting (K=25, N=20),
     % verbose mode (use option ['verbose',true]), and uses a deterministic value of
     % sigma (option ['sigma',1.116]) (otherwise sigma is learnt in pure speckle).
     % Exactly 5404 targets will be extracted from u0.
     [w0,s0,d0,targets] = sardecomp(u0,'verbose',true,'sigma',1.116);

     %% ALTERNATIVELY (incomplete extraction) you can limit the process to a given number of extractions using option 'nmax'
     % e.g., for a limit of 100 extractions, using [w0,s0,d0,targets] = sardecomp(u0,'verbose',true,'nmax',100);

     %% have a look to the speckle component 'w0'
     % remark that after extracting the most strongest bright targets, the image has
     % a small dynamic and saturation is no longer mandatory (this is especially true
     % if we do not limit the extraction to 100 targets)
     figure(); cpsarimshow(abs(w0),'Name','speckle component w0');
     figure(); cpsarimshow(abs(w0),'saturation',0.1,'Name','speckle component w0 (0.1%% saturation)');

     %% have a look to the target component 's0' (linear combination of sampled cardinal sines)
     % and the recombination of the targets into discrete Diracs 'd0'
     figure(); cpsarimshow(abs(s0),'saturation',0.1,'Name','target component s0 (0.1%% saturation)');

     %% compute R(u0) = w0 + d0 (replace the cardinal sines by discrete Diracs to remove sidelobes)
     figure(); cpsarimshow(abs(w0+d0),'saturation',0.5,'Name','R(u0) = w0 + d0 (0.5%% saturation)');

     %% also check the additive decomposition: u0 = w0 + s0
     diff = u0-(w0+s0);
     fprintf('average l2 distance between u0 and (w0+s0) = %g\n',sqrt(mean(abs(diff(:)).^2)));

   #+end_src

   ADDITIONAL EXPERIMENT (in Matlab):

   Let us estimate the average number of false detections done by our
   algorithm. To that aim, we run the algorithm on synthetic images of fully
   developed speckle and we count the average number of detections done in those
   images (each detection done in pure speckle noise is a false detection).

   According to our NFA model, the average number of false detections should be
   less than the 'epsilon' input parameter of our algorithm (i.e. setting
   epsilon = 1 should yields an average number of false detections smaller than
   1). In practice, the average number of false detections obtained using our
   algorithm can be slightly above espsilon, but still very close to epsilon. We
   refer to the research paper [2] for explanations on this matter.

   #+begin_src matlab

     %% Add the toolbox to Matlab search path
     path_cpsar = 'PATH_TO_THE_ROOT_DIRECTORY_OF_THIS_TOOLBOX' % you must change this line (enter absolute or relative path to the root directory of this software)
     addpath(genpath(path_cpsar)); % add CPSAR-TOOLS (and subdirectory) to your current search path

     %% run the algorithm on pure speckle images with size 128 x 128 (do it 100 times)
     ny = 128; nx = 128;
     epsilon = 1;
     Nsimu = 100;
     Ndetections = zeros(Nsimu,1);
     for id = 1:Nsimu
       u0 = randn(ny,nx) + 1i*randn(ny,nx);
       [~,~,~,targets] = sardecomp(u0,'epsilon',epsilon);
       Ndetections(id) = size(targets,1);
       fprintf(' + simulation %d/%d : number of false detections = %d\n',id,Nsimu,Ndetections(id));
     end
     fprintf('Average number of false detections = %g, input value of epsilon = %g\n',mean(Ndetections),epsilon);

   #+end_src

* Citations

  Feel free to use this code. However, if you use it for your research or in
  software code, please be so kind as to cite the papers [1] and [2] (see
  section 'Footnotes').

* Versioning

  + version 1.0 (11/2017): first public release.

* Contact

  If you use, adapt or improve this software, we would love to hear about it!
  Feel free to contact us using the following email adresses:

  [[mailto:Remy.Abergel%20AT%20parisdescartes.fr][Rémy Abergel]] : Remy.Abergel AT parisdescartes.fr
  [[mailto:Loic.Denis%20AT%20univ-st-etienne.fr][Loïc Denis]] : Loic.Denis AT univ-st-etienne.fr
  [[mailto:Said.Ladjal%20AT%20telecom-paristech.fr][Saïd Ladjal]] : Said.Ladjal AT telecom-paristech.fr
  [[mailto:Florence.Tupin%20AT%20telecom-paristech.fr][Florence Tupin]] : Florence.Tupin AT telecom-paristech.fr

* Footnotes

[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).

[3] MegaWave, http://megawave.cmla.ens-cachan.fr/.
