Changeset 16107 for trunk/Ohana/src/opihi/cmd.data/fft2d.c

Timestamp:

Jan 17, 2008, 8:41:08 AM (19 years ago)

Author:

eugene

Message:

converting fft and fftn to code from Jones, with Magnier version of multi-dim

File:

• trunk/Ohana/src/opihi/cmd.data/fft2d.c (modified) (3 diffs)

trunk/Ohana/src/opihi/cmd.data/fft2d.c

@@ -3 +3 @@
 int fft2d (int argc, char **argv) {
   
-  int i, N, Nx, Ny, Naxis[2];
-  int Npix, ZeroImaginary, isign;
-  float *t1, *t2, *out, *temp;
+  int ZeroImaginary, forward;
+  int N, Nx, Ny;
   Buffer *Ire, *Iim, *Ore, *Oim;;
 
-  isign = 1;
+  forward = TRUE;
   if ((N = get_argument (argc, argv, "-inverse"))) {
     remove_argument (N, &argc, argv);
-    isign = -1;
+    forward = FALSE;
   }
 
@@ -28 +27 @@
   }    
   if ((Ire = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
+
+  Nx = Ire[0].header.Naxis[0];
+  Ny = Ire[0].header.Naxis[1];
+
+  // check input image dimensions (match lengths? binary lengths?)
+  if (!ZeroImaginary) {
+    if ((Nx != Iim[0].header.Naxis[0]) || 
+        (Ny != Iim[0].header.Naxis[1])) {
+      gprint (GP_ERR, "image size mismatch\n");
+      return (FALSE);
+    }
+  }
+  if (!IsBinary (Nx, NULL)) {
+    gprint (GP_ERR, "Nx is not a binary number!\n");
+    return (FALSE);
+  }
+  if (!IsBinary (Ny, NULL)) {
+    gprint (GP_ERR, "Ny is not a binary number!\n");
+    return (FALSE);
+  }
+
   if ((Ore = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE);
   if ((Oim = SelectBuffer (argv[5], ANYBUFFER, TRUE)) == NULL) return (FALSE);
@@ -36 +56 @@
   gfits_free_matrix (&Oim[0].matrix);
   gfits_free_header (&Oim[0].header);
+  
+  /* fix up output headers (real) & allocate data buffer */
+  CreateBuffer (Ore, Nx, Ny, -32, 0.0, 1.0);
+  CreateBuffer (Oim, Nx, Ny, -32, 0.0, 1.0);
 
-  /* get image dimensions, check value */
-  Npix = Ire[0].header.Naxis[0]*Ire[0].header.Naxis[1];
-  Nx = Ire[0].header.Naxis[0];
-  Ny = Ire[0].header.Naxis[1];
-  Naxis[0] = Ny; Naxis[1] = Nx;
-  if (!IsBinaryOld (Npix)) {
-    gprint (GP_ERR, "dimensions are not binary!\n");
-    return (FALSE);
+  gfits_copy_header (&Ire[0].header, &Ore[0].header);
+  gfits_copy_header (&Ire[0].header, &Oim[0].header);
+
+  // copy data to output buffers (fft is done in place)
+  memcpy (Ore[0].matrix.buffer, Ire[0].matrix.buffer, Nx*Ny*sizeof(float));
+
+  if (ZeroImaginary) {
+    memset (Oim[0].matrix.buffer, 0, Nx*Ny*sizeof(float));
+  } else {
+    memcpy (Oim[0].matrix.buffer, Iim[0].matrix.buffer, Nx*Ny*sizeof(float));
   }
-  
-  /* create working space */
-  t1 = (float *) Ire[0].matrix.buffer;
-  ALLOCATE (temp, float, 2*Npix);
-  out = temp;
 
-  /* copy input to working space */
-  if (ZeroImaginary) {
-    for (i = 0; i < Npix; i++, t1++) {
-      *out = *t1;
-      out++;
-      *out = 0;
-      out++;
-    }
-  } else {
-    t2 = (float *) Iim[0].matrix.buffer;
-    for (i = 0; i < Npix; i++, t1++, t2++) {
-      *out = *t1;
-      out++;
-      *out = *t2;
-      out++;
-    }
-  } 
-    
   /* run the fft */
-  fftN (temp, Naxis, 2, isign);
-
-  /* fix up output headers (real) */
-  gfits_copy_header (&Ire[0].header, &Ore[0].header);
-  gfits_modify (&Ore[0].header, "NAXIS1", "%d", 1, Nx);
-  gfits_modify (&Ore[0].header, "NAXIS2", "%d", 1, Ny);
-  Ore[0].header.Naxis[0] = Nx;
-  Ore[0].header.Naxis[1] = Ny;
-  Ore[0].bitpix = Ire[0].bitpix;
-  Ore[0].unsign = Ire[0].unsign;
-  Ore[0].bscale = Ire[0].bscale;
-  Ore[0].bzero  = Ire[0].bzero;
-  gfits_create_matrix (&Ore[0].header, &Ore[0].matrix);
-
-  /* fix up output headers (imaginary) */
-  gfits_copy_header (&Ire[0].header, &Oim[0].header);
-  gfits_modify (&Oim[0].header, "NAXIS1", "%d", 1, Nx);
-  gfits_modify (&Oim[0].header, "NAXIS2", "%d", 1, Ny);
-  Oim[0].header.Naxis[0] = Nx;
-  Oim[0].header.Naxis[1] = Ny;
-  Oim[0].bitpix = Ire[0].bitpix;
-  Oim[0].unsign = Ire[0].unsign;
-  Oim[0].bscale = Ire[0].bscale;
-  Oim[0].bzero  = Ire[0].bzero;
-  gfits_create_matrix (&Oim[0].header, &Oim[0].matrix);
-
-  /* move data from working space to output buffers */
-  out = temp;
-  t1 = (float *) Ore[0].matrix.buffer;
-  t2 = (float *) Oim[0].matrix.buffer;
-  for (i = 0; i < Npix; i++, t1++, t2++) {
-    *t1 = *out / Npix;
-    out ++;
-    *t2 = *out / Npix;
-    out ++;
-  }    
-
-  free (temp);
+  fftND ((float *)Ore[0].matrix.buffer, (float *)Oim[0].matrix.buffer, 2, Ire[0].header.Naxis, forward);
 
   return (TRUE);