Index: /trunk/psLib/src/dataManip/psFFT.c
===================================================================
--- /trunk/psLib/src/dataManip/psFFT.c	(revision 813)
+++ /trunk/psLib/src/dataManip/psFFT.c	(revision 814)
@@ -5,6 +5,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-05-28 21:07:17 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-05-29 01:11:01 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -25,5 +25,5 @@
 static bool p_fftwWisdomImported = false;
 
-psImage* psImageFFT(const psImage *image, int flags)
+psImage* psImageFFT(psImage* out, const psImage* in, int flags)
 {
     static fftwf_plan lastForwardPlan = NULL;
@@ -39,8 +39,7 @@
     psElemType type;
     fftwf_plan plan;
-    psImage* out = NULL;
 
     /* got good image data? */
-    if (image==NULL) {
+    if (in==NULL) {
         return NULL;
     }
@@ -51,5 +50,5 @@
     }
 
-    type = image->type.type;
+    type = in->type.type;
 
     if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) {
@@ -72,6 +71,6 @@
     }
 
-    numRows = image->numRows;
-    numCols = (image->numCols-1)*2;
+    numRows = in->numRows;
+    numCols = in->numCols;
     // n.b., numCols needs adjustment for c->r, see below.
 
@@ -80,7 +79,7 @@
             // 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);
+            numCols = (numCols-1)*2;
+
+            out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
 
             if (lastReversePlanRows != numRows || lastReversePlanCols != numCols
@@ -88,15 +87,17 @@
                     || lastReversePlan == NULL) {
                 lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows,
-                                                        (fftwf_complex*)image->data.C32[0],out->data.F32[0],
+                                                        (fftwf_complex*)in->data.C32[0],
+                                                        out->data.F32[0],
                                                         P_FFTW_PLAN_RIGOR);
             }
         } else { // complex result desired
-            out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+            out = psImageRecycle(out,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,
+                                                    (fftwf_complex*)in->data.C32[0],
+                                                    (fftwf_complex*)out->data.C32[0],
+                                                    FFTW_BACKWARD,
                                                     P_FFTW_PLAN_RIGOR);
             }
@@ -106,5 +107,5 @@
         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);
+            out = psImageRecycle(out,numCols/2+1,numRows,PS_TYPE_C32);
 
             if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols
@@ -112,16 +113,18 @@
                 /* need new plan created */
                 lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows,
-                                                        image->data.F32[0],(fftwf_complex*)out->data.C32[0],
+                                                        in->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);
+            out = psImageRecycle(out,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*)in->data.C32[0],
                                                     (fftwf_complex*)out->data.C32[0],
-                                                    FFTW_FORWARD,P_FFTW_PLAN_RIGOR);
+                                                    FFTW_FORWARD,
+                                                    P_FFTW_PLAN_RIGOR);
             }
         }
@@ -171,5 +174,5 @@
         psC32* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F32[row];
@@ -184,5 +187,5 @@
         psC64* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F64[row];
@@ -223,5 +226,5 @@
         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);
+        out = psImageRecycle(out,numCols,numRows,type);
         memset(out->data.v[0],0,PSELEMTYPE_SIZEOF(type)*numCols*numRows);
         return out;
@@ -232,5 +235,5 @@
         psC32* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F32[row];
@@ -245,5 +248,5 @@
         psC64* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F64[row];
@@ -264,10 +267,9 @@
 }
 
-psImage *psImageComplex(psImage *real, const psImage *imag)
+psImage *psImageComplex(psImage* out, psImage *real, const psImage *imag)
 {
     psElemType type;
     unsigned int numCols;
     unsigned int numRows;
-    psImage* out = NULL;
 
 
@@ -282,4 +284,10 @@
     if (imag->type.type != type) {
         psError(__func__,"The inputs to psImageComplex must be the same type.");
+        return NULL;
+    }
+
+    if (imag->numCols != numCols ||
+            imag->numRows != numRows) {
+        psError(__func__,"The inputs to psImageComplex must be the same dimensions.");
         return NULL;
     }
@@ -300,5 +308,5 @@
         psF32* imagRow;
 
-        out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32);
 
         for (unsigned int row=0;row<numRows;row++) {
@@ -311,10 +319,10 @@
             }
         }
-    } else if (type == PS_TYPE_C64) {
+    } else if (type == PS_TYPE_F64) {
         psC64* outRow;
         psF64* realRow;
         psF64* imagRow;
 
-        out = psImageAlloc(numCols,numRows,PS_TYPE_C64);
+        out = psImageRecycle(out, numCols,numRows,PS_TYPE_C64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.C64[row];
@@ -363,5 +371,5 @@
         psC32* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_C32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.C32[row];
@@ -376,5 +384,5 @@
         psC64* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_C64);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_C64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.C64[row];
@@ -395,5 +403,5 @@
 }
 
-psImage *psImagePowerSpectrum(const psImage* in)
+psImage *psImagePowerSpectrum(psImage* out, const psImage* in)
 {
     psElemType type;
@@ -401,6 +409,4 @@
     unsigned int numRows;
     int numElementsSquared;
-    psImage* out = NULL;
-
 
     if (in == NULL) {
@@ -426,5 +432,5 @@
 
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F32[row];
@@ -444,5 +450,5 @@
 
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F64[row];
@@ -465,2 +471,347 @@
 
 }
+
+/************************************** Vector Functions ***************************************/
+
+psVector* psVectorFFT(psVector* out, const psVector* in, int flags)
+{
+    static fftwf_plan lastForwardPlan = NULL;
+    static int lastForwardPlanType = 0;
+    static int lastForwardPlanElements = 0;
+    static fftwf_plan lastReversePlan = NULL;
+    static int lastReversePlanResultType = 0;
+    static int lastReversePlanElements = 0;
+    unsigned int numElements;
+    psElemType type;
+    fftwf_plan plan;
+
+    /* got good image data? */
+    if (in==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 = in->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) {
+        fftwf_import_system_wisdom();
+        p_fftwWisdomImported = true;
+    }
+
+    numElements = in->n;
+    // n.b., numElements 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)
+            numElements = (numElements-1)*2;
+
+            out = psVectorRecycle(out,PS_TYPE_F32,numElements);
+
+            if (lastReversePlanElements != numElements ||
+                    lastReversePlanResultType != PS_TYPE_F32 ||
+                    lastReversePlan == NULL) {
+                lastReversePlan = fftwf_plan_dft_c2r_1d(numElements,
+                                                        (fftwf_complex*)in->vec.cf,
+                                                        out->vec.f,
+                                                        P_FFTW_PLAN_RIGOR);
+            }
+        } else { // complex result desired
+            out = psVectorRecycle(out,PS_TYPE_C32,numElements);
+            if (lastReversePlanElements != numElements ||
+                    lastReversePlanResultType != PS_TYPE_C32 ||
+                    lastReversePlan == NULL) {
+                lastReversePlan = fftwf_plan_dft_1d(numElements,
+                                                    (fftwf_complex*)in->vec.cf,
+                                                    (fftwf_complex*)out->vec.cf,
+                                                    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 = psVectorRecycle(out,PS_TYPE_C32,numElements);
+
+            if (lastForwardPlanElements != numElements ||
+                    lastForwardPlanType != type ||
+                    lastForwardPlan == NULL)  {
+                /* need new plan created */
+                lastForwardPlan = fftwf_plan_dft_r2c_1d(numElements,
+                                                        in->vec.f,
+                                                        (fftwf_complex*)out->vec.cf,
+                                                        P_FFTW_PLAN_RIGOR);
+            }
+        } else if (type == PS_TYPE_C32) { // complex data input
+            out = psVectorRecycle(out,PS_TYPE_C32,numElements);
+
+            if (lastForwardPlanElements != numElements ||
+                    lastForwardPlanType != type ||
+                    lastForwardPlan == NULL)  {
+                lastForwardPlan = fftwf_plan_dft_1d(numElements,
+                                                    (fftwf_complex*)in->vec.cf,
+                                                    (fftwf_complex*)out->vec.cf,
+                                                    FFTW_FORWARD,
+                                                    P_FFTW_PLAN_RIGOR);
+            }
+        }
+        plan = lastForwardPlan;
+    }
+
+    /* check if a plan exists now*/
+    if (plan == NULL) {
+        psError(__func__,"Failed to create FFTW plan.");
+        psVectorFree(out);
+        return NULL;
+    }
+
+    /* finally, call FFTW with the plan made above */
+    fftwf_execute(plan);
+
+    return out;
+
+}
+
+
+psVector *psVectorReal(psVector *out, const psVector* in)
+{
+    psElemType type;
+    unsigned int numElements;
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numElements = in->n;
+
+    /* if not a complex number, this is logically just a copy */
+    if (! PS_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 a vector copy was performed.");
+        out = psVectorRecycle(out,type,numElements);
+        memcpy(out->vec.v,in->vec.v,numElements*PSELEMTYPE_SIZEOF(type));
+        return out;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outVec;
+        psC32* inVec = in->vec.cf;
+
+        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
+        outVec = out->vec.f;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            outVec[i] = crealf(inVec[i]);
+        }
+    } else {
+        psError(__func__,"Can not extract real component from given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorImaginary(psVector *out, const psVector* in)
+{
+    psElemType type;
+    unsigned int numElements;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numElements = in->n;
+
+    /* if not a complex number, this is logically just zeroed image of same size */
+    if (! PS_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 zeroed vector was returned.");
+        out = psVectorRecycle(out,numElements,type);
+        memset(out->vec.v,0,PSELEMTYPE_SIZEOF(type)*numElements);
+        return out;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outVec;
+        psC32* inVec = in->vec.cf;
+
+        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
+        outVec = out->vec.f;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            outVec[i] = cimagf(inVec[i]);
+        }
+    } else {
+        psError(__func__,"Can not extract imaginary component from given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorComplex(psVector* out, psVector *real, const psVector *imag)
+{
+    psElemType type;
+    unsigned int numElements;
+
+
+    if (real == NULL || imag == NULL) {
+        return NULL;
+    }
+
+    type = real->type.type;
+    if (real->n >= imag->n) {
+        numElements = real->n;
+    } else {
+        numElements = imag->n;
+    }
+
+    if (imag->type.type != type) {
+        psError(__func__,"The inputs to psVectorComplex must be the same type.");
+        return NULL;
+    }
+
+    if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__,"The inputs to psVectorComplex can not be complex.");
+        return NULL;
+    }
+
+    if (type == PS_TYPE_F32) {
+        psC32* outVec;
+        psF32* realVec = real->vec.f;
+        psF32* imagVec = imag->vec.f;
+
+        out = psVectorRecycle(out,PS_TYPE_C32,numElements);
+        outVec = out->vec.cf;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            outVec[i] = realVec[i] + I*imagVec[i];
+        }
+    } else {
+        psError(__func__,"Can not merge real and imaginary portions for given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorConjugate(psVector *out, const psVector *in)
+{
+    psElemType type;
+    unsigned int numElements;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numElements = in->n;
+
+    /* if not a complex number, this is logically just a image copy */
+    if (! PS_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. "
+                 "Vector copy was performed instead.");
+
+        out = psVectorRecycle(out,type,numElements);
+        memcpy(out->vec.v,in->vec.v,PSELEMTYPE_SIZEOF(type)*numElements);
+        return out;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psC32* outVec;
+        psC32* inVec = in->vec.cf;
+
+        out = psVectorRecycle(out,PS_TYPE_C32,numElements);
+        outVec = out->vec.cf;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            outVec[i] = crealf(inVec[i]) - I*cimagf(inVec[i]);
+        }
+    } else {
+        psError(__func__,"Can not compute complex conjugate for given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorPowerSpectrum(psVector* out, const psVector* in)
+{
+    psElemType type;
+    unsigned int numElements;
+    int numElementsSquared;
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numElements = in->n;
+    numElementsSquared = numElements*numElements;
+
+    /* if not a complex number, this is not implemented */
+    if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__,"Power Spectrum for non-complex inputs is not implemented.");
+        return NULL;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outVec;
+        psC32* inVec = in->vec.cf;
+        psF32 real;
+        psF32 imag;
+
+
+        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
+        outVec = out->vec.f;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            real = crealf(inVec[i]);
+            imag = cimagf(inVec[i]);
+            outVec[i] = (real*real+imag*imag)/numElementsSquared;
+        }
+    } else {
+        psError(__func__,"Can not power spectrum for given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+
+}
Index: /trunk/psLib/src/dataManip/psFFT.h
===================================================================
--- /trunk/psLib/src/dataManip/psFFT.h	(revision 813)
+++ /trunk/psLib/src/dataManip/psFFT.h	(revision 814)
@@ -5,6 +5,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-05-28 03:15:40 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-05-29 01:11:01 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23,19 +23,30 @@
 
 typedef enum {
-    PS_FFT_FORWARD = 0,                 ///< psImageFFT/psVectorFFT should perform a forward FFT.
-    PS_FFT_REVERSE = 1,                 ///< psImageFFT/psVectorFFT should perform a reverse FFT.
+    /// psImageFFT/psVectorFFT should perform a forward FFT.
+    PS_FFT_FORWARD = 0,
+
+    ///< psImageFFT/psVectorFFT should perform a reverse FFT.
+    PS_FFT_REVERSE = 1,
+
+    /** 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.
+    */
     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);
+psImage* psImageFFT(psImage* out, const psImage* in, int flags);
+psImage* psImageReal(psImage *out, const psImage* in);
+psImage* psImageImaginary(psImage *out, const psImage* in);
+psImage* psImageComplex(psImage* out, psImage *real, const psImage *imag);
+psImage* psImageConjugate(psImage *out, const psImage *in);
+psImage* psImagePowerSpectrum(psImage* out, const psImage* in);
+
+psVector* psVectorFFT(psVector* out, const psVector* in, int flags);
+psVector* psVectorReal(psVector* out, const psVector* in);
+psVector* psVectorImaginary(psVector* out, const psVector* in);
+psVector* psVectorComplex(psVector* out, psVector* real, const psVector* imag);
+psVector* psVectorConjugate(psVector* out, const psVector* in);
+psVector* psVectorPowerSpectrum(psVector* out, const psVector* in);
 
 
Index: /trunk/psLib/src/dataManip/psVectorFFT.c
===================================================================
--- /trunk/psLib/src/dataManip/psVectorFFT.c	(revision 813)
+++ /trunk/psLib/src/dataManip/psVectorFFT.c	(revision 814)
@@ -5,6 +5,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-05-28 21:07:17 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-05-29 01:11:01 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -25,5 +25,5 @@
 static bool p_fftwWisdomImported = false;
 
-psImage* psImageFFT(const psImage *image, int flags)
+psImage* psImageFFT(psImage* out, const psImage* in, int flags)
 {
     static fftwf_plan lastForwardPlan = NULL;
@@ -39,8 +39,7 @@
     psElemType type;
     fftwf_plan plan;
-    psImage* out = NULL;
 
     /* got good image data? */
-    if (image==NULL) {
+    if (in==NULL) {
         return NULL;
     }
@@ -51,5 +50,5 @@
     }
 
-    type = image->type.type;
+    type = in->type.type;
 
     if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) {
@@ -72,6 +71,6 @@
     }
 
-    numRows = image->numRows;
-    numCols = (image->numCols-1)*2;
+    numRows = in->numRows;
+    numCols = in->numCols;
     // n.b., numCols needs adjustment for c->r, see below.
 
@@ -80,7 +79,7 @@
             // 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);
+            numCols = (numCols-1)*2;
+
+            out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
 
             if (lastReversePlanRows != numRows || lastReversePlanCols != numCols
@@ -88,15 +87,17 @@
                     || lastReversePlan == NULL) {
                 lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows,
-                                                        (fftwf_complex*)image->data.C32[0],out->data.F32[0],
+                                                        (fftwf_complex*)in->data.C32[0],
+                                                        out->data.F32[0],
                                                         P_FFTW_PLAN_RIGOR);
             }
         } else { // complex result desired
-            out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+            out = psImageRecycle(out,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,
+                                                    (fftwf_complex*)in->data.C32[0],
+                                                    (fftwf_complex*)out->data.C32[0],
+                                                    FFTW_BACKWARD,
                                                     P_FFTW_PLAN_RIGOR);
             }
@@ -106,5 +107,5 @@
         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);
+            out = psImageRecycle(out,numCols/2+1,numRows,PS_TYPE_C32);
 
             if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols
@@ -112,16 +113,18 @@
                 /* need new plan created */
                 lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows,
-                                                        image->data.F32[0],(fftwf_complex*)out->data.C32[0],
+                                                        in->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);
+            out = psImageRecycle(out,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*)in->data.C32[0],
                                                     (fftwf_complex*)out->data.C32[0],
-                                                    FFTW_FORWARD,P_FFTW_PLAN_RIGOR);
+                                                    FFTW_FORWARD,
+                                                    P_FFTW_PLAN_RIGOR);
             }
         }
@@ -171,5 +174,5 @@
         psC32* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F32[row];
@@ -184,5 +187,5 @@
         psC64* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F64[row];
@@ -223,5 +226,5 @@
         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);
+        out = psImageRecycle(out,numCols,numRows,type);
         memset(out->data.v[0],0,PSELEMTYPE_SIZEOF(type)*numCols*numRows);
         return out;
@@ -232,5 +235,5 @@
         psC32* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F32[row];
@@ -245,5 +248,5 @@
         psC64* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F64[row];
@@ -264,10 +267,9 @@
 }
 
-psImage *psImageComplex(psImage *real, const psImage *imag)
+psImage *psImageComplex(psImage* out, psImage *real, const psImage *imag)
 {
     psElemType type;
     unsigned int numCols;
     unsigned int numRows;
-    psImage* out = NULL;
 
 
@@ -282,4 +284,10 @@
     if (imag->type.type != type) {
         psError(__func__,"The inputs to psImageComplex must be the same type.");
+        return NULL;
+    }
+
+    if (imag->numCols != numCols ||
+            imag->numRows != numRows) {
+        psError(__func__,"The inputs to psImageComplex must be the same dimensions.");
         return NULL;
     }
@@ -300,5 +308,5 @@
         psF32* imagRow;
 
-        out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32);
 
         for (unsigned int row=0;row<numRows;row++) {
@@ -311,10 +319,10 @@
             }
         }
-    } else if (type == PS_TYPE_C64) {
+    } else if (type == PS_TYPE_F64) {
         psC64* outRow;
         psF64* realRow;
         psF64* imagRow;
 
-        out = psImageAlloc(numCols,numRows,PS_TYPE_C64);
+        out = psImageRecycle(out, numCols,numRows,PS_TYPE_C64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.C64[row];
@@ -363,5 +371,5 @@
         psC32* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_C32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.C32[row];
@@ -376,5 +384,5 @@
         psC64* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_C64);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_C64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.C64[row];
@@ -395,5 +403,5 @@
 }
 
-psImage *psImagePowerSpectrum(const psImage* in)
+psImage *psImagePowerSpectrum(psImage* out, const psImage* in)
 {
     psElemType type;
@@ -401,6 +409,4 @@
     unsigned int numRows;
     int numElementsSquared;
-    psImage* out = NULL;
-
 
     if (in == NULL) {
@@ -426,5 +432,5 @@
 
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F32[row];
@@ -444,5 +450,5 @@
 
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F64[row];
@@ -465,2 +471,347 @@
 
 }
+
+/************************************** Vector Functions ***************************************/
+
+psVector* psVectorFFT(psVector* out, const psVector* in, int flags)
+{
+    static fftwf_plan lastForwardPlan = NULL;
+    static int lastForwardPlanType = 0;
+    static int lastForwardPlanElements = 0;
+    static fftwf_plan lastReversePlan = NULL;
+    static int lastReversePlanResultType = 0;
+    static int lastReversePlanElements = 0;
+    unsigned int numElements;
+    psElemType type;
+    fftwf_plan plan;
+
+    /* got good image data? */
+    if (in==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 = in->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) {
+        fftwf_import_system_wisdom();
+        p_fftwWisdomImported = true;
+    }
+
+    numElements = in->n;
+    // n.b., numElements 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)
+            numElements = (numElements-1)*2;
+
+            out = psVectorRecycle(out,PS_TYPE_F32,numElements);
+
+            if (lastReversePlanElements != numElements ||
+                    lastReversePlanResultType != PS_TYPE_F32 ||
+                    lastReversePlan == NULL) {
+                lastReversePlan = fftwf_plan_dft_c2r_1d(numElements,
+                                                        (fftwf_complex*)in->vec.cf,
+                                                        out->vec.f,
+                                                        P_FFTW_PLAN_RIGOR);
+            }
+        } else { // complex result desired
+            out = psVectorRecycle(out,PS_TYPE_C32,numElements);
+            if (lastReversePlanElements != numElements ||
+                    lastReversePlanResultType != PS_TYPE_C32 ||
+                    lastReversePlan == NULL) {
+                lastReversePlan = fftwf_plan_dft_1d(numElements,
+                                                    (fftwf_complex*)in->vec.cf,
+                                                    (fftwf_complex*)out->vec.cf,
+                                                    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 = psVectorRecycle(out,PS_TYPE_C32,numElements);
+
+            if (lastForwardPlanElements != numElements ||
+                    lastForwardPlanType != type ||
+                    lastForwardPlan == NULL)  {
+                /* need new plan created */
+                lastForwardPlan = fftwf_plan_dft_r2c_1d(numElements,
+                                                        in->vec.f,
+                                                        (fftwf_complex*)out->vec.cf,
+                                                        P_FFTW_PLAN_RIGOR);
+            }
+        } else if (type == PS_TYPE_C32) { // complex data input
+            out = psVectorRecycle(out,PS_TYPE_C32,numElements);
+
+            if (lastForwardPlanElements != numElements ||
+                    lastForwardPlanType != type ||
+                    lastForwardPlan == NULL)  {
+                lastForwardPlan = fftwf_plan_dft_1d(numElements,
+                                                    (fftwf_complex*)in->vec.cf,
+                                                    (fftwf_complex*)out->vec.cf,
+                                                    FFTW_FORWARD,
+                                                    P_FFTW_PLAN_RIGOR);
+            }
+        }
+        plan = lastForwardPlan;
+    }
+
+    /* check if a plan exists now*/
+    if (plan == NULL) {
+        psError(__func__,"Failed to create FFTW plan.");
+        psVectorFree(out);
+        return NULL;
+    }
+
+    /* finally, call FFTW with the plan made above */
+    fftwf_execute(plan);
+
+    return out;
+
+}
+
+
+psVector *psVectorReal(psVector *out, const psVector* in)
+{
+    psElemType type;
+    unsigned int numElements;
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numElements = in->n;
+
+    /* if not a complex number, this is logically just a copy */
+    if (! PS_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 a vector copy was performed.");
+        out = psVectorRecycle(out,type,numElements);
+        memcpy(out->vec.v,in->vec.v,numElements*PSELEMTYPE_SIZEOF(type));
+        return out;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outVec;
+        psC32* inVec = in->vec.cf;
+
+        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
+        outVec = out->vec.f;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            outVec[i] = crealf(inVec[i]);
+        }
+    } else {
+        psError(__func__,"Can not extract real component from given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorImaginary(psVector *out, const psVector* in)
+{
+    psElemType type;
+    unsigned int numElements;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numElements = in->n;
+
+    /* if not a complex number, this is logically just zeroed image of same size */
+    if (! PS_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 zeroed vector was returned.");
+        out = psVectorRecycle(out,numElements,type);
+        memset(out->vec.v,0,PSELEMTYPE_SIZEOF(type)*numElements);
+        return out;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outVec;
+        psC32* inVec = in->vec.cf;
+
+        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
+        outVec = out->vec.f;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            outVec[i] = cimagf(inVec[i]);
+        }
+    } else {
+        psError(__func__,"Can not extract imaginary component from given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorComplex(psVector* out, psVector *real, const psVector *imag)
+{
+    psElemType type;
+    unsigned int numElements;
+
+
+    if (real == NULL || imag == NULL) {
+        return NULL;
+    }
+
+    type = real->type.type;
+    if (real->n >= imag->n) {
+        numElements = real->n;
+    } else {
+        numElements = imag->n;
+    }
+
+    if (imag->type.type != type) {
+        psError(__func__,"The inputs to psVectorComplex must be the same type.");
+        return NULL;
+    }
+
+    if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__,"The inputs to psVectorComplex can not be complex.");
+        return NULL;
+    }
+
+    if (type == PS_TYPE_F32) {
+        psC32* outVec;
+        psF32* realVec = real->vec.f;
+        psF32* imagVec = imag->vec.f;
+
+        out = psVectorRecycle(out,PS_TYPE_C32,numElements);
+        outVec = out->vec.cf;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            outVec[i] = realVec[i] + I*imagVec[i];
+        }
+    } else {
+        psError(__func__,"Can not merge real and imaginary portions for given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorConjugate(psVector *out, const psVector *in)
+{
+    psElemType type;
+    unsigned int numElements;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numElements = in->n;
+
+    /* if not a complex number, this is logically just a image copy */
+    if (! PS_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. "
+                 "Vector copy was performed instead.");
+
+        out = psVectorRecycle(out,type,numElements);
+        memcpy(out->vec.v,in->vec.v,PSELEMTYPE_SIZEOF(type)*numElements);
+        return out;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psC32* outVec;
+        psC32* inVec = in->vec.cf;
+
+        out = psVectorRecycle(out,PS_TYPE_C32,numElements);
+        outVec = out->vec.cf;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            outVec[i] = crealf(inVec[i]) - I*cimagf(inVec[i]);
+        }
+    } else {
+        psError(__func__,"Can not compute complex conjugate for given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorPowerSpectrum(psVector* out, const psVector* in)
+{
+    psElemType type;
+    unsigned int numElements;
+    int numElementsSquared;
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numElements = in->n;
+    numElementsSquared = numElements*numElements;
+
+    /* if not a complex number, this is not implemented */
+    if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__,"Power Spectrum for non-complex inputs is not implemented.");
+        return NULL;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outVec;
+        psC32* inVec = in->vec.cf;
+        psF32 real;
+        psF32 imag;
+
+
+        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
+        outVec = out->vec.f;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            real = crealf(inVec[i]);
+            imag = cimagf(inVec[i]);
+            outVec[i] = (real*real+imag*imag)/numElementsSquared;
+        }
+    } else {
+        psError(__func__,"Can not power spectrum for given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+
+}
Index: /trunk/psLib/src/dataManip/psVectorFFT.h
===================================================================
--- /trunk/psLib/src/dataManip/psVectorFFT.h	(revision 813)
+++ /trunk/psLib/src/dataManip/psVectorFFT.h	(revision 814)
@@ -5,6 +5,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-05-28 03:15:40 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-05-29 01:11:01 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23,19 +23,30 @@
 
 typedef enum {
-    PS_FFT_FORWARD = 0,                 ///< psImageFFT/psVectorFFT should perform a forward FFT.
-    PS_FFT_REVERSE = 1,                 ///< psImageFFT/psVectorFFT should perform a reverse FFT.
+    /// psImageFFT/psVectorFFT should perform a forward FFT.
+    PS_FFT_FORWARD = 0,
+
+    ///< psImageFFT/psVectorFFT should perform a reverse FFT.
+    PS_FFT_REVERSE = 1,
+
+    /** 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.
+    */
     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);
+psImage* psImageFFT(psImage* out, const psImage* in, int flags);
+psImage* psImageReal(psImage *out, const psImage* in);
+psImage* psImageImaginary(psImage *out, const psImage* in);
+psImage* psImageComplex(psImage* out, psImage *real, const psImage *imag);
+psImage* psImageConjugate(psImage *out, const psImage *in);
+psImage* psImagePowerSpectrum(psImage* out, const psImage* in);
+
+psVector* psVectorFFT(psVector* out, const psVector* in, int flags);
+psVector* psVectorReal(psVector* out, const psVector* in);
+psVector* psVectorImaginary(psVector* out, const psVector* in);
+psVector* psVectorComplex(psVector* out, psVector* real, const psVector* imag);
+psVector* psVectorConjugate(psVector* out, const psVector* in);
+psVector* psVectorPowerSpectrum(psVector* out, const psVector* in);
 
 
Index: /trunk/psLib/src/fft/psVectorFFT.c
===================================================================
--- /trunk/psLib/src/fft/psVectorFFT.c	(revision 813)
+++ /trunk/psLib/src/fft/psVectorFFT.c	(revision 814)
@@ -5,6 +5,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-05-28 21:07:17 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-05-29 01:11:01 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -25,5 +25,5 @@
 static bool p_fftwWisdomImported = false;
 
-psImage* psImageFFT(const psImage *image, int flags)
+psImage* psImageFFT(psImage* out, const psImage* in, int flags)
 {
     static fftwf_plan lastForwardPlan = NULL;
@@ -39,8 +39,7 @@
     psElemType type;
     fftwf_plan plan;
-    psImage* out = NULL;
 
     /* got good image data? */
-    if (image==NULL) {
+    if (in==NULL) {
         return NULL;
     }
@@ -51,5 +50,5 @@
     }
 
-    type = image->type.type;
+    type = in->type.type;
 
     if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) {
@@ -72,6 +71,6 @@
     }
 
-    numRows = image->numRows;
-    numCols = (image->numCols-1)*2;
+    numRows = in->numRows;
+    numCols = in->numCols;
     // n.b., numCols needs adjustment for c->r, see below.
 
@@ -80,7 +79,7 @@
             // 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);
+            numCols = (numCols-1)*2;
+
+            out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
 
             if (lastReversePlanRows != numRows || lastReversePlanCols != numCols
@@ -88,15 +87,17 @@
                     || lastReversePlan == NULL) {
                 lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows,
-                                                        (fftwf_complex*)image->data.C32[0],out->data.F32[0],
+                                                        (fftwf_complex*)in->data.C32[0],
+                                                        out->data.F32[0],
                                                         P_FFTW_PLAN_RIGOR);
             }
         } else { // complex result desired
-            out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+            out = psImageRecycle(out,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,
+                                                    (fftwf_complex*)in->data.C32[0],
+                                                    (fftwf_complex*)out->data.C32[0],
+                                                    FFTW_BACKWARD,
                                                     P_FFTW_PLAN_RIGOR);
             }
@@ -106,5 +107,5 @@
         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);
+            out = psImageRecycle(out,numCols/2+1,numRows,PS_TYPE_C32);
 
             if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols
@@ -112,16 +113,18 @@
                 /* need new plan created */
                 lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows,
-                                                        image->data.F32[0],(fftwf_complex*)out->data.C32[0],
+                                                        in->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);
+            out = psImageRecycle(out,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*)in->data.C32[0],
                                                     (fftwf_complex*)out->data.C32[0],
-                                                    FFTW_FORWARD,P_FFTW_PLAN_RIGOR);
+                                                    FFTW_FORWARD,
+                                                    P_FFTW_PLAN_RIGOR);
             }
         }
@@ -171,5 +174,5 @@
         psC32* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F32[row];
@@ -184,5 +187,5 @@
         psC64* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F64[row];
@@ -223,5 +226,5 @@
         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);
+        out = psImageRecycle(out,numCols,numRows,type);
         memset(out->data.v[0],0,PSELEMTYPE_SIZEOF(type)*numCols*numRows);
         return out;
@@ -232,5 +235,5 @@
         psC32* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F32[row];
@@ -245,5 +248,5 @@
         psC64* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F64[row];
@@ -264,10 +267,9 @@
 }
 
-psImage *psImageComplex(psImage *real, const psImage *imag)
+psImage *psImageComplex(psImage* out, psImage *real, const psImage *imag)
 {
     psElemType type;
     unsigned int numCols;
     unsigned int numRows;
-    psImage* out = NULL;
 
 
@@ -282,4 +284,10 @@
     if (imag->type.type != type) {
         psError(__func__,"The inputs to psImageComplex must be the same type.");
+        return NULL;
+    }
+
+    if (imag->numCols != numCols ||
+            imag->numRows != numRows) {
+        psError(__func__,"The inputs to psImageComplex must be the same dimensions.");
         return NULL;
     }
@@ -300,5 +308,5 @@
         psF32* imagRow;
 
-        out = psImageAlloc(numCols,numRows,PS_TYPE_C32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32);
 
         for (unsigned int row=0;row<numRows;row++) {
@@ -311,10 +319,10 @@
             }
         }
-    } else if (type == PS_TYPE_C64) {
+    } else if (type == PS_TYPE_F64) {
         psC64* outRow;
         psF64* realRow;
         psF64* imagRow;
 
-        out = psImageAlloc(numCols,numRows,PS_TYPE_C64);
+        out = psImageRecycle(out, numCols,numRows,PS_TYPE_C64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.C64[row];
@@ -363,5 +371,5 @@
         psC32* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_C32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.C32[row];
@@ -376,5 +384,5 @@
         psC64* inRow;
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_C64);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_C64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.C64[row];
@@ -395,5 +403,5 @@
 }
 
-psImage *psImagePowerSpectrum(const psImage* in)
+psImage *psImagePowerSpectrum(psImage* out, const psImage* in)
 {
     psElemType type;
@@ -401,6 +409,4 @@
     unsigned int numRows;
     int numElementsSquared;
-    psImage* out = NULL;
-
 
     if (in == NULL) {
@@ -426,5 +432,5 @@
 
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F32);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F32[row];
@@ -444,5 +450,5 @@
 
 
-        out = psImageRealloc(out,numCols,numRows,PS_TYPE_F64);
+        out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64);
         for (unsigned int row=0;row<numRows;row++) {
             outRow = out->data.F64[row];
@@ -465,2 +471,347 @@
 
 }
+
+/************************************** Vector Functions ***************************************/
+
+psVector* psVectorFFT(psVector* out, const psVector* in, int flags)
+{
+    static fftwf_plan lastForwardPlan = NULL;
+    static int lastForwardPlanType = 0;
+    static int lastForwardPlanElements = 0;
+    static fftwf_plan lastReversePlan = NULL;
+    static int lastReversePlanResultType = 0;
+    static int lastReversePlanElements = 0;
+    unsigned int numElements;
+    psElemType type;
+    fftwf_plan plan;
+
+    /* got good image data? */
+    if (in==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 = in->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) {
+        fftwf_import_system_wisdom();
+        p_fftwWisdomImported = true;
+    }
+
+    numElements = in->n;
+    // n.b., numElements 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)
+            numElements = (numElements-1)*2;
+
+            out = psVectorRecycle(out,PS_TYPE_F32,numElements);
+
+            if (lastReversePlanElements != numElements ||
+                    lastReversePlanResultType != PS_TYPE_F32 ||
+                    lastReversePlan == NULL) {
+                lastReversePlan = fftwf_plan_dft_c2r_1d(numElements,
+                                                        (fftwf_complex*)in->vec.cf,
+                                                        out->vec.f,
+                                                        P_FFTW_PLAN_RIGOR);
+            }
+        } else { // complex result desired
+            out = psVectorRecycle(out,PS_TYPE_C32,numElements);
+            if (lastReversePlanElements != numElements ||
+                    lastReversePlanResultType != PS_TYPE_C32 ||
+                    lastReversePlan == NULL) {
+                lastReversePlan = fftwf_plan_dft_1d(numElements,
+                                                    (fftwf_complex*)in->vec.cf,
+                                                    (fftwf_complex*)out->vec.cf,
+                                                    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 = psVectorRecycle(out,PS_TYPE_C32,numElements);
+
+            if (lastForwardPlanElements != numElements ||
+                    lastForwardPlanType != type ||
+                    lastForwardPlan == NULL)  {
+                /* need new plan created */
+                lastForwardPlan = fftwf_plan_dft_r2c_1d(numElements,
+                                                        in->vec.f,
+                                                        (fftwf_complex*)out->vec.cf,
+                                                        P_FFTW_PLAN_RIGOR);
+            }
+        } else if (type == PS_TYPE_C32) { // complex data input
+            out = psVectorRecycle(out,PS_TYPE_C32,numElements);
+
+            if (lastForwardPlanElements != numElements ||
+                    lastForwardPlanType != type ||
+                    lastForwardPlan == NULL)  {
+                lastForwardPlan = fftwf_plan_dft_1d(numElements,
+                                                    (fftwf_complex*)in->vec.cf,
+                                                    (fftwf_complex*)out->vec.cf,
+                                                    FFTW_FORWARD,
+                                                    P_FFTW_PLAN_RIGOR);
+            }
+        }
+        plan = lastForwardPlan;
+    }
+
+    /* check if a plan exists now*/
+    if (plan == NULL) {
+        psError(__func__,"Failed to create FFTW plan.");
+        psVectorFree(out);
+        return NULL;
+    }
+
+    /* finally, call FFTW with the plan made above */
+    fftwf_execute(plan);
+
+    return out;
+
+}
+
+
+psVector *psVectorReal(psVector *out, const psVector* in)
+{
+    psElemType type;
+    unsigned int numElements;
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numElements = in->n;
+
+    /* if not a complex number, this is logically just a copy */
+    if (! PS_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 a vector copy was performed.");
+        out = psVectorRecycle(out,type,numElements);
+        memcpy(out->vec.v,in->vec.v,numElements*PSELEMTYPE_SIZEOF(type));
+        return out;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outVec;
+        psC32* inVec = in->vec.cf;
+
+        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
+        outVec = out->vec.f;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            outVec[i] = crealf(inVec[i]);
+        }
+    } else {
+        psError(__func__,"Can not extract real component from given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorImaginary(psVector *out, const psVector* in)
+{
+    psElemType type;
+    unsigned int numElements;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numElements = in->n;
+
+    /* if not a complex number, this is logically just zeroed image of same size */
+    if (! PS_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 zeroed vector was returned.");
+        out = psVectorRecycle(out,numElements,type);
+        memset(out->vec.v,0,PSELEMTYPE_SIZEOF(type)*numElements);
+        return out;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outVec;
+        psC32* inVec = in->vec.cf;
+
+        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
+        outVec = out->vec.f;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            outVec[i] = cimagf(inVec[i]);
+        }
+    } else {
+        psError(__func__,"Can not extract imaginary component from given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorComplex(psVector* out, psVector *real, const psVector *imag)
+{
+    psElemType type;
+    unsigned int numElements;
+
+
+    if (real == NULL || imag == NULL) {
+        return NULL;
+    }
+
+    type = real->type.type;
+    if (real->n >= imag->n) {
+        numElements = real->n;
+    } else {
+        numElements = imag->n;
+    }
+
+    if (imag->type.type != type) {
+        psError(__func__,"The inputs to psVectorComplex must be the same type.");
+        return NULL;
+    }
+
+    if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__,"The inputs to psVectorComplex can not be complex.");
+        return NULL;
+    }
+
+    if (type == PS_TYPE_F32) {
+        psC32* outVec;
+        psF32* realVec = real->vec.f;
+        psF32* imagVec = imag->vec.f;
+
+        out = psVectorRecycle(out,PS_TYPE_C32,numElements);
+        outVec = out->vec.cf;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            outVec[i] = realVec[i] + I*imagVec[i];
+        }
+    } else {
+        psError(__func__,"Can not merge real and imaginary portions for given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorConjugate(psVector *out, const psVector *in)
+{
+    psElemType type;
+    unsigned int numElements;
+
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numElements = in->n;
+
+    /* if not a complex number, this is logically just a image copy */
+    if (! PS_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. "
+                 "Vector copy was performed instead.");
+
+        out = psVectorRecycle(out,type,numElements);
+        memcpy(out->vec.v,in->vec.v,PSELEMTYPE_SIZEOF(type)*numElements);
+        return out;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psC32* outVec;
+        psC32* inVec = in->vec.cf;
+
+        out = psVectorRecycle(out,PS_TYPE_C32,numElements);
+        outVec = out->vec.cf;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            outVec[i] = crealf(inVec[i]) - I*cimagf(inVec[i]);
+        }
+    } else {
+        psError(__func__,"Can not compute complex conjugate for given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorPowerSpectrum(psVector* out, const psVector* in)
+{
+    psElemType type;
+    unsigned int numElements;
+    int numElementsSquared;
+
+    if (in == NULL) {
+        return NULL;
+    }
+
+    type = in->type.type;
+    numElements = in->n;
+    numElementsSquared = numElements*numElements;
+
+    /* if not a complex number, this is not implemented */
+    if (! PS_IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__,"Power Spectrum for non-complex inputs is not implemented.");
+        return NULL;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32* outVec;
+        psC32* inVec = in->vec.cf;
+        psF32 real;
+        psF32 imag;
+
+
+        out = psVectorRecycle(out,PS_TYPE_F32,numElements);
+        outVec = out->vec.f;
+
+        for (unsigned int i=0;i<numElements;i++) {
+            real = crealf(inVec[i]);
+            imag = cimagf(inVec[i]);
+            outVec[i] = (real*real+imag*imag)/numElementsSquared;
+        }
+    } else {
+        psError(__func__,"Can not power spectrum for given vector type (%d).",
+                type);
+        psVectorFree(out);
+        return NULL;
+    }
+
+    return out;
+
+}
Index: /trunk/psLib/src/fft/psVectorFFT.h
===================================================================
--- /trunk/psLib/src/fft/psVectorFFT.h	(revision 813)
+++ /trunk/psLib/src/fft/psVectorFFT.h	(revision 814)
@@ -5,6 +5,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-05-28 03:15:40 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-05-29 01:11:01 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23,19 +23,30 @@
 
 typedef enum {
-    PS_FFT_FORWARD = 0,                 ///< psImageFFT/psVectorFFT should perform a forward FFT.
-    PS_FFT_REVERSE = 1,                 ///< psImageFFT/psVectorFFT should perform a reverse FFT.
+    /// psImageFFT/psVectorFFT should perform a forward FFT.
+    PS_FFT_FORWARD = 0,
+
+    ///< psImageFFT/psVectorFFT should perform a reverse FFT.
+    PS_FFT_REVERSE = 1,
+
+    /** 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.
+    */
     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);
+psImage* psImageFFT(psImage* out, const psImage* in, int flags);
+psImage* psImageReal(psImage *out, const psImage* in);
+psImage* psImageImaginary(psImage *out, const psImage* in);
+psImage* psImageComplex(psImage* out, psImage *real, const psImage *imag);
+psImage* psImageConjugate(psImage *out, const psImage *in);
+psImage* psImagePowerSpectrum(psImage* out, const psImage* in);
+
+psVector* psVectorFFT(psVector* out, const psVector* in, int flags);
+psVector* psVectorReal(psVector* out, const psVector* in);
+psVector* psVectorImaginary(psVector* out, const psVector* in);
+psVector* psVectorComplex(psVector* out, psVector* real, const psVector* imag);
+psVector* psVectorConjugate(psVector* out, const psVector* in);
+psVector* psVectorPowerSpectrum(psVector* out, const psVector* in);