Logo Search packages:      
Sourcecode: paul version File versions  Download package

fft.c

/*****************************************************************************
 *                                                                           *
 * Program:   paul                                                           *
 *            (P)rogramm zur (A)uswertung und (U)mformung von                *
 *            (L)aserbildern                                                 *
 * Uses:      fftw (Fastest Fourier Transform in the West from MIT)          *
 * Modul:     fft.c                                                          *
 *            fast fourier transform                                         *
 * Author:    Andreas Tille                                                  *
 * Date:      04.10.1998                                                     *
 * Copyright: Andreas Tille, 1999; GNU Public License                        *
 *                                                                           *
 *****************************************************************************/

#ifdef HAVE_DATABOX 
#include <gtkdatabox.h>
#endif

#include "paul.h" 

#ifdef HAVE_LIBFFTW

/*
#define FFT_ABS(c, i, j) log(1.0 + (c).re * (c).re + (c).im * (c).im)
*/

#define FFT_ABS(c)            sqrt((c).re * (c).re + (c).im * (c).im)
#define RADIUS(i, j)          sqrt((i)*(i) + (j)*(j))

#include <math.h>
#include <srfftw.h>

#ifdef HAVE_DATABOX 
static int ShowFFTProfile(GList *piclist);
#endif

int FastFourierTransform(PAUL *p)
/* calculate fast fourier transforms of images
 * It converts image data of a gray (i.e. green for paul) image into real
 * there will be no check, whether the picture is real gray, 
 * the green component is simply used
 * --- Parameter: ---
 * PAUL *p                      : list of images
 * --- Return: ---
 * int   FastFourierTransform() : RET_OK or RET_ERR
 */
{
  static    rfftwnd_plan  plan = NULL;
  PICTURE                *bild;
  GList                  *pl;
  fftw_real              *fft, *fftprof;
  register int            storepix;
  register unsigned char *ap, *aip, *bp, *bip;
  unsigned char          *fip, *fipp;
  register fftw_real     *fftp, *fftpp, *fftip, *fftipp;
  fftw_complex           *fftop, *fftopp;
  fftw_real               max_val, radius, abs;
  int                     ow, i, j, n_prof, *prof;
  char                   *desc;
   
  g_return_val_if_fail ( IS_PAUL(p), RET_ERR ) ;
  g_return_val_if_fail ( CheckPicList(p->piclist), RET_ERR );

  for ( bild = BILD(pl = p->piclist); pl; bild = BILD(pl = pl->next) ) {
    storepix = bild->storepix;
    ow  = ((bild->W)>>1) + 1;
    fft = g_new0(fftw_real, bild->H*(ow<<1));

    for (fip = ( aip = ap = bild->DATA + (storepix == 3 ? 1 : 0) ) + storepix * bild->size, fftp = fft; 
       ap < fip; ap += storepix*bild->W, fftp += (ow<<1))
      for (fipp = ap + storepix * bild->W, fftpp = fftp; aip < fipp; aip += storepix, fftpp++)
        *fftpp = *aip;

    if ( plan && (plan->n[0] != bild->H || plan->n[1] != ow) ) {
      rfftwnd_destroy_plan(plan);
      plan = NULL;
    }
    if ( !plan ) 
      plan  = rfftw2d_create_plan(bild->H, bild->W, 
                                  FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_IN_PLACE);
      /* it is pure intention that bild->H and bild->W used in this order because the  *
       * FFTW people call the vertical distance width and the horizontal distance      *
       * height :). Anyway it works this way as test with simple shapes will show you. */

    rfftwnd_one_real_to_complex(plan, fft, NULL);

    n_prof = MIN(bild->W, bild->H);
    FREE(bild->fftprof);
    bild->fftprof = g_new0(fftw_real, n_prof);
    prof          = g_new0(int, n_prof);
    fftprof       = bild->fftprof;
    
    /* rearange complex field to real and find maximum value for scaling purposes: */
    max_val = 0.0;
    for ( fftip = (fftp = fft) + bild->H*ow, fftop = (fftw_complex *) fft, i = 0; 
          fftp < fftip; fftp += ow, fftop += ow, i++ )
      for (fftipp = (fftpp = fftp) + ow, fftopp = fftop, j = 0; fftpp < fftipp; fftpp++, fftopp++, j++) {
      register int ii = ((i <= bild->H/2) ? i : bild->H - i), ind;
        ind = (int)((radius = RADIUS(ii, j)) + 0.5);
        abs = FFT_ABS(*fftopp);
        if ( ind < n_prof ) {
          prof[ind]++;
        fftprof[ind] = fftprof[ind] + abs;
      }
      if ( (*fftpp = abs * radius) > max_val) max_val = *fftpp;
      }

    for ( i = 0; i < n_prof; i++ ){
      if ( !prof[i] ) continue;
      fftprof[i] = fftprof[i] / (fftw_real)(prof[i]);
    }
    FREE(prof);

    /* write to output image, scaled to the range 0..255 (grayscale): */
    max_val = 1.0 / max_val;

    memset(bild->DATA, 0, storepix*bild->size);  /* especially required for color images */
    /* Move left top of output matrix to bottom right of the image */
    for (ap = bild->DATA + storepix*(bild->W*((bild->H-1)/2) + (bild->W-1)/2) + (storepix == 3 ? 1 : 0), 
         fip = bild->DATA + storepix*bild->size, fftp = fft; ap < fip; ap += storepix*bild->W)
      for (fipp = (aip = ap) + storepix * ow; aip < fipp; aip += storepix, fftp++)
      *aip = (unsigned char)(*fftp * max_val * 255);
    /* Move left bottom of output matrix to top right of the image */
    for (ap = bild->DATA + storepix*((bild->W-1)/2) + (storepix == 3 ? 1 : 0), 
         fip = bild->DATA + storepix*bild->W*((bild->H-1)/2); ap < fip; ap += storepix * bild->W)
      for (fipp = (aip = ap) + storepix * ow; aip < fipp; aip += storepix, fftp++)
      *aip = (unsigned char)(*fftp * max_val * 255);
  /* rotate at middle 180 degrees from right to left */
    for (ap = bild->DATA + storepix * ((bild->W+1)/2) + (storepix == 3 ? 1 : 0), 
         bp = ap - 2*storepix + storepix * bild->W * (bild->H-(bild->H % 2 ? 1 : 2)); 
         bp > bild->DATA; ap += storepix*bild->W, bp -= storepix*bild->W)
      for (fipp = (aip = ap) + storepix * ((bild->W-1)/2), bip = bp; aip < fipp; 
           aip += storepix, bip -=storepix)
      *bip = *aip;

    FREE(fft);

    desc = g_strdup_printf("%s of %s", TypFFT, ImgFileName(bild));
    ImgChunksUpdate(bild, TypFFT, desc, APPFFT, FFT);
    FREE(desc);
  }
#ifdef HAVE_DATABOX 
  if ( p->filelist ) ShowFFTProfile(p->piclist);
#endif
  
  p->opt->f  &= ~FFT;
   
  return RET_OK;
}

#ifdef HAVE_DATABOX 
static int ShowFFTProfile(GList *piclist)
/* display radial FFT profile of all images in piclist
 * --- Parameter: ---
 * GList *piclist          : list with image data
 * --- Return: ---
 * int    ShowFFTProfile() : RET_ERR or RET_OK
 */
{
  static GtkWidget *window = NULL;
  GtkWidget        *curve;
  char             *buf;
  int               col, colstep, i, numsets = 0, maxlen = 0, len, index = 0;
  gfloat           *X, *Y;
  PICTURE          *bild;
  GList            *pl;
  GdkColor          color;
  
  g_return_val_if_fail ( BILD(piclist), RET_ERR );

  for ( bild = BILD(pl = piclist); pl; bild = BILD(pl = pl->next) ) {
    if ( bild->fftprof ) numsets++;
    if ( maxlen < (len = MIN(bild->W, bild->H) - 1) ) maxlen = len;
  }

  if ( !numsets ) {
     g_warning(_("No FFT profiles available"));
     return RET_ERR;
  }
  
  if ( window ) {
    if (!GTK_WIDGET_VISIBLE(window) ) gtk_widget_show_all(window);
    else {
      CurveWidgetDestroy(window, NULL);
      window = NULL;
    }
    return RET_OK;
  }

  if ( NBILDER(piclist) == 1 ) 
    buf = g_strdup_printf(_("FFT profile of %s"), BILD(piclist)->file);
  else
    buf = g_strdup(_("FFT profiles"));

  curve = GTK_WIDGET(gtk_databox_new());
  gtk_signal_connect(GTK_OBJECT(curve), "destroy", GTK_SIGNAL_FUNC(gtk_databox_data_destroy_all),
                     NULL);
  CreateDataboxFrame(curve, buf);
  FREE(buf);
  
  X = g_new0(gfloat, maxlen);

  for (i = 0; i < maxlen; i++ ) X[i] = i;
  
  colstep = (col = 0xFFFF) / numsets;

  for ( bild = BILD(pl = piclist); pl; bild = BILD(pl = pl->next) ) {
    if ( !bild->fftprof ) continue;
    Y   = g_new0(gfloat, maxlen);
    len = MIN(bild->W, bild->H) - 1;
    for (i = 0; i < len; i++ ) 
      Y[i] = log(1 + bild->fftprof[i]);
    color.red = color.green = color.blue = col;
    if ( ! index ) 
      index = gtk_databox_data_add_x_y(GTK_DATABOX(curve), len, X, Y, color, GTK_DATABOX_POINTS, 0);
    else
      gtk_databox_data_add_y(GTK_DATABOX(curve), len, Y, index, color, GTK_DATABOX_POINTS, 0);
    col -= colstep;
  }

  gtk_databox_rescale(GTK_DATABOX(curve));
  gtk_widget_show(curve);

  return RET_OK;
}
#endif /* Do not show FFT profile if libgtkdatabox does not exist */
 
#else

int FastFourierTransform(PAUL *p)
{
  g_warning(_("Sorry. %s was compiled without FFT support.\nInstall FFTW from ftp://theory.lcs.mit.edu/pub/fftw/\n and recompile %s to make it work."), exename, exename);
  return RET_OK;
}

#endif


Generated by  Doxygen 1.6.0   Back to index