Rewrite i.fft, i.ifft

Output maps are FP
Don't use cell_misc/<map>/fft{real,imag} files, use the output maps
Don't rescale to powers of two (not needed with FFTW)
Use fft2() with FFTW-style interleaved array, not fft() with separate real/imaginary arrays
Miscellaneous clean-up


git-svn-id: https://svn.osgeo.org/grass/grass/trunk@34353 15284696-431f-4ddb-bdfa-cd5b030d7da7

imagery/i.fft/Makefile

@@ -2,8 +2,8 @@ MODULE_TOPDIR = ../..
 
 PGM = i.fft
 
-LIBES     = $(IMAGERYLIB) $(GMATHLIB) $(GISLIB) $(FFTWLIB)
-DEPENDENCIES= $(IMAGERYDEP) $(GMATHDEP) $(GISDEP)
+LIBES     = $(GMATHLIB) $(GISLIB) $(FFTWLIB) $(MATHLIB)
+DEPENDENCIES= $(GMATHDEP) $(GISDEP)
 
 include $(MODULE_TOPDIR)/include/Make/Module.make
 

imagery/i.fft/fft_colors.c

@@ -1,24 +0,0 @@
-#include <grass/gis.h>
-#include "globals.h"
-
-int fft_colors(void)
-{
-    struct Colors colors;
-    struct Range range;
-    CELL min, max;
-
-    /* make a real component color table */
-    G_read_range(Cellmap_real, G_mapset(), &range);
-    G_get_range_min_max(&range, &min, &max);
-    G_make_wave_colors(&colors, min, max);
-    G_write_colors(Cellmap_real, G_mapset(), &colors);
-    G_free_colors(&colors);
-
-    /* make a imag component color table */
-    G_read_range(Cellmap_imag, G_mapset(), &range);
-    G_get_range_min_max(&range, &min, &max);
-    G_make_wave_colors(&colors, min, max);
-    G_write_colors(Cellmap_imag, G_mapset(), &colors);
-
-    return 0;
-}

imagery/i.fft/globals.h

@@ -1,6 +0,0 @@
-#ifndef __L_GLOBALS_H__
-#define __L_GLOBALS_H__
-
-extern char Cellmap_real[], Cellmap_imag[];
-
-#endif /* __L_GLOBALS_H__ */

imagery/i.fft/i.fft.html

@@ -8,35 +8,20 @@ components in frequency space.
 
 <H2>NOTES</H2>
 
-The real and imaginary components are stored as arrays of
-doubles in the <EM>cell_misc</EM> directory (for use in the
-inverse transform program,
-
-<EM><A HREF="i.ifft.html">i.ifft</A></EM>),
-
-and are also scaled and formatted into the
+The real and imaginary components are stored into the
 <B>real_image</B> and <B>imaginary_image</B> raster map
-layers for inspection, masking, etc.  In these raster map
+layers.  In these raster map
 layers the low frequency components are in the center and
 the high frequency components are toward the edges.  The
-<B>input_image</B> need not be square;  before
-processing, the X and Y dimensions of the
-<B>input_image</B> are padded with zeroes to the next
-highest power of two in extent (i.e., 256 x 256 is
-processed at that size, but 200 x 400 is padded to 256 x
-512).  The cell category values for viewing, etc., are
-calculated by taking the natural log of the actual values
-then rescaling to 255, or whatever optional range is given
-on the command line, as suggested by Richards (1986).  A
+<B>input_image</B> need not be square. A
 color table is assigned to the resultant map layer.
 
 
 <P>
 
-The current geographic region and mask settings are
-respected when reading the input file.  The presence of a
-mask will, in general, make the resulting fast Fourier
-transform invalid, or at least difficult to interpret.
+The current geographic region and mask settings are respected when
+reading the input file. The presence of nulls or a mask will make the
+resulting fast Fourier transform invalid.
 
 <H2>SEE ALSO</H2>
 

imagery/i.fft/local_proto.h

@@ -1,15 +0,0 @@
-#ifndef __LOCAL_PROTO_H__
-#define __LOCAL_PROTO_H__
-
-#include <grass/gis.h>
-
-/* fft_colors.c */
-int fft_colors(void);
-
-/* orig_wind.c */
-int put_orig_window(struct Cell_head *);
-
-/* save_fft.c */
-int save_fft(int, double *[2], double *, double *);
-
-#endif /* __LOCAL_PROTO_H__ */

imagery/i.fft/main.c

@@ -39,31 +39,37 @@
 #include <grass/gis.h>
 #include <grass/gmath.h>
 #include <grass/glocale.h>
-#include "globals.h"
-#include "local_proto.h"
 
-char Cellmap_real[GNAME_MAX], Cellmap_imag[GNAME_MAX];
+static void fft_colors(const char *name)
+{
+    struct Colors wave, colors;
+    struct FPRange range;
+    DCELL min, max;
+
+    G_read_fp_range(name, G_mapset(), &range);
+    G_get_fp_range_min_max(&range, &min, &max);
+    G_make_wave_colors(&wave, min, max);
+    G_abs_log_colors(&colors, &wave, 100);
+    G_write_colors(name, G_mapset(), &colors);
+    G_free_colors(&colors);
+}
 
 int main(int argc, char *argv[])
 {
     /* Global variable & function declarations */
-    int Range;
-    char Cellmap_orig[GNAME_MAX];
+    struct GModule *module;
+    struct {
+	struct Option *orig, *real, *imag;
+    } opt;
+    const char *Cellmap_real, *Cellmap_imag;
+    const char *Cellmap_orig;
     int inputfd, realfd, imagfd;	/* the input and output file descriptors */
-    char *inmapset;		/* the input mapset name */
     struct Cell_head window;
-    CELL *cell_row, *cell_row2;
-    double max, min, scale, temp;
-
-    int i, j;			/* Loop control variables */
-    int or, oc;			/* Original dimensions of image */
-    int rows, cols;		/* Smallest powers of 2 >= number of rows & columns */
+    DCELL *cell_real, *cell_imag;
+    int rows, cols;		/* number of rows & columns */
     long totsize;		/* Total number of data points */
-    double *data[2];		/* Data structure containing real & complex values of FFT */
-    int save_args();		/* function to stash the command line arguments */
-    struct GModule *module;
-    struct Option *op1, *op2, *op3, *op4;
-    int maskfd;
+    double (*data)[2];		/* Data structure containing real & complex values of FFT */
+    int i, j;			/* Loop control variables */
 
     G_gisinit(argv[0]);
 
@@ -73,199 +79,134 @@ int main(int argc, char *argv[])
 	_("Fast Fourier Transform (FFT) for image processing.");
 
     /* define options */
-    op1 = G_define_option();
-    op1->key = "input_image";
-    op1->type = TYPE_STRING;
-    op1->required = YES;
-    op1->multiple = NO;
-    op1->gisprompt = "old,cell,raster";
-    op1->description = _("Input raster map being fft");
-
-    op2 = G_define_option();
-    op2->key = "real_image";
-    op2->type = TYPE_STRING;
-    op2->required = YES;
-    op2->multiple = NO;
-    op2->gisprompt = "new,cell,raster";
-    op2->description = _("Output real part arrays stored as raster map");
-
-    op3 = G_define_option();
-    op3->key = "imaginary_image";
-    op3->type = TYPE_STRING;
-    op3->required = YES;
-    op3->multiple = NO;
-    op3->gisprompt = "new,cell,raster";
-    op3->description = _("Output imaginary part arrays stored as raster map");
-
-    op4 = G_define_option();
-    op4->key = "range";
-    op4->type = TYPE_INTEGER;
-    op4->required = NO;
-    op4->multiple = NO;
-    op4->answer = "255";
-    op4->description = _("Range of values in output display files");
-
-    /*call parser */
+    opt.orig = G_define_option();
+    opt.orig->key = "input_image";
+    opt.orig->type = TYPE_STRING;
+    opt.orig->required = YES;
+    opt.orig->multiple = NO;
+    opt.orig->gisprompt = "old,cell,raster";
+    opt.orig->description = _("Input raster map being fft");
+
+    opt.real = G_define_option();
+    opt.real->key = "real_image";
+    opt.real->type = TYPE_STRING;
+    opt.real->required = YES;
+    opt.real->multiple = NO;
+    opt.real->gisprompt = "new,cell,raster";
+    opt.real->description = _("Output real part arrays stored as raster map");
+
+    opt.imag = G_define_option();
+    opt.imag->key = "imaginary_image";
+    opt.imag->type = TYPE_STRING;
+    opt.imag->required = YES;
+    opt.imag->multiple = NO;
+    opt.imag->gisprompt = "new,cell,raster";
+    opt.imag->description = _("Output imaginary part arrays stored as raster map");
+
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
-    strcpy(Cellmap_orig, op1->answer);
-    strcpy(Cellmap_real, op2->answer);
-    strcpy(Cellmap_imag, op3->answer);
-
-    /* open input cell map */
-    if ((inmapset = G_find_cell(Cellmap_orig, "")) == NULL)
-	G_fatal_error(_("Raster map <%s> not found"), Cellmap_orig);
+    Cellmap_orig = opt.orig->answer;
+    Cellmap_real = opt.real->answer;
+    Cellmap_imag = opt.imag->answer;
 
-    inputfd = G_open_cell_old(Cellmap_orig, inmapset);
+    inputfd = G_open_cell_old(Cellmap_orig, "");
     if (inputfd < 0)
-	exit(EXIT_FAILURE);
+	G_fatal_error(_("Unable to open input map <%s>"), Cellmap_orig);
 
-    if ((maskfd = G_maskfd()) >= 0)
+    if (G_maskfd() >= 0)
 	G_warning(_("Raster MASK found, consider to remove "
 		    "(see man-page). Will continue..."));
 
-    sscanf(op4->answer, "%d", &Range);
-    if (Range <= 0)
-	G_fatal_error(_("Range less than or equal to zero not allowed."));
-
     G_get_set_window(&window);	/* get the current window for later */
-    put_orig_window(&window);
 
     /* get the rows and columns in the current window */
-    or = G_window_rows();
-    oc = G_window_cols();
-    rows = G_math_max_pow2(or);
-    cols = G_math_max_pow2(oc);
+    rows = G_window_rows();
+    cols = G_window_cols();
     totsize = rows * cols;
 
-    /*  fprintf(stderr,"Power 2 values : %d rows %d columns\n",rows,cols); */
-
     /* Allocate appropriate memory for the structure containing
-       the real and complex components of the FFT.  DATA[0] will
-       contain the real, and DATA[1] the complex component.
+       the real and complex components of the FFT.  data[...][0] will
+       contain the real, and data[...][1] the complex component.
      */
-    data[0] = (double *)G_malloc((rows * cols) * sizeof(double));
-    data[1] = (double *)G_malloc((rows * cols) * sizeof(double));
-
-    /* Initialize real & complex components to zero */
-    G_message(_("Initializing data...\n"));
-    {
-	register double *dptr1, *dptr0;
-
-	dptr0 = data[0];
-	dptr1 = data[1];
-	for (i = 0; i < totsize; i++) {
-	    *dptr0++ = *dptr1++ = 0.0;
-	}
-    }
+    data = G_malloc(rows * cols * 2 * sizeof(double));
 
     /* allocate the space for one row of cell map data */
-    cell_row = G_allocate_cell_buf();
+    cell_real = G_allocate_d_raster_buf();
+    cell_imag = G_allocate_d_raster_buf();
+
+#define C(i, j) ((i) * cols + (j))
 
     /* Read in cell map values */
     G_message(_("Reading the raster map..."));
-    for (i = 0; i < or; i++) {
-	if (G_get_map_row(inputfd, cell_row, i) < 0)
+    for (i = 0; i < rows; i++) {
+	if (G_get_d_raster_row(inputfd, cell_real, i) < 0)
 	    G_fatal_error(_("Error while reading input raster map."));
-	for (j = 0; j < oc; j++)
-	    *(data[0] + (i * cols) + j) = (double)cell_row[j];
+	for (j = 0; j < cols; j++) {
+	    data[C(i, j)][0] = cell_real[j];
+	    data[C(i, j)][1] = 0.0;
+	}
     }
-    /* close input cell map and release the row buffer */
+
+    /* close input cell map */
     G_close_cell(inputfd);
-    G_free(cell_row);
 
     /* perform FFT */
     G_message(_("Starting FFT..."));
-    fft(-1, data, totsize, cols, rows);
+    fft2(-1, data, totsize, cols, rows);
     G_message(_("FFT completed..."));
 
-    /* set up a window for the transform cell map */
-    window.rows = rows;
-    window.cols = cols;
-    window.south = window.north - window.rows * window.ns_res;
-    window.east = window.cols * window.ew_res + window.west;
-    G_set_window(&window);
-
-    /* open the output cell maps and allocate cell row buffers */
-    if ((realfd = G_open_cell_new(Cellmap_real)) < 0)
-	exit(1);
-    if ((imagfd = G_open_cell_new(Cellmap_imag)) < 0)
-	exit(1);
-    cell_row = G_allocate_cell_buf();
-    cell_row2 = G_allocate_cell_buf();
+    /* open the output cell maps */
+    if ((realfd = G_open_fp_cell_new(Cellmap_real)) < 0)
+	G_fatal_error(_("Unable to open real output map <%s>"), Cellmap_real);
+    if ((imagfd = G_open_fp_cell_new(Cellmap_imag)) < 0)
+	G_fatal_error(_("Unable to open imaginary output map <%s>"), Cellmap_imag);
+
+#define SWAP1(a, b)				\
+    do {					\
+	double temp = (a);			\
+	(a) = (b);				\
+	(b) = temp;				\
+    } while (0)
+
+#define SWAP2(a, b)				\
+    do {					\
+	SWAP1(data[(a)][0], data[(b)][0]);	\
+	SWAP1(data[(a)][1], data[(b)][1]);	\
+    } while (0)
 
     /* rotate the data array for standard display */
     G_message(_("Rotating data..."));
+    for (i = 0; i < rows; i++)
+	for (j = 0; j < cols / 2; j++)
+	    SWAP2(C(i, j), C(i, j + cols / 2));
+    for (i = 0; i < rows / 2; i++)
+	for (j = 0; j < cols; j++)
+	    SWAP2(C(i, j), C(i + rows / 2, j));
+
+    G_message(_("Writing transformed data..."));
+
     for (i = 0; i < rows; i++) {
-	for (j = 0; j < cols / 2; j++) {
-	    temp = *(data[0] + i * cols + j);
-	    *(data[0] + i * cols + j) = *(data[0] + i * cols + j + cols / 2);
-	    *(data[0] + i * cols + j + cols / 2) = temp;
-	    temp = *(data[1] + i * cols + j);
-	    *(data[1] + i * cols + j) = *(data[1] + i * cols + j + cols / 2);
-	    *(data[1] + i * cols + j + cols / 2) = temp;
-	}
-    }
-    for (i = 0; i < rows / 2; i++) {
 	for (j = 0; j < cols; j++) {
-	    temp = *(data[0] + i * cols + j);
-	    *(data[0] + i * cols + j) =
-		*(data[0] + (i + rows / 2) * cols + j);
-	    *(data[0] + (i + rows / 2) * cols + j) = temp;
-	    temp = *(data[1] + i * cols + j);
-	    *(data[1] + i * cols + j) =
-		*(data[1] + (i + rows / 2) * cols + j);
-	    *(data[1] + (i + rows / 2) * cols + j) = temp;
+	    cell_real[j] = data[C(i, j)][0];
+	    cell_imag[j] = data[C(i, j)][1];
 	}
+	G_put_d_raster_row(realfd, cell_real);
+	G_put_d_raster_row(imagfd, cell_imag);
     }
 
-    G_message(_("Writing transformed data to file..."));
-    /* write out the double arrays to cell_misc/file/FFTREAL and FFTIMAG */
-    max = 0.0;
-    min = 0.0;
-    save_fft(totsize, data, &max, &min);
-
-    G_message(_("Writing viewable versions of transformed data to files..."));
-    /* Write out result to a new cell map */
-    /*
-       for (i=0; i<rows; i++) {
-       for (j=0; j<cols; j++) {
-       *(cell_row+j) = (CELL) (log(1.0+fabs(*(data[0]+i*cols+j)
-       )) * scale);
-       *(cell_row2+j) = (CELL) (log(1.0+fabs(*(data[1]+i*cols+j
-       ))) * scale);
-       }
-     */
-    scale = (double)Range / log(1.0 + max > -min ? max : -min);
-    {
-	register double *data0, *data1;
-	register CELL *cptr1, *cptr2;
-
-	for (i = 0; i < rows; i++) {
-	    data0 = data[0] + i * cols;
-	    data1 = data[1] + i * cols;
-	    cptr1 = cell_row;
-	    cptr2 = cell_row2;
-	    for (j = 0; j < cols; j++) {
-		*cptr1++ = (CELL) (log(1.0 + fabs(*data0++)) * scale);
-		*cptr2++ = (CELL) (log(1.0 + fabs(*data1++)) * scale);
-	    }
-	    G_put_raster_row(realfd, cell_row, CELL_TYPE);
-	    G_put_raster_row(imagfd, cell_row2, CELL_TYPE);
-	}
-    }
     G_close_cell(realfd);
     G_close_cell(imagfd);
-    G_free(cell_row);
-    G_free(cell_row2);
 
-    /* set up the color tables for histogram streched grey scale */
-    fft_colors();
+    G_free(cell_real);
+    G_free(cell_imag);
+
+    /* set up the color tables */
+    fft_colors(Cellmap_real);
+    fft_colors(Cellmap_imag);
 
     /* Release memory resources */
-    for (i = 0; i < 2; i++)
-	G_free(data[i]);
+    G_free(data);
 
     G_done_msg(_("Transform successful."));
 

imagery/i.fft/orig_wind.c

@@ -1,20 +0,0 @@
-
-#include <grass/gis.h>
-#include <grass/glocale.h>
-#include "globals.h"
-
-#define FFTWINDOW "fftwindow"
-
-
-int put_orig_window(struct Cell_head *hd)
-{
-    char buffer[100];
-
-    /* save the window */
-    sprintf(buffer, "cell_misc/%s", Cellmap_real);
-    G__put_window(hd, buffer, FFTWINDOW);
-    sprintf(buffer, "cell_misc/%s", Cellmap_imag);
-    G__put_window(hd, buffer, FFTWINDOW);
-
-    return 0;
-}

imagery/i.fft/save_fft.c

@@ -1,42 +0,0 @@
-#include <stdio.h>
-#include <grass/gis.h>
-#include <grass/glocale.h>
-#include "globals.h"
-
-
-int save_fft(int total, double *data[2], double *maximum, double *minimum)
-{
-    int i;
-    double max, min, *temp;
-    FILE *fp;
-
-    max = *maximum;
-    min = *minimum;
-
-    if ((fp = G_fopen_new_misc("cell_misc", "fftreal", Cellmap_real)) == NULL)
-	G_fatal_error(_("Unable to open file in the cell_misc directory."));
-    fwrite((char *)data[0], sizeof(double), (size_t) total, fp);
-    fclose(fp);
-
-    if ((fp = G_fopen_new_misc("cell_misc", "fftimag", Cellmap_imag)) == NULL)
-	G_fatal_error(_("Unable to open file in the cell_misc directory."));
-    fwrite((char *)data[1], sizeof(double), (size_t) total, fp);
-    fclose(fp);
-
-    temp = data[0];
-    for (i = 0; i < total; i++, temp++) {
-	max = (max > *temp) ? max : *temp;
-	min = (min < *temp) ? min : *temp;
-    }
-
-    temp = data[1];
-    for (i = 0; i < total; i++, temp++) {
-	max = (max > *temp) ? max : *temp;
-	min = (min < *temp) ? min : *temp;
-    }
-
-    *maximum = max;
-    *minimum = min;
-
-    return 0;
-}

imagery/i.ifft/Makefile

@@ -2,8 +2,8 @@ MODULE_TOPDIR = ../..
 
 PGM = i.ifft
 
-LIBES     = $(IMAGERYLIB) $(GMATHLIB) $(GISLIB) $(FFTWLIB)
-DEPENDENCIES= $(IMAGERYDEP) $(GMATHDEP) $(GISDEP)
+LIBES     = $(GMATHLIB) $(GISLIB) $(FFTWLIB) $(MATHLIB)
+DEPENDENCIES= $(GMATHDEP) $(GISDEP)
 
 include $(MODULE_TOPDIR)/include/Make/Module.make
 

imagery/i.ifft/globals.h

@@ -1,6 +0,0 @@
-#ifndef __L_GLOBALS_H__
-#define __L_GLOBALS_H__
-
-extern char Cellmap_real[], Cellmap_imag[];
-
-#endif /* __L_GLOBALS_H__ */

imagery/i.ifft/i.ifft.html

@@ -20,22 +20,6 @@ When <EM>i.ifft</EM> is executed, it (automatically) uses the same GRASS
 region definition setting that was used during the original transformation
 done with <EM><A HREF="i.fft.html">i.fft</A></EM>.
 
-
-<P>
-
-The real and imaginary components are read from arrays of
-doubles in the <EM>cell_misc</EM> directory (produced by
-the forward transform program,
-
-<EM><A HREF="i.fft.html">i.fft</A></EM>),
-
-and the reconstructed image will preserve the cell value
-scaling of the original image processed by
-
-<EM><A HREF="i.fft.html">i.fft</A></EM>.  No color
-table is assigned to the output map;  one should be created
-before viewing the <EM>output_image</EM>.
-
 <H2>SEE ALSO</H2>
 
 M. Frigo and S. G. Johnson (1998): "FFTW: An Adaptive Software Architecture

imagery/i.ifft/local_proto.h

@@ -1,5 +0,0 @@
-/* orig_wind.c */
-int get_orig_window(struct Cell_head *, char *, char *);
-
-/* save_fft.c */
-int save_fft(int, double *[2], double *, double *);

imagery/i.ifft/main.c

@@ -18,29 +18,38 @@
 #include <grass/gis.h>
 #include <grass/glocale.h>
 #include <grass/gmath.h>
-#include "globals.h"
-#include "local_proto.h"
 
-char Cellmap_real[GNAME_MAX], Cellmap_imag[GNAME_MAX];
+static void fft_colors(const char *name)
+{
+    struct Colors colors;
+    struct FPRange range;
+    DCELL min, max;
+
+    /* make a real component color table */
+    G_read_fp_range(name, G_mapset(), &range);
+    G_get_fp_range_min_max(&range, &min, &max);
+    G_make_grey_scale_fp_colors(&colors, min, max);
+    G_write_colors(name, G_mapset(), &colors);
+}
 
 int main(int argc, char *argv[])
 {
     /* Global variable & function declarations */
-    char Cellmap_orig[GNAME_MAX];
-    FILE *realfp, *imagfp;	/* the input and output file descriptors */
-    int outputfd, maskfd;	/* the input and output file descriptors */
-    char *realmapset, *imagmapset;	/* the input mapset names */
-    struct Cell_head orig_wind, realhead;
-    CELL *cell_row, *maskbuf = NULL;
+    struct GModule *module;
+    struct {
+	struct Option *orig, *real, *imag;
+    } opt;
+    const char *Cellmap_real, *Cellmap_imag;
+    const char *Cellmap_orig;
+    int realfd, imagfd,  outputfd, maskfd;	/* the input and output file descriptors */
+    struct Cell_head realhead, imaghead;
+    DCELL *cell_real, *cell_imag;
+    CELL *maskbuf;
 
     int i, j;			/* Loop control variables */
-    int or, oc;			/* Original dimensions of image */
-    int rows, cols;		/* Smallest powers of 2 >= number of rows & columns */
+    int rows, cols;		/* number of rows & columns */
     long totsize;		/* Total number of data points */
-    int halfrows, halfcols;
-    double *data[2];		/* Data structure containing real & complex values of FFT */
-    struct Option *op1, *op2, *op3;
-    struct GModule *module;
+    double (*data)[2];		/* Data structure containing real & complex values of FFT */
 
     G_gisinit(argv[0]);
 
@@ -51,188 +60,165 @@ int main(int argc, char *argv[])
 	_("Inverse Fast Fourier Transform (IFFT) for image processing.");
 
     /* define options */
-    op1 = G_define_option();
-    op1->key = "real_image";
-    op1->type = TYPE_STRING;
-    op1->required = YES;
-    op1->multiple = NO;
-    op1->gisprompt = "old,cell,raster";
-    op1->description = _("Input raster map (image fft, real part)");
-
-    op2 = G_define_option();
-    op2->key = "imaginary_image";
-    op2->type = TYPE_STRING;
-    op2->required = YES;
-    op2->multiple = NO;
-    op2->gisprompt = "old,cell,raster";
-    op2->description = _("Input raster map (image fft, imaginary part");
-
-    op3 = G_define_option();
-    op3->key = "output_image";
-    op3->type = TYPE_STRING;
-    op3->required = YES;
-    op3->multiple = NO;
-    op3->gisprompt = "new,cell,raster";
-    op3->description = _("Output inverse raster map after IFFT");
+    opt.real = G_define_option();
+    opt.real->key = "real_image";
+    opt.real->type = TYPE_STRING;
+    opt.real->required = YES;
+    opt.real->multiple = NO;
+    opt.real->gisprompt = "old,cell,raster";
+    opt.real->description = _("Input raster map (image fft, real part)");
+
+    opt.imag = G_define_option();
+    opt.imag->key = "imaginary_image";
+    opt.imag->type = TYPE_STRING;
+    opt.imag->required = YES;
+    opt.imag->multiple = NO;
+    opt.imag->gisprompt = "old,cell,raster";
+    opt.imag->description = _("Input raster map (image fft, imaginary part");
+
+    opt.orig = G_define_option();
+    opt.orig->key = "output_image";
+    opt.orig->type = TYPE_STRING;
+    opt.orig->required = YES;
+    opt.orig->multiple = NO;
+    opt.orig->gisprompt = "new,cell,raster";
+    opt.orig->description = _("Output inverse raster map after IFFT");
 
     /*call parser */
     if (G_parser(argc, argv))
 	exit(EXIT_FAILURE);
 
-    strcpy(Cellmap_real, op1->answer);
-    strcpy(Cellmap_imag, op2->answer);
-    strcpy(Cellmap_orig, op3->answer);
+    Cellmap_real = opt.real->answer;
+    Cellmap_imag = opt.imag->answer;
+    Cellmap_orig = opt.orig->answer;
 
     /* open input raster map */
-    if ((realmapset = G_find_cell(Cellmap_real, "")) == NULL)
-	G_fatal_error(_("Unable to find the real-image map <%s>"),
-		      Cellmap_real);
-
-    if ((realfp =
-	 G_fopen_old_misc("cell_misc", "fftreal", Cellmap_real,
-			  realmapset)) == NULL)
-	G_fatal_error(_("Unable to open real-image in the cell_misc directory.\nInput map probably wasn't created by i.fft"));
+    realfd = G_open_cell_old(Cellmap_real, "");
+    if (realfd < 0)
+	G_fatal_error(_("Unable to open real input map <%s>"), Cellmap_real);
 
-    if ((imagmapset = G_find_cell(Cellmap_imag, "")) == NULL)
-	G_fatal_error(_("Unable to find the imaginary-image <%s>"),
-		      Cellmap_imag);
-
-    if ((imagfp =
-	 G_fopen_old_misc("cell_misc", "fftimag", Cellmap_imag,
-			  imagmapset)) == NULL)
-	G_fatal_error(_("Unable to open imaginary-image in the cell_misc directory.\nInput map probably wasn't created by i.fft"));
+    imagfd = G_open_cell_old(Cellmap_imag, "");
+    if (imagfd < 0)
+	G_fatal_error(_("Unable to open imaginary input map <%s>"), Cellmap_imag);
 
     /* get and compare the original window data */
-    get_orig_window(&orig_wind, realmapset, imagmapset);
+    G_get_cellhd(Cellmap_real, "", &realhead);
+    G_get_cellhd(Cellmap_imag, "", &imaghead);
+
+    if (realhead.proj   != imaghead.proj   ||
+	realhead.zone   != imaghead.zone   ||
+	realhead.north  != imaghead.north  ||
+	realhead.south  != imaghead.south  ||
+	realhead.east   != imaghead.east   ||
+	realhead.west   != imaghead.west   ||
+	realhead.ew_res != imaghead.ew_res ||
+	realhead.ns_res != imaghead.ns_res)
+	G_fatal_error(_("The real and imaginary original windows did not match."));
 
-    or = orig_wind.rows;
-    oc = orig_wind.cols;
-    G_get_cellhd(Cellmap_real, realmapset, &realhead);
     G_set_window(&realhead);	/* set the window to the whole cell map */
 
     /* get the rows and columns in the current window */
     rows = G_window_rows();
     cols = G_window_cols();
     totsize = rows * cols;
-    halfrows = rows / 2;
-    halfcols = cols / 2;
-
-    G_message(_("Power 2 values : [%d] rows [%d] columns."), rows, cols);
 
     /* Allocate appropriate memory for the structure containing
        the real and complex components of the FFT.  DATA[0] will
        contain the real, and DATA[1] the complex component.
      */
-    data[0] = (double *)G_malloc((rows * cols) * sizeof(double));
-    data[1] = (double *)G_malloc((rows * cols) * sizeof(double));
+    data = G_malloc(rows * cols * 2 * sizeof(double));
+
+    /* allocate the space for one row of cell map data */
+    cell_real = G_allocate_d_raster_buf();
+    cell_imag = G_allocate_d_raster_buf();
 
-    /* Initialize real & complex components to zero */
+#define C(i, j) ((i) * cols + (j))
+
+    /* Read in cell map values */
     G_message(_("Reading the raster maps..."));
-    {
-	fread((char *)data[0], sizeof(double), totsize, realfp);
-	fread((char *)data[1], sizeof(double), totsize, imagfp);
+    for (i = 0; i < rows; i++) {
+	if (G_get_d_raster_row(realfd, cell_real, i) < 0)
+	    G_fatal_error(_("Error while reading real input raster map."));
+	if (G_get_d_raster_row(imagfd, cell_imag, i) < 0)
+	    G_fatal_error(_("Error while reading imaginary input raster map."));
+	for (j = 0; j < cols; j++) {
+	    data[C(i, j)][0] = cell_real[j];
+	    data[C(i, j)][1] = cell_imag[j];
+	}
     }
 
+    /* close input cell maps */
+    G_close_cell(realfd);
+    G_close_cell(imagfd);
+
     /* Read in cell map values */
     G_message(_("Masking the raster maps..."));
     maskfd = G_maskfd();
-    if (maskfd >= 0)
+    if (maskfd >= 0) {
 	maskbuf = G_allocate_cell_buf();
 
-    if (maskfd >= 0) {
 	for (i = 0; i < rows; i++) {
-	    double *data0, *data1;
-
-	    data0 = data[0] + i * cols;
-	    data1 = data[1] + i * cols;
 	    G_get_map_row(maskfd, maskbuf, i);
-	    for (j = 0; j < cols; j++, data0++, data1++) {
-		if (maskbuf[j] == (CELL) 0) {
-		    *(data0) = 0.0;
-		    *(data1) = 0.0;
+	    for (j = 0; j < cols; j++) {
+		if (maskbuf[j] == 0) {
+		    data[C(i, j)][0] = 0.0;
+		    data[C(i, j)][1] = 0.0;
 		}
 	    }
 	}
-    }
 
-    G_message(_("Rotating data arrays..."));
-    /* rotate the data array for standard display */
-    for (i = 0; i < rows; i++) {
-	double temp;
-
-	for (j = 0; j < halfcols; j++) {
-	    temp = *(data[0] + i * cols + j);
-	    *(data[0] + i * cols + j) = *(data[0] + i * cols + j + halfcols);
-	    *(data[0] + i * cols + j + halfcols) = temp;
-	    temp = *(data[1] + i * cols + j);
-	    *(data[1] + i * cols + j) = *(data[1] + i * cols + j + halfcols);
-	    *(data[1] + i * cols + j + halfcols) = temp;
-	}
+	G_close_cell(maskfd);
+	G_free(maskbuf);
     }
-    for (i = 0; i < halfrows; i++) {
-	double temp;
 
-	for (j = 0; j < cols; j++) {
-	    temp = *(data[0] + i * cols + j);
-	    *(data[0] + i * cols + j) =
-		*(data[0] + (i + halfrows) * cols + j);
-	    *(data[0] + (i + halfrows) * cols + j) = temp;
-	    temp = *(data[1] + i * cols + j);
-	    *(data[1] + i * cols + j) =
-		*(data[1] + (i + halfrows) * cols + j);
-	    *(data[1] + (i + halfrows) * cols + j) = temp;
-	}
-    }
+#define SWAP1(a, b)				\
+    do {					\
+	double temp = (a);			\
+	(a) = (b);				\
+	(b) = temp;				\
+    } while (0)
 
+#define SWAP2(a, b)				\
+    do {					\
+	SWAP1(data[(a)][0], data[(b)][0]);	\
+	SWAP1(data[(a)][1], data[(b)][1]);	\
+    } while (0)
 
-    /* close input cell maps and release the row buffers */
-    fclose(realfp);
-    fclose(imagfp);
-    if (maskfd >= 0) {
-	G_close_cell(maskfd);
-	G_free(maskbuf);
-    }
+    /* rotate the data array for standard display */
+    G_message(_("Rotating data arrays..."));
+    for (i = 0; i < rows; i++)
+	for (j = 0; j < cols / 2; j++)
+	    SWAP2(C(i, j), C(i, j + cols / 2));
+    for (i = 0; i < rows / 2; i++)
+	for (j = 0; j < cols; j++)
+	    SWAP2(C(i, j), C(i + rows / 2, j));
 
     /* perform inverse FFT */
     G_message(_("Starting Inverse FFT..."));
-    fft(1, data, totsize, cols, rows);
+    fft2(1, data, totsize, cols, rows);
     G_message(_("Inverse FFT completed..."));
 
-    /* set up a window for the transform cell map */
-    G_set_window(&orig_wind);
-
-    /* open the output cell map and allocate a cell row buffer */
-    if ((outputfd = G_open_cell_new(Cellmap_orig)) < 0)
-	G_fatal_error(_("Unable to open output file."));
-
-    cell_row = G_allocate_cell_buf();
+    /* open the output cell map */
+    if ((outputfd = G_open_fp_cell_new(Cellmap_orig)) < 0)
+	G_fatal_error(_("Unable to open output map <%s>"), Cellmap_orig);
 
     /* Write out result to a new cell map */
-    G_message(_("Writing data to file..."));
-    for (i = 0; i < or; i++) {
-	for (j = 0; j < oc; j++) {
-	    *(cell_row + j) = (CELL) (*(data[0] + i * cols + j) + 0.5);
-	}
-	G_put_raster_row(outputfd, cell_row, CELL_TYPE);
+    G_message(_("Writing output map..."));
+    for (i = 0; i < rows; i++) {
+	for (j = 0; j < cols; j++)
+	    cell_real[j] = data[C(i, j)][0];
+	G_put_d_raster_row(outputfd, cell_real);
     }
+
     G_close_cell(outputfd);
 
-    G_free(cell_row);
-    {
-	struct Colors colors;
-	struct Range range;
-	CELL min, max;
-
-	/* make a real component color table */
-	G_read_range(Cellmap_orig, G_mapset(), &range);
-	G_get_range_min_max(&range, &min, &max);
-	G_make_grey_scale_colors(&colors, min, max);
-	G_write_colors(Cellmap_orig, G_mapset(), &colors);
-    }
+    G_free(cell_real);
+    G_free(cell_imag);
+
+    fft_colors(Cellmap_orig);
 
     /* Release memory resources */
-    G_free(data[0]);
-    G_free(data[1]);
+    G_free(data);
 
     G_done_msg(_("Transform successful."));
 

imagery/i.ifft/orig_wind.c

@@ -1,30 +0,0 @@
-#include <grass/gis.h>
-#include <grass/glocale.h>
-#include "globals.h"
-
-#define FFTWINDOW "fftwindow"
-
-
-int get_orig_window(struct Cell_head *hd, char *rmapset, char *imapset)
-{
-    struct Cell_head tmphd;
-    char buffer[100];
-
-    /* get the windows for both the real and imaginary parts */
-    sprintf(buffer, "cell_misc/%s", Cellmap_real);
-    G__get_window(hd, buffer, FFTWINDOW, rmapset);
-    sprintf(buffer, "cell_misc/%s", Cellmap_imag);
-    G__get_window(&tmphd, buffer, FFTWINDOW, imapset);
-
-    /* compare them */
-    if (hd->proj != tmphd.proj ||
-	hd->zone != tmphd.zone ||
-	hd->north != tmphd.north ||
-	hd->south != tmphd.south ||
-	hd->east != tmphd.east ||
-	hd->west != tmphd.west ||
-	hd->ew_res != tmphd.ew_res || hd->ns_res != tmphd.ns_res)
-	G_fatal_error(_("The real and imaginary original windows did not match."));
-
-    return 0;
-}