Index: /trunk/psLib/src/dataManip/psFFT.c
===================================================================
--- /trunk/psLib/src/dataManip/psFFT.c	(revision 801)
+++ /trunk/psLib/src/dataManip/psFFT.c	(revision 801)
@@ -0,0 +1,466 @@
+/** @file  psFFT.c
+ *
+ *  @brief Contains FFT transforms functions
+ *
+ *  @author Robert DeSonia, MHPCC
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-05-28 03:15:40 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
+
+#include <unistd.h>
+#include <stdbool.h>
+#include <string.h>
+#include <complex.h>
+
+#include "psFFT.h"
+#include "psError.h"
+#include "psMemory.h"
+#include "psLogMsg.h"
+
+#define P_FFTW_PLAN_RIGOR FFTW_ESTIMATE
+
+static bool p_fftwWisdomImported = false;
+
+psImage* psImageFFT(const psImage *image, int flags)
+{
+    static fftwf_plan lastForwardPlan = NULL;
+    static int lastForwardPlanType = 0;
+    static int lastForwardPlanRows = 0;
+    static int lastForwardPlanCols = 0;
+    static fftwf_plan lastReversePlan = NULL;
+    static int lastReversePlanResultType = 0;
+    static int lastReversePlanRows = 0;
+    static int lastReversePlanCols = 0;
+    unsigned int numCols;
+    unsigned int numRows;
+    psElemType type;
+    fftwf_plan plan;
+    psImage* out = NULL;
+
+    /* got good image data? */
+    if (image==NULL) {
+        return NULL;
+    }
+
+    if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) {
+        psError(__func__,"Can not perform a forward FFT to real result.");
+        return NULL;
+    }
+
+    type = image->type.type;
+
+    if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) {
+        psError(__func__,"Input image must be a 32-bit float or complex image (type=%d)",
+                type);
+        return NULL;
+    }
+
+    if ( (type != PS_TYPE_C32) && ((flags & PS_FFT_REVERSE) != 0) ) {
+        psError(__func__,"Input image must be complex image for reverse FFT (type=%d).",
+                type);
+        return NULL;
+
+    }
+
+    /* make sure the system-level wisdom information is imported. */
+    if (! p_fftwWisdomImported) {
+        fftw_import_system_wisdom();
+        p_fftwWisdomImported = true;
+    }
+
+    numRows = image->numRows;
+    numCols = (image->numCols-1)*2;
+    // n.b., numCols needs adjustment for c->r, see below.
+
+    if (flags & PS_FFT_REVERSE) { // reverse transform
+        if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r)
+            // Assuming image is from a r2c transform.
+            // (FFTW only uses nx/2+1 for r->c transform results)
+            numCols = (image->numCols-1)*2;
+
+            out = psImageAlloc(numCols,numRows,PS_TYPE_F32);
+
+            if (lastReversePlanRows != numRows || lastReversePlanCols != numCols
+                    || lastReversePlanResultType != PS_TYPE_F32
+                    || lastReversePlan == NULL) {
+                lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows,
+                                                        (fftwf_complex*)image->data.C32[0],out->data.F32[0],
+                                                        P_FFTW_PLAN_RIGOR);
+            }
+        } else { // complex result desired
+            out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+            if (lastReversePlanRows != numRows || lastReversePlanCols != numCols
+                    || lastReversePlanResultType != PS_TYPE_C32
+                    || lastReversePlan == NULL) {
+                lastReversePlan = fftwf_plan_dft_2d(numCols, numRows,
+                                                    (fftwf_complex*)image->data.C32[0],
+                                                    (fftwf_complex*)out->data.C32[0],FFTW_BACKWARD,
+                                                    P_FFTW_PLAN_RIGOR);
+            }
+        }
+        plan = lastReversePlan;
+    } else { // forward transform
+        if (type == PS_TYPE_F32) { // real data input
+            // (FFTW only uses nx/2+1 for r->c transform results)
+            out = psImageAlloc(numCols/2+1,numRows,PS_TYPE_C32);
+
+            if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols
+                    || lastForwardPlanType != type || lastForwardPlan == NULL)  {
+                /* need new plan created */
+                lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows,
+                                                        image->data.F32[0],(fftwf_complex*)out->data.C32[0],
+                                                        P_FFTW_PLAN_RIGOR);
+            }
+        } else if (type == PS_TYPE_C32) { // complex data input
+            out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+
+            if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols
+                    || lastForwardPlanType != type || lastForwardPlan == NULL)  {
+                lastForwardPlan = fftwf_plan_dft_2d(numCols, numRows,
+                                                    (fftwf_complex*)image->data.C32[0],
+                                                    (fftwf_complex*)out->data.C32[0],
+                                                    FFTW_FORWARD,P_FFTW_PLAN_RIGOR);
+            }
+        }
+        plan = lastForwardPlan;
+    }
+
+    /* check if a plan exists now*/
+    if (plan == NULL) {
+        psError(__func__,"Failed to create FFTW plan.");
+        psImageFree(out);
+        return NULL;
+    }
+
+    /* finally, call FFTW with the plan made above */
+    fftwf_execute(plan);
+
+    return out;
+
+}
+
+
+psImage *psImageReal(psImage *out, const psImage* in)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numCols = in->numCols;
+    numRows = in->numRows;
+
+    /* if not a complex number, this is logically just a copy */
+    if (! IS_PSELEMTYPE_COMPLEX(type)) {
+        // Warn user, as this is probably not expected
+        psLogMsg(__func__,PS_LOG_WARN,"Real portion of a non-Complex type called called for. "
+                 "Just an image copy was performed.");
+        return psImageCopy(out,in,type);
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outRow;
+        psC32* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F32[row];
+            inRow = in->data.C32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = crealf(inRow[col]);
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psF64* outRow;
+        psC64* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F64[row];
+            inRow = in->data.C64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = creal(inRow[col]);
+            }
+        }
+    } else {
+        psError(__func__,"Can not extract real component from given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImageImaginary(psImage *out, const psImage* in)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numCols = in->numCols;
+    numRows = in->numRows;
+
+    /* if not a complex number, this is logically just zeroed image of same size */
+    if (! IS_PSELEMTYPE_COMPLEX(type)) {
+        // Warn user, as this is probably not expected
+        psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. "
+                 "A zero image was returned.");
+        out = psImageRealloc(out,numCols,numRows,type);
+        memset(out->data.v[0],0,PSELEMTYPE_SIZEOF(type)*numCols*numRows);
+        return out;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outRow;
+        psC32* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F32[row];
+            inRow = in->data.C32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = cimagf(inRow[col]);
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psF64* outRow;
+        psC64* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F64[row];
+            inRow = in->data.C64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = cimag(inRow[col]);
+            }
+        }
+    } else {
+        psError(__func__,"Can not extract imaginary component from given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImageComplex(psImage *real, const psImage *imag)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+    psImage* out = NULL;
+
+
+    if (real == NULL || imag == NULL) {
+        return NULL;
+    }
+
+    type = real->type.type;
+    numCols = real->numCols;
+    numRows = real->numRows;
+
+    if (imag->type.type != type) {
+        psError(__func__,"The inputs to psImageComplex must be the same type.");
+        return NULL;
+    }
+
+    if (IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__,"The inputs to psImageComplex can not be complex.");
+        return NULL;
+    }
+
+    if (type != PS_TYPE_F32 && type != PS_TYPE_F64) {
+        psError(__func__,"The input type to psImageComplex must be a floating point.");
+        return NULL;
+    }
+
+    if (type == PS_TYPE_F32) {
+        psC32* outRow;
+        psF32* realRow;
+        psF32* imagRow;
+
+        out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.C32[row];
+            realRow = real->data.F32[row];
+            imagRow = imag->data.F32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = realRow[col] + I*imagRow[col];
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psC64* outRow;
+        psF64* realRow;
+        psF64* imagRow;
+
+        out = psImageAlloc(numCols,numRows,PS_TYPE_C64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.C64[row];
+            realRow = real->data.F64[row];
+            imagRow = imag->data.F64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = realRow[col] + I*imagRow[col];
+            }
+        }
+    } else {
+        psError(__func__,"Can not merge real and imaginary portions for given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImageConjugate(psImage *out, const psImage *in)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numCols = in->numCols;
+    numRows = in->numRows;
+
+    /* if not a complex number, this is logically just a image copy */
+    if (! IS_PSELEMTYPE_COMPLEX(type)) {
+        // Warn user, as this is probably not expected
+        psLogMsg(__func__,PS_LOG_WARN,"Complex Conjugate of a non-Complex type called for. "
+                 "Image copy was performed instead.");
+        return psImageCopy(out,in,type);
+    }
+
+    if (type == PS_TYPE_C32) {
+        psC32* outRow;
+        psC32* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_C32);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.C32[row];
+            inRow = in->data.C32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = crealf(inRow[col]) - I*cimagf(inRow[col]);
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psC64* outRow;
+        psC64* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_C64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.C64[row];
+            inRow = in->data.C64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = creal(inRow[col]) - I*cimag(inRow[col]);
+            }
+        }
+    } else {
+        psError(__func__,"Can not compute complex conjugate for given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImagePowerSpectrum(const psImage* in)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+    int numElementsSquared;
+    psImage* out = NULL;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numCols = in->numCols;
+    numRows = in->numRows;
+    numElementsSquared = numCols*numCols*numRows*numRows;
+
+    /* if not a complex number, this is not implemented */
+    if (! IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__,"Power Spectrum for non-complex inputs is not implemented.");
+        return NULL;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outRow;
+        psC32* inRow;
+        psF32 real;
+        psF32 imag;
+
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F32[row];
+            inRow = in->data.C32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                real = crealf(inRow[col]);
+                imag = cimagf(inRow[col]);
+                outRow[col] = (real*real+imag*imag)/numElementsSquared;
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psF64* outRow;
+        psC64* inRow;
+        psF64 real;
+        psF64 imag;
+
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F64[row];
+            inRow = in->data.C64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                real = crealf(inRow[col]);
+                imag = cimagf(inRow[col]);
+                outRow[col] = real*real+imag*imag/numElementsSquared;
+            }
+        }
+    } else {
+        psError(__func__,"Can not power spectrum for given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+
+}
Index: /trunk/psLib/src/dataManip/psFFT.h
===================================================================
--- /trunk/psLib/src/dataManip/psFFT.h	(revision 801)
+++ /trunk/psLib/src/dataManip/psFFT.h	(revision 801)
@@ -0,0 +1,43 @@
+/** @file  psFFT.h
+ *
+ *  @brief Contains FFT transforms functions
+ *
+ *  @author Robert DeSonia, MHPCC
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-05-28 03:15:40 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
+
+#ifndef PS_FFT_H
+#define PS_FFT_H
+
+#include <fftw3.h>
+
+#include "psImage.h"
+#include "psVector.h"
+
+/** Details on FFT implementation (private). */
+
+
+typedef enum {
+    PS_FFT_FORWARD = 0,                 ///< psImageFFT/psVectorFFT should perform a forward FFT.
+    PS_FFT_REVERSE = 1,                 ///< psImageFFT/psVectorFFT should perform a reverse FFT.
+    PS_FFT_REAL_RESULT = 2
+                         /**< psImageFFT/psVectorFFT should return a real image.  This is valid for
+                          *   only reverse FFT, i.e., the psImageFFT/psVectorFFT flag parameter
+                          *   should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT.
+                          */
+} psFftFlags;
+
+psImage *psImageFFT(const psImage *image, int flags);
+psImage *psImageReal(psImage *out, const psImage* in);
+psImage *psImageImaginary(psImage *out, const psImage* in);
+psImage *psImageComplex(psImage *real, const psImage *imag);
+psImage *psImageConjugate(psImage *out, const psImage *in);
+psImage *psImagePowerSpectrum(const psImage* in);
+
+
+#endif
+
Index: /trunk/psLib/src/dataManip/psVectorFFT.c
===================================================================
--- /trunk/psLib/src/dataManip/psVectorFFT.c	(revision 801)
+++ /trunk/psLib/src/dataManip/psVectorFFT.c	(revision 801)
@@ -0,0 +1,466 @@
+/** @file  psFFT.c
+ *
+ *  @brief Contains FFT transforms functions
+ *
+ *  @author Robert DeSonia, MHPCC
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-05-28 03:15:40 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
+
+#include <unistd.h>
+#include <stdbool.h>
+#include <string.h>
+#include <complex.h>
+
+#include "psFFT.h"
+#include "psError.h"
+#include "psMemory.h"
+#include "psLogMsg.h"
+
+#define P_FFTW_PLAN_RIGOR FFTW_ESTIMATE
+
+static bool p_fftwWisdomImported = false;
+
+psImage* psImageFFT(const psImage *image, int flags)
+{
+    static fftwf_plan lastForwardPlan = NULL;
+    static int lastForwardPlanType = 0;
+    static int lastForwardPlanRows = 0;
+    static int lastForwardPlanCols = 0;
+    static fftwf_plan lastReversePlan = NULL;
+    static int lastReversePlanResultType = 0;
+    static int lastReversePlanRows = 0;
+    static int lastReversePlanCols = 0;
+    unsigned int numCols;
+    unsigned int numRows;
+    psElemType type;
+    fftwf_plan plan;
+    psImage* out = NULL;
+
+    /* got good image data? */
+    if (image==NULL) {
+        return NULL;
+    }
+
+    if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) {
+        psError(__func__,"Can not perform a forward FFT to real result.");
+        return NULL;
+    }
+
+    type = image->type.type;
+
+    if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) {
+        psError(__func__,"Input image must be a 32-bit float or complex image (type=%d)",
+                type);
+        return NULL;
+    }
+
+    if ( (type != PS_TYPE_C32) && ((flags & PS_FFT_REVERSE) != 0) ) {
+        psError(__func__,"Input image must be complex image for reverse FFT (type=%d).",
+                type);
+        return NULL;
+
+    }
+
+    /* make sure the system-level wisdom information is imported. */
+    if (! p_fftwWisdomImported) {
+        fftw_import_system_wisdom();
+        p_fftwWisdomImported = true;
+    }
+
+    numRows = image->numRows;
+    numCols = (image->numCols-1)*2;
+    // n.b., numCols needs adjustment for c->r, see below.
+
+    if (flags & PS_FFT_REVERSE) { // reverse transform
+        if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r)
+            // Assuming image is from a r2c transform.
+            // (FFTW only uses nx/2+1 for r->c transform results)
+            numCols = (image->numCols-1)*2;
+
+            out = psImageAlloc(numCols,numRows,PS_TYPE_F32);
+
+            if (lastReversePlanRows != numRows || lastReversePlanCols != numCols
+                    || lastReversePlanResultType != PS_TYPE_F32
+                    || lastReversePlan == NULL) {
+                lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows,
+                                                        (fftwf_complex*)image->data.C32[0],out->data.F32[0],
+                                                        P_FFTW_PLAN_RIGOR);
+            }
+        } else { // complex result desired
+            out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+            if (lastReversePlanRows != numRows || lastReversePlanCols != numCols
+                    || lastReversePlanResultType != PS_TYPE_C32
+                    || lastReversePlan == NULL) {
+                lastReversePlan = fftwf_plan_dft_2d(numCols, numRows,
+                                                    (fftwf_complex*)image->data.C32[0],
+                                                    (fftwf_complex*)out->data.C32[0],FFTW_BACKWARD,
+                                                    P_FFTW_PLAN_RIGOR);
+            }
+        }
+        plan = lastReversePlan;
+    } else { // forward transform
+        if (type == PS_TYPE_F32) { // real data input
+            // (FFTW only uses nx/2+1 for r->c transform results)
+            out = psImageAlloc(numCols/2+1,numRows,PS_TYPE_C32);
+
+            if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols
+                    || lastForwardPlanType != type || lastForwardPlan == NULL)  {
+                /* need new plan created */
+                lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows,
+                                                        image->data.F32[0],(fftwf_complex*)out->data.C32[0],
+                                                        P_FFTW_PLAN_RIGOR);
+            }
+        } else if (type == PS_TYPE_C32) { // complex data input
+            out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+
+            if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols
+                    || lastForwardPlanType != type || lastForwardPlan == NULL)  {
+                lastForwardPlan = fftwf_plan_dft_2d(numCols, numRows,
+                                                    (fftwf_complex*)image->data.C32[0],
+                                                    (fftwf_complex*)out->data.C32[0],
+                                                    FFTW_FORWARD,P_FFTW_PLAN_RIGOR);
+            }
+        }
+        plan = lastForwardPlan;
+    }
+
+    /* check if a plan exists now*/
+    if (plan == NULL) {
+        psError(__func__,"Failed to create FFTW plan.");
+        psImageFree(out);
+        return NULL;
+    }
+
+    /* finally, call FFTW with the plan made above */
+    fftwf_execute(plan);
+
+    return out;
+
+}
+
+
+psImage *psImageReal(psImage *out, const psImage* in)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numCols = in->numCols;
+    numRows = in->numRows;
+
+    /* if not a complex number, this is logically just a copy */
+    if (! IS_PSELEMTYPE_COMPLEX(type)) {
+        // Warn user, as this is probably not expected
+        psLogMsg(__func__,PS_LOG_WARN,"Real portion of a non-Complex type called called for. "
+                 "Just an image copy was performed.");
+        return psImageCopy(out,in,type);
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outRow;
+        psC32* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F32[row];
+            inRow = in->data.C32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = crealf(inRow[col]);
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psF64* outRow;
+        psC64* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F64[row];
+            inRow = in->data.C64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = creal(inRow[col]);
+            }
+        }
+    } else {
+        psError(__func__,"Can not extract real component from given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImageImaginary(psImage *out, const psImage* in)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numCols = in->numCols;
+    numRows = in->numRows;
+
+    /* if not a complex number, this is logically just zeroed image of same size */
+    if (! IS_PSELEMTYPE_COMPLEX(type)) {
+        // Warn user, as this is probably not expected
+        psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. "
+                 "A zero image was returned.");
+        out = psImageRealloc(out,numCols,numRows,type);
+        memset(out->data.v[0],0,PSELEMTYPE_SIZEOF(type)*numCols*numRows);
+        return out;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outRow;
+        psC32* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F32[row];
+            inRow = in->data.C32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = cimagf(inRow[col]);
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psF64* outRow;
+        psC64* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F64[row];
+            inRow = in->data.C64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = cimag(inRow[col]);
+            }
+        }
+    } else {
+        psError(__func__,"Can not extract imaginary component from given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImageComplex(psImage *real, const psImage *imag)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+    psImage* out = NULL;
+
+
+    if (real == NULL || imag == NULL) {
+        return NULL;
+    }
+
+    type = real->type.type;
+    numCols = real->numCols;
+    numRows = real->numRows;
+
+    if (imag->type.type != type) {
+        psError(__func__,"The inputs to psImageComplex must be the same type.");
+        return NULL;
+    }
+
+    if (IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__,"The inputs to psImageComplex can not be complex.");
+        return NULL;
+    }
+
+    if (type != PS_TYPE_F32 && type != PS_TYPE_F64) {
+        psError(__func__,"The input type to psImageComplex must be a floating point.");
+        return NULL;
+    }
+
+    if (type == PS_TYPE_F32) {
+        psC32* outRow;
+        psF32* realRow;
+        psF32* imagRow;
+
+        out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.C32[row];
+            realRow = real->data.F32[row];
+            imagRow = imag->data.F32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = realRow[col] + I*imagRow[col];
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psC64* outRow;
+        psF64* realRow;
+        psF64* imagRow;
+
+        out = psImageAlloc(numCols,numRows,PS_TYPE_C64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.C64[row];
+            realRow = real->data.F64[row];
+            imagRow = imag->data.F64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = realRow[col] + I*imagRow[col];
+            }
+        }
+    } else {
+        psError(__func__,"Can not merge real and imaginary portions for given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImageConjugate(psImage *out, const psImage *in)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numCols = in->numCols;
+    numRows = in->numRows;
+
+    /* if not a complex number, this is logically just a image copy */
+    if (! IS_PSELEMTYPE_COMPLEX(type)) {
+        // Warn user, as this is probably not expected
+        psLogMsg(__func__,PS_LOG_WARN,"Complex Conjugate of a non-Complex type called for. "
+                 "Image copy was performed instead.");
+        return psImageCopy(out,in,type);
+    }
+
+    if (type == PS_TYPE_C32) {
+        psC32* outRow;
+        psC32* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_C32);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.C32[row];
+            inRow = in->data.C32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = crealf(inRow[col]) - I*cimagf(inRow[col]);
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psC64* outRow;
+        psC64* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_C64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.C64[row];
+            inRow = in->data.C64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = creal(inRow[col]) - I*cimag(inRow[col]);
+            }
+        }
+    } else {
+        psError(__func__,"Can not compute complex conjugate for given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImagePowerSpectrum(const psImage* in)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+    int numElementsSquared;
+    psImage* out = NULL;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numCols = in->numCols;
+    numRows = in->numRows;
+    numElementsSquared = numCols*numCols*numRows*numRows;
+
+    /* if not a complex number, this is not implemented */
+    if (! IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__,"Power Spectrum for non-complex inputs is not implemented.");
+        return NULL;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outRow;
+        psC32* inRow;
+        psF32 real;
+        psF32 imag;
+
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F32[row];
+            inRow = in->data.C32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                real = crealf(inRow[col]);
+                imag = cimagf(inRow[col]);
+                outRow[col] = (real*real+imag*imag)/numElementsSquared;
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psF64* outRow;
+        psC64* inRow;
+        psF64 real;
+        psF64 imag;
+
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F64[row];
+            inRow = in->data.C64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                real = crealf(inRow[col]);
+                imag = cimagf(inRow[col]);
+                outRow[col] = real*real+imag*imag/numElementsSquared;
+            }
+        }
+    } else {
+        psError(__func__,"Can not power spectrum for given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+
+}
Index: /trunk/psLib/src/dataManip/psVectorFFT.h
===================================================================
--- /trunk/psLib/src/dataManip/psVectorFFT.h	(revision 801)
+++ /trunk/psLib/src/dataManip/psVectorFFT.h	(revision 801)
@@ -0,0 +1,43 @@
+/** @file  psFFT.h
+ *
+ *  @brief Contains FFT transforms functions
+ *
+ *  @author Robert DeSonia, MHPCC
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-05-28 03:15:40 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
+
+#ifndef PS_FFT_H
+#define PS_FFT_H
+
+#include <fftw3.h>
+
+#include "psImage.h"
+#include "psVector.h"
+
+/** Details on FFT implementation (private). */
+
+
+typedef enum {
+    PS_FFT_FORWARD = 0,                 ///< psImageFFT/psVectorFFT should perform a forward FFT.
+    PS_FFT_REVERSE = 1,                 ///< psImageFFT/psVectorFFT should perform a reverse FFT.
+    PS_FFT_REAL_RESULT = 2
+                         /**< psImageFFT/psVectorFFT should return a real image.  This is valid for
+                          *   only reverse FFT, i.e., the psImageFFT/psVectorFFT flag parameter
+                          *   should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT.
+                          */
+} psFftFlags;
+
+psImage *psImageFFT(const psImage *image, int flags);
+psImage *psImageReal(psImage *out, const psImage* in);
+psImage *psImageImaginary(psImage *out, const psImage* in);
+psImage *psImageComplex(psImage *real, const psImage *imag);
+psImage *psImageConjugate(psImage *out, const psImage *in);
+psImage *psImagePowerSpectrum(const psImage* in);
+
+
+#endif
+
Index: /trunk/psLib/src/fft/psVectorFFT.c
===================================================================
--- /trunk/psLib/src/fft/psVectorFFT.c	(revision 801)
+++ /trunk/psLib/src/fft/psVectorFFT.c	(revision 801)
@@ -0,0 +1,466 @@
+/** @file  psFFT.c
+ *
+ *  @brief Contains FFT transforms functions
+ *
+ *  @author Robert DeSonia, MHPCC
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-05-28 03:15:40 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
+
+#include <unistd.h>
+#include <stdbool.h>
+#include <string.h>
+#include <complex.h>
+
+#include "psFFT.h"
+#include "psError.h"
+#include "psMemory.h"
+#include "psLogMsg.h"
+
+#define P_FFTW_PLAN_RIGOR FFTW_ESTIMATE
+
+static bool p_fftwWisdomImported = false;
+
+psImage* psImageFFT(const psImage *image, int flags)
+{
+    static fftwf_plan lastForwardPlan = NULL;
+    static int lastForwardPlanType = 0;
+    static int lastForwardPlanRows = 0;
+    static int lastForwardPlanCols = 0;
+    static fftwf_plan lastReversePlan = NULL;
+    static int lastReversePlanResultType = 0;
+    static int lastReversePlanRows = 0;
+    static int lastReversePlanCols = 0;
+    unsigned int numCols;
+    unsigned int numRows;
+    psElemType type;
+    fftwf_plan plan;
+    psImage* out = NULL;
+
+    /* got good image data? */
+    if (image==NULL) {
+        return NULL;
+    }
+
+    if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) {
+        psError(__func__,"Can not perform a forward FFT to real result.");
+        return NULL;
+    }
+
+    type = image->type.type;
+
+    if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) {
+        psError(__func__,"Input image must be a 32-bit float or complex image (type=%d)",
+                type);
+        return NULL;
+    }
+
+    if ( (type != PS_TYPE_C32) && ((flags & PS_FFT_REVERSE) != 0) ) {
+        psError(__func__,"Input image must be complex image for reverse FFT (type=%d).",
+                type);
+        return NULL;
+
+    }
+
+    /* make sure the system-level wisdom information is imported. */
+    if (! p_fftwWisdomImported) {
+        fftw_import_system_wisdom();
+        p_fftwWisdomImported = true;
+    }
+
+    numRows = image->numRows;
+    numCols = (image->numCols-1)*2;
+    // n.b., numCols needs adjustment for c->r, see below.
+
+    if (flags & PS_FFT_REVERSE) { // reverse transform
+        if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r)
+            // Assuming image is from a r2c transform.
+            // (FFTW only uses nx/2+1 for r->c transform results)
+            numCols = (image->numCols-1)*2;
+
+            out = psImageAlloc(numCols,numRows,PS_TYPE_F32);
+
+            if (lastReversePlanRows != numRows || lastReversePlanCols != numCols
+                    || lastReversePlanResultType != PS_TYPE_F32
+                    || lastReversePlan == NULL) {
+                lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows,
+                                                        (fftwf_complex*)image->data.C32[0],out->data.F32[0],
+                                                        P_FFTW_PLAN_RIGOR);
+            }
+        } else { // complex result desired
+            out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+            if (lastReversePlanRows != numRows || lastReversePlanCols != numCols
+                    || lastReversePlanResultType != PS_TYPE_C32
+                    || lastReversePlan == NULL) {
+                lastReversePlan = fftwf_plan_dft_2d(numCols, numRows,
+                                                    (fftwf_complex*)image->data.C32[0],
+                                                    (fftwf_complex*)out->data.C32[0],FFTW_BACKWARD,
+                                                    P_FFTW_PLAN_RIGOR);
+            }
+        }
+        plan = lastReversePlan;
+    } else { // forward transform
+        if (type == PS_TYPE_F32) { // real data input
+            // (FFTW only uses nx/2+1 for r->c transform results)
+            out = psImageAlloc(numCols/2+1,numRows,PS_TYPE_C32);
+
+            if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols
+                    || lastForwardPlanType != type || lastForwardPlan == NULL)  {
+                /* need new plan created */
+                lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows,
+                                                        image->data.F32[0],(fftwf_complex*)out->data.C32[0],
+                                                        P_FFTW_PLAN_RIGOR);
+            }
+        } else if (type == PS_TYPE_C32) { // complex data input
+            out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+
+            if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols
+                    || lastForwardPlanType != type || lastForwardPlan == NULL)  {
+                lastForwardPlan = fftwf_plan_dft_2d(numCols, numRows,
+                                                    (fftwf_complex*)image->data.C32[0],
+                                                    (fftwf_complex*)out->data.C32[0],
+                                                    FFTW_FORWARD,P_FFTW_PLAN_RIGOR);
+            }
+        }
+        plan = lastForwardPlan;
+    }
+
+    /* check if a plan exists now*/
+    if (plan == NULL) {
+        psError(__func__,"Failed to create FFTW plan.");
+        psImageFree(out);
+        return NULL;
+    }
+
+    /* finally, call FFTW with the plan made above */
+    fftwf_execute(plan);
+
+    return out;
+
+}
+
+
+psImage *psImageReal(psImage *out, const psImage* in)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numCols = in->numCols;
+    numRows = in->numRows;
+
+    /* if not a complex number, this is logically just a copy */
+    if (! IS_PSELEMTYPE_COMPLEX(type)) {
+        // Warn user, as this is probably not expected
+        psLogMsg(__func__,PS_LOG_WARN,"Real portion of a non-Complex type called called for. "
+                 "Just an image copy was performed.");
+        return psImageCopy(out,in,type);
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outRow;
+        psC32* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F32[row];
+            inRow = in->data.C32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = crealf(inRow[col]);
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psF64* outRow;
+        psC64* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F64[row];
+            inRow = in->data.C64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = creal(inRow[col]);
+            }
+        }
+    } else {
+        psError(__func__,"Can not extract real component from given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImageImaginary(psImage *out, const psImage* in)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numCols = in->numCols;
+    numRows = in->numRows;
+
+    /* if not a complex number, this is logically just zeroed image of same size */
+    if (! IS_PSELEMTYPE_COMPLEX(type)) {
+        // Warn user, as this is probably not expected
+        psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. "
+                 "A zero image was returned.");
+        out = psImageRealloc(out,numCols,numRows,type);
+        memset(out->data.v[0],0,PSELEMTYPE_SIZEOF(type)*numCols*numRows);
+        return out;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outRow;
+        psC32* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F32[row];
+            inRow = in->data.C32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = cimagf(inRow[col]);
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psF64* outRow;
+        psC64* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F64[row];
+            inRow = in->data.C64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = cimag(inRow[col]);
+            }
+        }
+    } else {
+        psError(__func__,"Can not extract imaginary component from given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImageComplex(psImage *real, const psImage *imag)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+    psImage* out = NULL;
+
+
+    if (real == NULL || imag == NULL) {
+        return NULL;
+    }
+
+    type = real->type.type;
+    numCols = real->numCols;
+    numRows = real->numRows;
+
+    if (imag->type.type != type) {
+        psError(__func__,"The inputs to psImageComplex must be the same type.");
+        return NULL;
+    }
+
+    if (IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__,"The inputs to psImageComplex can not be complex.");
+        return NULL;
+    }
+
+    if (type != PS_TYPE_F32 && type != PS_TYPE_F64) {
+        psError(__func__,"The input type to psImageComplex must be a floating point.");
+        return NULL;
+    }
+
+    if (type == PS_TYPE_F32) {
+        psC32* outRow;
+        psF32* realRow;
+        psF32* imagRow;
+
+        out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.C32[row];
+            realRow = real->data.F32[row];
+            imagRow = imag->data.F32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = realRow[col] + I*imagRow[col];
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psC64* outRow;
+        psF64* realRow;
+        psF64* imagRow;
+
+        out = psImageAlloc(numCols,numRows,PS_TYPE_C64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.C64[row];
+            realRow = real->data.F64[row];
+            imagRow = imag->data.F64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = realRow[col] + I*imagRow[col];
+            }
+        }
+    } else {
+        psError(__func__,"Can not merge real and imaginary portions for given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImageConjugate(psImage *out, const psImage *in)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numCols = in->numCols;
+    numRows = in->numRows;
+
+    /* if not a complex number, this is logically just a image copy */
+    if (! IS_PSELEMTYPE_COMPLEX(type)) {
+        // Warn user, as this is probably not expected
+        psLogMsg(__func__,PS_LOG_WARN,"Complex Conjugate of a non-Complex type called for. "
+                 "Image copy was performed instead.");
+        return psImageCopy(out,in,type);
+    }
+
+    if (type == PS_TYPE_C32) {
+        psC32* outRow;
+        psC32* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_C32);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.C32[row];
+            inRow = in->data.C32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = crealf(inRow[col]) - I*cimagf(inRow[col]);
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psC64* outRow;
+        psC64* inRow;
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_C64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.C64[row];
+            inRow = in->data.C64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                outRow[col] = creal(inRow[col]) - I*cimag(inRow[col]);
+            }
+        }
+    } else {
+        psError(__func__,"Can not compute complex conjugate for given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImagePowerSpectrum(const psImage* in)
+{
+    psElemType type;
+    unsigned int numCols;
+    unsigned int numRows;
+    int numElementsSquared;
+    psImage* out = NULL;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numCols = in->numCols;
+    numRows = in->numRows;
+    numElementsSquared = numCols*numCols*numRows*numRows;
+
+    /* if not a complex number, this is not implemented */
+    if (! IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__,"Power Spectrum for non-complex inputs is not implemented.");
+        return NULL;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outRow;
+        psC32* inRow;
+        psF32 real;
+        psF32 imag;
+
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F32[row];
+            inRow = in->data.C32[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                real = crealf(inRow[col]);
+                imag = cimagf(inRow[col]);
+                outRow[col] = (real*real+imag*imag)/numElementsSquared;
+            }
+        }
+    } else if (type == PS_TYPE_C64) {
+        psF64* outRow;
+        psC64* inRow;
+        psF64 real;
+        psF64 imag;
+
+
+        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        for (unsigned int row=0;row<numRows;row++) {
+            outRow = out->data.F64[row];
+            inRow = in->data.C64[row];
+
+            for (unsigned int col=0;col<numCols;col++) {
+                real = crealf(inRow[col]);
+                imag = cimagf(inRow[col]);
+                outRow[col] = real*real+imag*imag/numElementsSquared;
+            }
+        }
+    } else {
+        psError(__func__,"Can not power spectrum for given image type (%d).",
+                type);
+        psImageFree(out);
+        return NULL;
+    }
+
+    return out;
+
+}
Index: /trunk/psLib/src/fft/psVectorFFT.h
===================================================================
--- /trunk/psLib/src/fft/psVectorFFT.h	(revision 801)
+++ /trunk/psLib/src/fft/psVectorFFT.h	(revision 801)
@@ -0,0 +1,43 @@
+/** @file  psFFT.h
+ *
+ *  @brief Contains FFT transforms functions
+ *
+ *  @author Robert DeSonia, MHPCC
+ *
+ *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-05-28 03:15:40 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ */
+
+#ifndef PS_FFT_H
+#define PS_FFT_H
+
+#include <fftw3.h>
+
+#include "psImage.h"
+#include "psVector.h"
+
+/** Details on FFT implementation (private). */
+
+
+typedef enum {
+    PS_FFT_FORWARD = 0,                 ///< psImageFFT/psVectorFFT should perform a forward FFT.
+    PS_FFT_REVERSE = 1,                 ///< psImageFFT/psVectorFFT should perform a reverse FFT.
+    PS_FFT_REAL_RESULT = 2
+                         /**< psImageFFT/psVectorFFT should return a real image.  This is valid for
+                          *   only reverse FFT, i.e., the psImageFFT/psVectorFFT flag parameter
+                          *   should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT.
+                          */
+} psFftFlags;
+
+psImage *psImageFFT(const psImage *image, int flags);
+psImage *psImageReal(psImage *out, const psImage* in);
+psImage *psImageImaginary(psImage *out, const psImage* in);
+psImage *psImageComplex(psImage *real, const psImage *imag);
+psImage *psImageConjugate(psImage *out, const psImage *in);
+psImage *psImagePowerSpectrum(const psImage* in);
+
+
+#endif
+