/*

   This program is freely available on the web page

   http://www.mi.parisdescartes.fr/~moisan/p+s

   I hope that you will find it useful.
   If you use it for a publication, please mention 
   this web page and the paper it makes reference to.
   If you modify this program, please indicate in the
   code that you did so and leave this message.
   You can also report bugs or suggestions to 
   lionel.moisan [AT] parisdescartes.fr

*/
/*---------------------------- MegaWave2 module ----------------------------*/
/* mwcommand
name = {perdecomp};
version = {"1.1"};
author = {"Lionel Moisan"};
function = {"Periodic/Smooth image decomposition: u=p+s"};
usage = {
  u->u              "input Fimage",
  p<-p              "periodic component (Fimage)",
  {
    s<-s            "smooth component (optional Fimage)"
  }
};
*/
/*----------------------------------------------------------------------
 v1.0 (04/2007): initial version
 v1.1 (10/2012): simplified version 
----------------------------------------------------------------------*/

#include <stdio.h>
#include <math.h>
#include "mw.h"

Fimage perdecomp(u,p,s)
     Fimage u,p,s;
{
  Fimage re,im,s_orig;
  int nx,ny,adr,x,y;
  double a,cx,cy;
  float b;

  /* memory allocations */
  nx = u->ncol;
  ny = u->nrow;
  s_orig = s;
  s = mw_change_fimage(s,ny,nx);
  p = mw_change_fimage(p,ny,nx);
  re = mw_new_fimage();
  im = mw_new_fimage();
  if (!s || !p || !re || !im) mwerror(FATAL,1,"Not enough memory\n");  

  /* fill "boundary image" */
  mw_clear_fimage(s,0.);
  for (x=0;x<nx;x++) {
    b = u->gray[x]-u->gray[(ny-1)*nx+x];
    s->gray[x] += b;
    s->gray[(ny-1)*nx+x] += -b;
  }
  for (y=0;y<ny;y++) {
    b = u->gray[y*nx]-u->gray[y*nx+nx-1];
    s->gray[y*nx] += b;
    s->gray[y*nx+nx-1] += -b;
  }

  /* Fourier transform */
  fft2d(s,NULL,re,im,0);
  
  /* (0,0) frequency */
  re->gray[0] = im->gray[0] = 0.;

  cx = 2.*M_PI/(double)nx;
  cy = 2.*M_PI/(double)ny;

  /* other frequencies */
  y = 1; /*to avoid (0,0) frequency */
  for (x=0;x<nx;x++) {
    a = cos(cx*(double)x);
    for (;y<ny;y++) {
      adr = y*nx+x;
      b = 0.5/(float)(2.-a-cos(cy*(double)y));
      re->gray[adr] *= b;
      im->gray[adr] *= b;
    }
    y=0;
  }

   /* inverse Fourier transform */
  fft2d(re,im,s,NULL,1);
  
  /* get p=u-s */
  for (adr=nx*ny;adr--;)
    p->gray[adr] = u->gray[adr]-s->gray[adr];

  mw_delete_fimage(im);
  mw_delete_fimage(re);
  if (!s_orig) mw_delete_fimage(s);

  return(p);
}