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

   This file is a part of the

             COMPLEX PSEUDO-RAW SAR IMAGE PROCESSING TOOLBOX

   Copyright (C) 2017 Rémy Abergel (Remy.Abergel AT parisdescartes.fr), Loïc
   Denis (Loic.Denis AT univ-st-etienne.fr), Saïd Ladjal (Said.Ladjal AT
   telecom-paristech.fr) and Florence Tupin (Florence.Tupin AT
   telecom-paristech.fr).

   This program is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful, but
   WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program.  If not, see <http://www.gnu.org/licenses/>.

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

/*

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

   gcc -O3 fftspectrum.c tiffread.c tiffwrite.c -lfftw3 -lm -ltiff -o fftspectrum


   Description:
   ============

   This module can be used to compute the image

                 out = log(1+|DFT(in)|)   or   out = |DFT(in)|

   noting 'in' the input image and 'DFT' the Discrete Fourier Transform
   operator.

   The output image can be saved in TIFF 8-bits (default) of in TIFF-float
   (option -f). When saved in 8-bits, a contrast-change is applied to the output
   to set its min value is 0 and its max value to 255.

   Note: the frequency coordinates are shifted so that the (0,0) frequency is
   placed at the center of the output image.

 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "fftw3.h"

#define ASSERT_ALLOC(cmd) if(NULL == (cmd)) { printf("Not enough memory\n"); exit(EXIT_FAILURE); }

extern int load_2pages_tiff_into_fftw_complex(fftw_complex**,int*,int*,char*,char); // see source file 'tiffread.c'
extern int save_monopage_striped_tiff_from_double(char*,double*,int,int,char,int); // see source file 'tiffwrite.c'

static void display_usage()
{
  printf("\nCompute DFT spectrum in log scale: out = log(1+|DFT(in)|).\n\n");
  printf("Usage: fftspectrum [-f] in out\n\n");
  printf("   -a  : avoid log scale, save out = |DFT(in)|\n");
  printf("   -c  : set output min/max values equal to 0 and 255 (affine contrast change)\n");
  printf("   -f  : save output as a tiff-float instead of unsigned char tiff\n");
  printf("   in  : input complex image, (stack of 2 TIFF images, first image=real part\n");
  printf("         of 'in', second image=imaginary part of 'in')\n");
  printf("   out : output image (TIFF)\n\n");
}


int main(int argc, char **argv)
{
  char *filenamein = NULL,*filenameout = NULL,aflag,vflag,cflag,fflag;
  double *v = NULL;
  fftw_complex *u, *dft_u = NULL;
  fftw_plan plan_u = NULL;
  double min,max,val,c1,c2;
  int x,y,a,b,adr,i,nx,ny;
  aflag = vflag = cflag = fflag = 0;

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

  /* main parser */
  for(i=1;i<argc;i++) {
    if (strcmp(argv[i],"-c") == 0) cflag = 1;
    else if (strcmp(argv[i],"-f") == 0) fflag = 1;
    else if (strcmp(argv[i],"-a") == 0) aflag = 1;
    else if (NULL == filenamein) filenamein = argv[i];
    else if (NULL == filenameout) filenameout = argv[i];
    else { display_usage(); printf("Error: could not retrieve properly input arguments.\n"); return EXIT_FAILURE; }
  }

  /* check consistency */
  if ((NULL == filenamein) || (NULL == filenameout)) {
    display_usage();
    if (NULL == filenamein) printf("Error: input 'in' is missing\n\n");
    else if(NULL == filenameout) printf("Error: output 'out' is missing\n\n");
    return EXIT_FAILURE;
  }

  /*******************/
  /* Open TIFF Image */
  /*******************/

  if(EXIT_FAILURE == load_2pages_tiff_into_fftw_complex(&u,&nx,&ny,filenamein,vflag)) {
    printf("Error: could not open input '%s'\n",filenamein);
   return EXIT_FAILURE;
  }

  /*********************/
  /* Memory allocation */
  /*********************/
  ASSERT_ALLOC(v = (double*) malloc (nx*ny * sizeof(double)));
  ASSERT_ALLOC(dft_u = (fftw_complex*) fftw_malloc (nx*ny * sizeof(fftw_complex)));
  ASSERT_ALLOC(plan_u = fftw_plan_dft_2d(ny,nx,u,dft_u,FFTW_FORWARD,FFTW_ESTIMATE));

  /*********************************************/
  /* KERNEL: compute DFT spectrum in log scale */
  /*********************************************/
  fftw_execute(plan_u); // set dft_u = DFT(u);
  min = 0.; max = -1;
  // set v = log(1+|DFT(u)|)
  for(y=0,b=-(ny/2);b<=(ny-1)/2;b++,y++)
    for(x=0,a=-(nx/2);a<=(nx-1)/2;a++,x++) {
      adr = ((b+ny)%ny)*nx + (a+nx)%nx; // location of frequency (a,b) in dft_u
      val = aflag ? hypot(dft_u[adr][0],dft_u[adr][1]) : log(1.+hypot(dft_u[adr][0],dft_u[adr][1]));
      v[y*nx+x] = val;
      if(min > val) min = val;
      if(max < val) max = val;
    }
  // apply contrast change (if needed)
  if(cflag && (min < max))
    for(c1=255./(max-min),c2=-c1*min,adr=0;adr<nx*ny;adr++)
      v[adr] = c1*v[adr]+c2;

  /*********************************************/
  /*               END OF KERNEL               */
  /*********************************************/

  /********************************/
  /* Write output as a TIFF Image */
  /********************************/
  if(EXIT_FAILURE == save_monopage_striped_tiff_from_double(filenameout,v,nx,ny,vflag,fflag?0:8)) {
    printf("Error: could not save output '%s'\n",filenameout);
    fftw_free(u);
    fftw_free(dft_u);
    fftw_destroy_plan(plan_u);
    free(v);
    return EXIT_FAILURE;
  }

  /* free memory */
  fftw_free(u);
  fftw_free(dft_u);
  fftw_destroy_plan(plan_u);
  free(v);

  return EXIT_SUCCESS;
}
