prolate.c

image processing including fourier,wavelet,segmentation etc.

/*--------------------------- MegaWave2 module  -----------------------------*/
/* mwcommand
 name = {prolate};
 author = {"Lionel Moisan"};
 version = {"1.1"};
 function = {"Create a prolate kernel"};
 usage = {
   'n':[n=256]->m  "image size (power of 2, default 256)",
   s->s            "size of kernel in spatial domain (odd, e.g. 3)",
   d->d            "relative diameter of kernel in Fourier domain (0<d<=1)",
   out<-out        "resulting kernel (sxs Fimage)"
};
*/
/*--------------------------------------------------------------------------*/

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

extern void fft2d();

Fimage prolate(s,d,m,out)
     int s;
     float d;
     Fimage out;
     int *m;
{
  Fimage re,im,tmp;
  int n,r,i,x,y,xx,yy,cont;
  float v,rad2;
  double norm,energy,last_energy;

  n = *m;
  rad2 = (float)(n*n)*d*d*0.25;

  tmp = mw_change_fimage(NULL,n,n);
  re = mw_new_fimage();
  im = mw_new_fimage();
  out = mw_change_fimage(out,s,s);
  mw_clear_fimage(out,1./(float)s);

  /* MAIN LOOP */
  i = 0;
  do {
    
    /* copy kernel */
    mw_clear_fimage(tmp,0.);
    for (x=s;x--;) for (y=s;y--;)
      tmp->gray[((y-s/2+n)%n)*n+(x-s/2+n)%n] = out->gray[y*s+x];
    
    /* take Fourier Transform */
    fft2d(tmp,NULL,re,im,0);
    
    /* project on Fourier essential support */
    norm = 0.;
    for (x=n;x--;) for (y=n;y--;) {
      xx = (x>n/2?n-x:x);
      yy = (y>n/2?n-y:y);
      if ((float)(xx*xx+yy*yy)>rad2) re->gray[y*n+x] = im->gray[y*n+x] = 0.;
      else norm += (double)re->gray[y*n+x]*(double)re->gray[y*n+x]
	     +(double)im->gray[y*n+x]*(double)im->gray[y*n+x];
    }
    energy = sqrt(norm)/(double)n;
    
    /* take inverse Fourier Transform */
    fft2d(re,im,tmp,NULL,1);
    
    /* copy new kernel and normalize */
    norm = 0.;
    for (x=s;x--;) for (y=s;y--;) {
      v = tmp->gray[((y-s/2+n)%n)*n+(x-s/2+n)%n];
      out->gray[y*s+x] = v;
      norm += (double)v*(double)v;
    }
    norm = sqrt(norm);
    for (x=s*s;x--;) out->gray[x] /= (float)norm;

    cont = (i==0 || (i<=100 && energy>last_energy*1.00001));
    last_energy = energy;
    i++;
  } while (cont);

  printf("%d iterations, energy = %f\n",i,(float)energy);

  /* normalize total weight */
  norm = 0.;
  for (x=s*s;x--;) {
    v = out->gray[x];
    if (v>0) norm += (double)v; else  norm -= (double)v;
  }
  for (x=s*s;x--;) out->gray[x] /= (float)norm;

  mw_delete_fimage(im);
  mw_delete_fimage(re);
  mw_delete_fimage(tmp);

  return(out);
}