[GRASS-SVN] r34353 - in grass/trunk/imagery: i.fft i.ifft
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Nov 17 20:41:17 EST 2008
Author: glynn
Date: 2008-11-17 20:41:17 -0500 (Mon, 17 Nov 2008)
New Revision: 34353
Removed:
grass/trunk/imagery/i.fft/fft_colors.c
grass/trunk/imagery/i.fft/globals.h
grass/trunk/imagery/i.fft/local_proto.h
grass/trunk/imagery/i.fft/orig_wind.c
grass/trunk/imagery/i.fft/save_fft.c
grass/trunk/imagery/i.ifft/globals.h
grass/trunk/imagery/i.ifft/local_proto.h
grass/trunk/imagery/i.ifft/orig_wind.c
Modified:
grass/trunk/imagery/i.fft/Makefile
grass/trunk/imagery/i.fft/i.fft.html
grass/trunk/imagery/i.fft/main.c
grass/trunk/imagery/i.ifft/Makefile
grass/trunk/imagery/i.ifft/i.ifft.html
grass/trunk/imagery/i.ifft/main.c
Log:
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
Modified: grass/trunk/imagery/i.fft/Makefile
===================================================================
--- grass/trunk/imagery/i.fft/Makefile 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.fft/Makefile 2008-11-18 01:41:17 UTC (rev 34353)
@@ -2,8 +2,8 @@
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
Deleted: grass/trunk/imagery/i.fft/fft_colors.c
===================================================================
--- grass/trunk/imagery/i.fft/fft_colors.c 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.fft/fft_colors.c 2008-11-18 01:41:17 UTC (rev 34353)
@@ -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;
-}
Deleted: grass/trunk/imagery/i.fft/globals.h
===================================================================
--- grass/trunk/imagery/i.fft/globals.h 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.fft/globals.h 2008-11-18 01:41:17 UTC (rev 34353)
@@ -1,6 +0,0 @@
-#ifndef __L_GLOBALS_H__
-#define __L_GLOBALS_H__
-
-extern char Cellmap_real[], Cellmap_imag[];
-
-#endif /* __L_GLOBALS_H__ */
Modified: grass/trunk/imagery/i.fft/i.fft.html
===================================================================
--- grass/trunk/imagery/i.fft/i.fft.html 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.fft/i.fft.html 2008-11-18 01:41:17 UTC (rev 34353)
@@ -8,35 +8,20 @@
<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>
Deleted: grass/trunk/imagery/i.fft/local_proto.h
===================================================================
--- grass/trunk/imagery/i.fft/local_proto.h 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.fft/local_proto.h 2008-11-18 01:41:17 UTC (rev 34353)
@@ -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__ */
Modified: grass/trunk/imagery/i.fft/main.c
===================================================================
--- grass/trunk/imagery/i.fft/main.c 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.fft/main.c 2008-11-18 01:41:17 UTC (rev 34353)
@@ -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 @@
_("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");
+ 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");
- 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");
+ 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");
- 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");
+ 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");
- 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 */
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
- strcpy(Cellmap_orig, op1->answer);
- strcpy(Cellmap_real, op2->answer);
- strcpy(Cellmap_imag, op3->answer);
+ Cellmap_orig = opt.orig->answer;
+ Cellmap_real = opt.real->answer;
+ Cellmap_imag = opt.imag->answer;
- /* open input cell map */
- if ((inmapset = G_find_cell(Cellmap_orig, "")) == NULL)
- G_fatal_error(_("Raster map <%s> not found"), Cellmap_orig);
-
- 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));
+ data = G_malloc(rows * cols * 2 * 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;
- }
- }
-
/* 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 */
+ 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);
- /* 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();
+#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."));
Deleted: grass/trunk/imagery/i.fft/orig_wind.c
===================================================================
--- grass/trunk/imagery/i.fft/orig_wind.c 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.fft/orig_wind.c 2008-11-18 01:41:17 UTC (rev 34353)
@@ -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;
-}
Deleted: grass/trunk/imagery/i.fft/save_fft.c
===================================================================
--- grass/trunk/imagery/i.fft/save_fft.c 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.fft/save_fft.c 2008-11-18 01:41:17 UTC (rev 34353)
@@ -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;
-}
Modified: grass/trunk/imagery/i.ifft/Makefile
===================================================================
--- grass/trunk/imagery/i.ifft/Makefile 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.ifft/Makefile 2008-11-18 01:41:17 UTC (rev 34353)
@@ -2,8 +2,8 @@
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
Deleted: grass/trunk/imagery/i.ifft/globals.h
===================================================================
--- grass/trunk/imagery/i.ifft/globals.h 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.ifft/globals.h 2008-11-18 01:41:17 UTC (rev 34353)
@@ -1,6 +0,0 @@
-#ifndef __L_GLOBALS_H__
-#define __L_GLOBALS_H__
-
-extern char Cellmap_real[], Cellmap_imag[];
-
-#endif /* __L_GLOBALS_H__ */
Modified: grass/trunk/imagery/i.ifft/i.ifft.html
===================================================================
--- grass/trunk/imagery/i.ifft/i.ifft.html 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.ifft/i.ifft.html 2008-11-18 01:41:17 UTC (rev 34353)
@@ -20,22 +20,6 @@
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
Deleted: grass/trunk/imagery/i.ifft/local_proto.h
===================================================================
--- grass/trunk/imagery/i.ifft/local_proto.h 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.ifft/local_proto.h 2008-11-18 01:41:17 UTC (rev 34353)
@@ -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 *);
Modified: grass/trunk/imagery/i.ifft/main.c
===================================================================
--- grass/trunk/imagery/i.ifft/main.c 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.ifft/main.c 2008-11-18 01:41:17 UTC (rev 34353)
@@ -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 @@
_("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)");
+ 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)");
- 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");
+ 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");
- 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.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);
+ realfd = G_open_cell_old(Cellmap_real, "");
+ if (realfd < 0)
+ G_fatal_error(_("Unable to open real input 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"));
+ imagfd = G_open_cell_old(Cellmap_imag, "");
+ if (imagfd < 0)
+ G_fatal_error(_("Unable to open imaginary input map <%s>"), Cellmap_imag);
- 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"));
-
/* 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);
- or = orig_wind.rows;
- oc = orig_wind.cols;
- G_get_cellhd(Cellmap_real, realmapset, &realhead);
+ 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."));
+
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));
- /* Initialize real & complex components to zero */
+ /* allocate the space for one row of cell map data */
+ 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 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 */
+ if ((outputfd = G_open_fp_cell_new(Cellmap_orig)) < 0)
+ G_fatal_error(_("Unable to open output map <%s>"), Cellmap_orig);
- /* 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();
-
/* 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;
+ G_free(cell_real);
+ G_free(cell_imag);
- /* 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);
- }
+ fft_colors(Cellmap_orig);
/* Release memory resources */
- G_free(data[0]);
- G_free(data[1]);
+ G_free(data);
G_done_msg(_("Transform successful."));
Deleted: grass/trunk/imagery/i.ifft/orig_wind.c
===================================================================
--- grass/trunk/imagery/i.ifft/orig_wind.c 2008-11-18 01:37:56 UTC (rev 34352)
+++ grass/trunk/imagery/i.ifft/orig_wind.c 2008-11-18 01:41:17 UTC (rev 34353)
@@ -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;
-}
More information about the grass-commit
mailing list