Index: trunk/psLib/src/fft/psVectorFFT.c
===================================================================
--- trunk/psLib/src/fft/psVectorFFT.c	(revision 1406)
+++ trunk/psLib/src/fft/psVectorFFT.c	(revision 1407)
@@ -1,2 +1,3 @@
+
 /** @file  psFFT.c
 *
@@ -5,6 +6,6 @@
 *  @author Robert DeSonia, MHPCC
 *
-*  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2004-08-06 22:34:05 $
+*  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-08-07 00:06:06 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -27,5 +28,5 @@
 static bool p_fftwWisdomImported = false;
 
-psImage* psImageFFT( psImage* out, const psImage* in, psFftDirection direction )
+psImage *psImageFFT(psImage * out, const psImage * in, psFftDirection direction)
 {
     unsigned int numCols;
@@ -35,35 +36,32 @@
 
     /* got good image data? */
-    if ( in == NULL ) {
-        psFree( out );
-        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 );
-        psFree( out );
-        return NULL;
-    }
-
-    if ( type != PS_TYPE_C32 && direction == PS_FFT_REVERSE ) {
-        psError( __func__, "Input image must be complex image for reverse FFT (type=%d).",
-                 type );
-        psFree( out );
-        return NULL;
-
-    }
-
-    if ( type != PS_TYPE_F32 && direction == PS_FFT_FORWARD ) {
-        psError( __func__, "Input image must be real image for forward FFT (type=%d).",
-                 type );
-        psFree( out );
+    if (in == NULL) {
+        psFree(out);
+        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);
+        psFree(out);
+        return NULL;
+    }
+
+    if (type != PS_TYPE_C32 && direction == PS_FFT_REVERSE) {
+        psError(__func__, "Input image must be complex image for reverse FFT (type=%d).", type);
+        psFree(out);
+        return NULL;
+
+    }
+
+    if (type != PS_TYPE_F32 && direction == PS_FFT_FORWARD) {
+        psError(__func__, "Input image must be real image for forward FFT (type=%d).", type);
+        psFree(out);
         return NULL;
     }
 
     /* make sure the system-level wisdom information is imported. */
-    if ( ! p_fftwWisdomImported ) {
+    if (!p_fftwWisdomImported) {
         fftwf_import_system_wisdom();
         p_fftwWisdomImported = true;
@@ -73,30 +71,27 @@
     numCols = in->numCols;
 
-    out = psImageCopy( out, in, PS_TYPE_C32 );
-
-    plan = fftwf_plan_dft_2d( numCols, numRows,
-                              ( fftwf_complex* ) out->data.C32[ 0 ],
-                              ( fftwf_complex* ) out->data.C32[ 0 ],
-                              direction,
-                              P_FFTW_PLAN_RIGOR );
-
-    /* check if a plan exists now*/
-    if ( plan == NULL ) {
-        psError( __func__, "Failed to create FFTW plan." );
-        psFree( out );
+    out = psImageCopy(out, in, PS_TYPE_C32);
+
+    plan = fftwf_plan_dft_2d(numCols, numRows,
+                             (fftwf_complex *) out->data.C32[0],
+                             (fftwf_complex *) out->data.C32[0], direction, P_FFTW_PLAN_RIGOR);
+
+    /* check if a plan exists now */
+    if (plan == NULL) {
+        psError(__func__, "Failed to create FFTW plan.");
+        psFree(out);
         return NULL;
     }
 
     /* finally, call FFTW with the plan made above */
-    fftwf_execute( plan );
-
-    fftwf_destroy_plan( plan );
-
-    return out;
-
-}
-
-
-psImage *psImageReal( psImage *out, const psImage* in )
+    fftwf_execute(plan);
+
+    fftwf_destroy_plan(plan);
+
+    return out;
+
+}
+
+psImage *psImageReal(psImage * out, const psImage * in)
 {
     psElemType type;
@@ -104,7 +99,6 @@
     unsigned int numRows;
 
-
-    if ( in == NULL ) {
-        psFree( out );
+    if (in == NULL) {
+        psFree(out);
         return NULL;
     }
@@ -115,49 +109,47 @@
 
     /* if not a complex number, this is logically just a copy */
-    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
+    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 an image copy was performed." );
-        return psImageCopy( out, in, type );
-    }
-
-    if ( type == PS_TYPE_C32 ) {
-        psF32 * outRow;
-        psC32* inRow;
-
-        out = psImageRecycle( 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 = psImageRecycle( 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 );
-            psFree( out );
-            return NULL;
-        }
-
-    return out;
-}
-
-psImage *psImageImaginary( psImage *out, const psImage* in )
+        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 = psImageRecycle(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 = psImageRecycle(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);
+        psFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImageImaginary(psImage * out, const psImage * in)
 {
     psElemType type;
@@ -165,7 +157,6 @@
     unsigned int numRows;
 
-
-    if ( in == NULL ) {
-        psFree( out );
+    if (in == NULL) {
+        psFree(out);
         return NULL;
     }
@@ -176,51 +167,49 @@
 
     /* if not a complex number, this is logically just zeroed image of same size */
-    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
+    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 zero image was returned." );
-        out = psImageRecycle( out, numCols, numRows, type );
-        memset( out->data.V[ 0 ], 0, PSELEMTYPE_SIZEOF( type ) * numCols * numRows );
+        psLogMsg(__func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
+                 "A zero image was returned.");
+        out = psImageRecycle(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 = psImageRecycle( 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 = psImageRecycle( 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 );
-            psFree( out );
-            return NULL;
-        }
-
-    return out;
-}
-
-psImage *psImageComplex( psImage* out, psImage *real, const psImage *imag )
+    if (type == PS_TYPE_C32) {
+        psF32 *outRow;
+        psC32 *inRow;
+
+        out = psImageRecycle(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 = psImageRecycle(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);
+        psFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImageComplex(psImage * out, psImage * real, const psImage * imag)
 {
     psElemType type;
@@ -228,7 +217,6 @@
     unsigned int numRows;
 
-
-    if ( real == NULL || imag == NULL ) {
-        psFree( out );
+    if (real == NULL || imag == NULL) {
+        psFree(out);
         return NULL;
     }
@@ -238,72 +226,69 @@
     numRows = real->numRows;
 
-    if ( imag->type.type != type ) {
-        psError( __func__, "The inputs to psImageComplex must be the same type." );
-        psFree( out );
-        return NULL;
-    }
-
-    if ( imag->numCols != numCols ||
-            imag->numRows != numRows ) {
-        psError( __func__, "The inputs to psImageComplex must be the same dimensions." );
-        psFree( out );
-        return NULL;
-    }
-
-    if ( PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
-        psError( __func__, "The inputs to psImageComplex can not be complex." );
-        psFree( out );
-        return NULL;
-    }
-
-    if ( type != PS_TYPE_F32 && type != PS_TYPE_F64 ) {
-        psError( __func__, "The input type to psImageComplex must be a floating point." );
-        psFree( out );
-        return NULL;
-    }
-
-    if ( type == PS_TYPE_F32 ) {
-        psC32 * outRow;
-        psF32* realRow;
-        psF32* imagRow;
-
-        out = psImageRecycle( out, 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_F64 ) {
-            psC64 * outRow;
-            psF64* realRow;
-            psF64* imagRow;
-
-            out = psImageRecycle( out, 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 );
-            psFree( out );
-            return NULL;
-        }
-
-    return out;
-}
-
-psImage *psImageConjugate( psImage *out, const psImage *in )
+    if (imag->type.type != type) {
+        psError(__func__, "The inputs to psImageComplex must be the same type.");
+        psFree(out);
+        return NULL;
+    }
+
+    if (imag->numCols != numCols || imag->numRows != numRows) {
+        psError(__func__, "The inputs to psImageComplex must be the same dimensions.");
+        psFree(out);
+        return NULL;
+    }
+
+    if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__, "The inputs to psImageComplex can not be complex.");
+        psFree(out);
+        return NULL;
+    }
+
+    if (type != PS_TYPE_F32 && type != PS_TYPE_F64) {
+        psError(__func__, "The input type to psImageComplex must be a floating point.");
+        psFree(out);
+        return NULL;
+    }
+
+    if (type == PS_TYPE_F32) {
+        psC32 *outRow;
+        psF32 *realRow;
+        psF32 *imagRow;
+
+        out = psImageRecycle(out, 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_F64) {
+        psC64 *outRow;
+        psF64 *realRow;
+        psF64 *imagRow;
+
+        out = psImageRecycle(out, 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);
+        psFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImageConjugate(psImage * out, const psImage * in)
 {
     psElemType type;
@@ -311,7 +296,6 @@
     unsigned int numRows;
 
-
-    if ( in == NULL ) {
-        psFree( out );
+    if (in == NULL) {
+        psFree(out);
         return NULL;
     }
@@ -322,49 +306,47 @@
 
     /* if not a complex number, this is logically just a image copy */
-    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
+    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. "
-                  "Image copy was performed instead." );
-        return psImageCopy( out, in, type );
-    }
-
-    if ( type == PS_TYPE_C32 ) {
-        psC32 * outRow;
-        psC32* inRow;
-
-        out = psImageRecycle( 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 = psImageRecycle( 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 );
-            psFree( out );
-            return NULL;
-        }
-
-    return out;
-}
-
-psImage *psImagePowerSpectrum( psImage* out, const psImage* in )
+        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 = psImageRecycle(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 = psImageRecycle(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);
+        psFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psImage *psImagePowerSpectrum(psImage * out, const psImage * in)
 {
     psElemType type;
@@ -373,6 +355,6 @@
     int numElementsSquared;
 
-    if ( in == NULL ) {
-        psFree( out );
+    if (in == NULL) {
+        psFree(out);
         return NULL;
     }
@@ -384,53 +366,49 @@
 
     /* 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." );
-        psFree( out );
-        return NULL;
-    }
-
-    if ( type == PS_TYPE_C32 ) {
-        psF32 * outRow;
-        psC32* inRow;
+    if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__, "Power Spectrum for non-complex inputs is not implemented.");
+        psFree(out);
+        return NULL;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32 *outRow;
+        psC32 *inRow;
         psF32 real;
         psF32 imag;
 
-
-        out = psImageRecycle( 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 = psImageRecycle( 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 );
-            psFree( out );
-            return NULL;
-        }
+        out = psImageRecycle(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 = psImageRecycle(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);
+        psFree(out);
+        return NULL;
+    }
 
     return out;
@@ -440,5 +418,5 @@
 /************************************** Vector Functions ***************************************/
 
-psVector* psVectorFFT( psVector* out, const psVector* in, psFftDirection direction )
+psVector *psVectorFFT(psVector * out, const psVector * in, psFftDirection direction)
 {
     unsigned int numElements;
@@ -447,35 +425,32 @@
 
     /* got good image data? */
-    if ( in == NULL ) {
-        psFree( out );
-        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 );
-        psFree( out );
-        return NULL;
-    }
-
-    if ( ( type != PS_TYPE_C32 ) && ( direction == PS_FFT_REVERSE ) ) {
-        psError( __func__, "Input image must be complex image for reverse FFT (type=%d).",
-                 type );
-        psFree( out );
-        return NULL;
-
-    }
-
-    if ( ( type != PS_TYPE_F32 ) && ( direction == PS_FFT_FORWARD ) ) {
-        psError( __func__, "Input image must be real image for forward FFT (type=%d).",
-                 type );
-        psFree( out );
+    if (in == NULL) {
+        psFree(out);
+        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);
+        psFree(out);
+        return NULL;
+    }
+
+    if ((type != PS_TYPE_C32) && (direction == PS_FFT_REVERSE)) {
+        psError(__func__, "Input image must be complex image for reverse FFT (type=%d).", type);
+        psFree(out);
+        return NULL;
+
+    }
+
+    if ((type != PS_TYPE_F32) && (direction == PS_FFT_FORWARD)) {
+        psError(__func__, "Input image must be real image for forward FFT (type=%d).", type);
+        psFree(out);
         return NULL;
     }
 
     /* make sure the system-level wisdom information is imported. */
-    if ( ! p_fftwWisdomImported ) {
+    if (!p_fftwWisdomImported) {
         fftwf_import_system_wisdom();
         p_fftwWisdomImported = true;
@@ -484,51 +459,50 @@
     numElements = in->n;
 
-    out = psVectorRecycle( out, numElements, PS_TYPE_C32 );
+    out = psVectorRecycle(out, numElements, PS_TYPE_C32);
     out->n = numElements;
 
-    if ( type == PS_TYPE_F32 ) {
+    if (type == PS_TYPE_F32) {
         // need to convert to complex
-        psC32 * outVec = out->data.C32;
-        psF32* inVec = in->data.F32;
-        for ( unsigned int i = 0;i < numElements;i++ ) {
-            outVec[ i ] = inVec[ i ];
-        }
-    } else {
-        psC32* outVec = out->data.C32;
-        psC32* inVec = in->data.C32;
-        for ( unsigned int i = 0;i < numElements;i++ ) {
-            outVec[ i ] = inVec[ i ];
-        }
-    }
-
-    plan = fftwf_plan_dft_1d( numElements,
-                              ( fftwf_complex* ) out->data.C32,
-                              ( fftwf_complex* ) out->data.C32,
-                              direction,
-                              P_FFTW_PLAN_RIGOR );
-
-    /* check if a plan exists now*/
-    if ( plan == NULL ) {
-        psError( __func__, "Failed to create FFTW plan." );
-        psFree( out );
+        psC32 *outVec = out->data.C32;
+        psF32 *inVec = in->data.F32;
+
+        for (unsigned int i = 0; i < numElements; i++) {
+            outVec[i] = inVec[i];
+        }
+    } else {
+        psC32 *outVec = out->data.C32;
+        psC32 *inVec = in->data.C32;
+
+        for (unsigned int i = 0; i < numElements; i++) {
+            outVec[i] = inVec[i];
+        }
+    }
+
+    plan = fftwf_plan_dft_1d(numElements,
+                             (fftwf_complex *) out->data.C32,
+                             (fftwf_complex *) out->data.C32, direction, P_FFTW_PLAN_RIGOR);
+
+    /* check if a plan exists now */
+    if (plan == NULL) {
+        psError(__func__, "Failed to create FFTW plan.");
+        psFree(out);
         return NULL;
     }
 
     /* finally, call FFTW with the plan made above */
-    fftwf_execute( plan );
-
-    fftwf_destroy_plan( plan );
-
-    return out;
-}
-
-
-psVector *psVectorReal( psVector *out, const psVector* in )
+    fftwf_execute(plan);
+
+    fftwf_destroy_plan(plan);
+
+    return out;
+}
+
+psVector *psVectorReal(psVector * out, const psVector * in)
 {
     psElemType type;
     unsigned int numElements;
 
-    if ( in == NULL ) {
-        psFree( out );
+    if (in == NULL) {
+        psFree(out);
         return NULL;
     }
@@ -538,43 +512,41 @@
 
     /* if not a complex number, this is logically just a copy */
-    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
+    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, numElements, type );
+        psLogMsg(__func__, PS_LOG_WARN, "Real portion of a non-Complex type called called for. "
+                 "Just a vector copy was performed.");
+        out = psVectorRecycle(out, numElements, type);
         out->n = numElements;
-        memcpy( out->data.V, in->data.V, numElements * PSELEMTYPE_SIZEOF( type ) );
+        memcpy(out->data.V, in->data.V, numElements * PSELEMTYPE_SIZEOF(type));
         return out;
     }
 
-    if ( type == PS_TYPE_C32 ) {
-        psF32 * outVec;
-        psC32* inVec = in->data.C32;
-
-        out = psVectorRecycle( out, numElements, PS_TYPE_F32 );
+    if (type == PS_TYPE_C32) {
+        psF32 *outVec;
+        psC32 *inVec = in->data.C32;
+
+        out = psVectorRecycle(out, numElements, PS_TYPE_F32);
         out->n = numElements;
         outVec = out->data.F32;
 
-        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 );
-        psFree( out );
-        return NULL;
-    }
-
-    return out;
-}
-
-psVector *psVectorImaginary( psVector *out, const psVector* in )
+        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);
+        psFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorImaginary(psVector * out, const psVector * in)
 {
     psElemType type;
     unsigned int numElements;
 
-
-    if ( in == NULL ) {
-        psFree( out );
+    if (in == NULL) {
+        psFree(out);
         return NULL;
     }
@@ -584,48 +556,46 @@
 
     /* if not a complex number, this is logically just zeroed image of same size */
-    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
+    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 );
+        psLogMsg(__func__, PS_LOG_WARN, "Imaginary portion of a non-Complex type called for. "
+                 "A zeroed vector was returned.");
+        out = psVectorRecycle(out, numElements, type);
         out->n = numElements;
-        memset( out->data.V, 0, PSELEMTYPE_SIZEOF( type ) * numElements );
+        memset(out->data.V, 0, PSELEMTYPE_SIZEOF(type) * numElements);
         return out;
     }
 
-    if ( type == PS_TYPE_C32 ) {
-        psF32 * outVec;
-        psC32* inVec = in->data.C32;
-
-        out = psVectorRecycle( out, numElements, PS_TYPE_F32 );
+    if (type == PS_TYPE_C32) {
+        psF32 *outVec;
+        psC32 *inVec = in->data.C32;
+
+        out = psVectorRecycle(out, numElements, PS_TYPE_F32);
         out->n = numElements;
         outVec = out->data.F32;
 
-        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 );
-        psFree( out );
-        return NULL;
-    }
-
-    return out;
-}
-
-psVector *psVectorComplex( psVector* out, psVector *real, const psVector *imag )
+        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);
+        psFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorComplex(psVector * out, psVector * real, const psVector * imag)
 {
     psElemType type;
     unsigned int numElements;
 
-
-    if ( real == NULL || imag == NULL ) {
-        psFree( out );
+    if (real == NULL || imag == NULL) {
+        psFree(out);
         return NULL;
     }
 
     type = real->type.type;
-    if ( real->n < imag->n ) {
+    if (real->n < imag->n) {
         numElements = real->n;
     } else {
@@ -633,46 +603,44 @@
     }
 
-    if ( imag->type.type != type ) {
-        psError( __func__, "The inputs to psVectorComplex must be the same type." );
-        psFree( out );
-        return NULL;
-    }
-
-    if ( PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
-        psError( __func__, "The inputs to psVectorComplex can not be complex." );
-        psFree( out );
-        return NULL;
-    }
-
-    if ( type == PS_TYPE_F32 ) {
-        psC32 * outVec;
-        psF32* realVec = real->data.F32;
-        psF32* imagVec = imag->data.F32;
-
-        out = psVectorRecycle( out, numElements, PS_TYPE_C32 );
+    if (imag->type.type != type) {
+        psError(__func__, "The inputs to psVectorComplex must be the same type.");
+        psFree(out);
+        return NULL;
+    }
+
+    if (PS_IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__, "The inputs to psVectorComplex can not be complex.");
+        psFree(out);
+        return NULL;
+    }
+
+    if (type == PS_TYPE_F32) {
+        psC32 *outVec;
+        psF32 *realVec = real->data.F32;
+        psF32 *imagVec = imag->data.F32;
+
+        out = psVectorRecycle(out, numElements, PS_TYPE_C32);
         out->n = numElements;
         outVec = out->data.C32;
 
-        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 );
-        psFree( out );
-        return NULL;
-    }
-
-    return out;
-}
-
-psVector *psVectorConjugate( psVector *out, const psVector *in )
+        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);
+        psFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorConjugate(psVector * out, const psVector * in)
 {
     psElemType type;
     unsigned int numElements;
 
-
-    if ( in == NULL ) {
-        psFree( out );
+    if (in == NULL) {
+        psFree(out);
         return NULL;
     }
@@ -682,37 +650,36 @@
 
     /* if not a complex number, this is logically just a image copy */
-    if ( ! PS_IS_PSELEMTYPE_COMPLEX( type ) ) {
+    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, numElements, type );
+        psLogMsg(__func__, PS_LOG_WARN, "Complex Conjugate of a non-Complex type called for. "
+                 "Vector copy was performed instead.");
+
+        out = psVectorRecycle(out, numElements, type);
         out->n = numElements;
-        memcpy( out->data.V, in->data.V, PSELEMTYPE_SIZEOF( type ) * numElements );
+        memcpy(out->data.V, in->data.V, PSELEMTYPE_SIZEOF(type) * numElements);
         return out;
     }
 
-    if ( type == PS_TYPE_C32 ) {
-        psC32 * outVec;
-        psC32* inVec = in->data.C32;
-
-        out = psVectorRecycle( out, numElements, PS_TYPE_C32 );
+    if (type == PS_TYPE_C32) {
+        psC32 *outVec;
+        psC32 *inVec = in->data.C32;
+
+        out = psVectorRecycle(out, numElements, PS_TYPE_C32);
         out->n = numElements;
         outVec = out->data.C32;
 
-        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 );
-        psFree( out );
-        return NULL;
-    }
-
-    return out;
-}
-
-psVector *psVectorPowerSpectrum( psVector* out, const psVector* in )
+        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);
+        psFree(out);
+        return NULL;
+    }
+
+    return out;
+}
+
+psVector *psVectorPowerSpectrum(psVector * out, const psVector * in)
 {
     psElemType type;
@@ -722,6 +689,6 @@
     unsigned int inNumElementsSquared;
 
-    if ( in == NULL ) {
-        psFree( out );
+    if (in == NULL) {
+        psFree(out);
         return NULL;
     }
@@ -734,42 +701,41 @@
 
     /* 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." );
-        psFree( out );
-        return NULL;
-    }
-
-    if ( type == PS_TYPE_C32 ) {
-        psF32 * outVec;
-        psC32* inVec = in->data.C32;
+    if (!PS_IS_PSELEMTYPE_COMPLEX(type)) {
+        psError(__func__, "Power Spectrum for non-complex inputs is not implemented.");
+        psFree(out);
+        return NULL;
+    }
+
+    if (type == PS_TYPE_C32) {
+        psF32 *outVec;
+        psC32 *inVec = in->data.C32;
         psF32 inAbs1;
         psF32 inAbs2;
 
-        out = psVectorRecycle( out, outNumElements, PS_TYPE_F32 );
+        out = psVectorRecycle(out, outNumElements, PS_TYPE_F32);
         out->n = outNumElements;
         outVec = out->data.F32;
 
         // from ADD: P_0 = |C_0|^2/N^2
-        inAbs1 = cabsf( inVec[ 0 ] );
-        outVec[ 0 ] = inAbs1 * inAbs1 / inNumElementsSquared;
+        inAbs1 = cabsf(inVec[0]);
+        outVec[0] = inAbs1 * inAbs1 / inNumElementsSquared;
 
         // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1)
-        for ( unsigned int i = 1;i < inHalfNumElements;i++ ) {
-            inAbs1 = cabsf( inVec[ i ] );
-            inAbs2 = cabsf( inVec[ inNumElements - i ] );
-            outVec[ i ] = ( inAbs1 * inAbs1 + inAbs2 * inAbs2 ) / inNumElementsSquared;
+        for (unsigned int i = 1; i < inHalfNumElements; i++) {
+            inAbs1 = cabsf(inVec[i]);
+            inAbs2 = cabsf(inVec[inNumElements - i]);
+            outVec[i] = (inAbs1 * inAbs1 + inAbs2 * inAbs2) / inNumElementsSquared;
         }
 
         // from ADD: P_N/2 = |C_N/2|^2/N^2
-        inAbs1 = cabsf( inVec[ inHalfNumElements ] );
-        outVec[ inHalfNumElements ] = inAbs1 * inAbs1 / inNumElementsSquared;
-    } else {
-        psError( __func__, "Can not power spectrum for given vector type (%d).",
-                 type );
-        psFree( out );
-        return NULL;
-    }
-
-    return out;
-
-}
+        inAbs1 = cabsf(inVec[inHalfNumElements]);
+        outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared;
+    } else {
+        psError(__func__, "Can not power spectrum for given vector type (%d).", type);
+        psFree(out);
+        return NULL;
+    }
+
+    return out;
+
+}