Index: trunk/psLib/src/dataManip/psFFT.c
===================================================================
--- trunk/psLib/src/dataManip/psFFT.c	(revision 1264)
+++ trunk/psLib/src/dataManip/psFFT.c	(revision 1349)
@@ -1,13 +1,13 @@
 /** @file  psFFT.c
- *
- *  @brief Contains FFT transforms functions
- *
- *  @author Robert DeSonia, MHPCC
- *
- *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-07-22 20:42:54 $
- *
- *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
- */
+*
+*  @brief Contains FFT transforms functions
+*
+*  @author Robert DeSonia, MHPCC
+*
+*  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2004-07-30 04:03:51 $
+*
+*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+*/
 
 #include <unistd.h>
@@ -27,5 +27,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;
@@ -33,327 +33,327 @@
     psElemType type;
     fftwf_plan plan;
-
+    
     /* 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 ( 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;
-
-    }
-
+            psError( __func__, "Input image must be complex image for reverse FFT (type=%d).",
+                     type );
+            psFree( out );
+            return NULL;
+            
+        }
+        
     /* make sure the system-level wisdom information is imported. */
-    if (! p_fftwWisdomImported) {
-        fftwf_import_system_wisdom();
-        p_fftwWisdomImported = true;
-    }
-
+    if ( ! p_fftwWisdomImported ) {
+            fftwf_import_system_wisdom();
+            p_fftwWisdomImported = true;
+        }
+        
     numRows = in->numRows;
     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);
-
+    
+    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;
-    }
-
+    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;
     unsigned int numCols;
     unsigned int numRows;
-
-
-    if (in == NULL) {
-        psFree(out);
-        return NULL;
-    }
-
+    
+    
+    if ( in == NULL ) {
+            psFree( out );
+            return NULL;
+        }
+        
     type = in->type.type;
     numCols = in->numCols;
     numRows = in->numRows;
-
+    
     /* 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 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)
+    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 )
 {
     psElemType type;
     unsigned int numCols;
     unsigned int numRows;
-
-
-    if (in == NULL) {
-        psFree(out);
-        return NULL;
-    }
-
+    
+    
+    if ( in == NULL ) {
+            psFree( out );
+            return NULL;
+        }
+        
     type = in->type.type;
     numCols = in->numCols;
     numRows = in->numRows;
-
+    
     /* if not a complex number, this is logically just zeroed image of same size */
-    if (! 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);
-        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 ( ! 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 );
+            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 )
 {
     psElemType type;
     unsigned int numCols;
     unsigned int numRows;
-
-
-    if (real == NULL || imag == NULL) {
-        psFree(out);
-        return NULL;
-    }
-
+    
+    
+    if ( real == NULL || imag == NULL ) {
+            psFree( out );
+            return NULL;
+        }
+        
     type = real->type.type;
     numCols = real->numCols;
     numRows = real->numRows;
-
-    if (imag->type.type != type) {
-        psError(__func__,"The inputs to psImageComplex must be the same type.");
-        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;
     unsigned int numCols;
     unsigned int numRows;
-
-
-    if (in == NULL) {
-        psFree(out);
-        return NULL;
-    }
-
+    
+    
+    if ( in == NULL ) {
+            psFree( out );
+            return NULL;
+        }
+        
     type = in->type.type;
     numCols = in->numCols;
     numRows = in->numRows;
-
+    
     /* if not a complex number, this is logically just a image copy */
-    if (! 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)
+    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 )
 {
     psElemType type;
@@ -361,339 +361,346 @@
     unsigned int numRows;
     int numElementsSquared;
-
-    if (in == NULL) {
-        psFree(out);
-        return NULL;
-    }
-
+    
+    if ( in == NULL ) {
+            psFree( out );
+            return NULL;
+        }
+        
     type = in->type.type;
     numCols = in->numCols;
     numRows = in->numRows;
-    numElementsSquared = numCols*numCols*numRows*numRows;
-
+    numElementsSquared = numCols * numCols * numRows * numRows;
+    
     /* 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;
-        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;
-    }
-
-    return out;
-
+    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;
+        }
+        
+    return out;
+    
 }
 
 /************************************** Vector Functions ***************************************/
 
-psVector* psVectorFFT(psVector* out, const psVector* in, psFftDirection direction)
+psVector* psVectorFFT( psVector* out, const psVector* in, psFftDirection direction )
 {
     unsigned int numElements;
     psElemType type;
     fftwf_plan plan;
-
+    
     /* 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 ( 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) {
-        fftwf_import_system_wisdom();
-        p_fftwWisdomImported = true;
-    }
-
+    if ( ! p_fftwWisdomImported ) {
+            fftwf_import_system_wisdom();
+            p_fftwWisdomImported = true;
+        }
+        
     numElements = in->n;
-
-    out = psVectorRecycle(out, PS_TYPE_C32, numElements);
+    
+    out = psVectorRecycle( out, PS_TYPE_C32, numElements );
     out->n = numElements;
-
-    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);
-
+    
+    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);
-        return NULL;
-    }
-
+    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);
-        return NULL;
-    }
-
+    
+    if ( in == NULL ) {
+            psFree( out );
+            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);
-        out->n = numElements;
-        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,PS_TYPE_F32,numElements);
-        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)
+    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 );
+            out->n = numElements;
+            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, PS_TYPE_F32, numElements );
+            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 )
 {
     psElemType type;
     unsigned int numElements;
-
-
-    if (in == NULL) {
-        psFree(out);
-        return NULL;
-    }
-
+    
+    
+    if ( in == NULL ) {
+            psFree( out );
+            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,type,numElements);
-        out->n = 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,PS_TYPE_F32,numElements);
-        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)
+    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, type, numElements );
+            out->n = 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, PS_TYPE_F32, numElements );
+            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 )
 {
     psElemType type;
     unsigned int numElements;
-
-
-    if (real == NULL || imag == NULL) {
-        psFree(out);
-        return NULL;
-    }
-
+    
+    
+    if ( real == NULL || imag == NULL ) {
+            psFree( out );
+            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.");
-        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,PS_TYPE_C32,numElements);
-        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)
+    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." );
+            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, PS_TYPE_C32, numElements );
+            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 )
 {
     psElemType type;
     unsigned int numElements;
-
-
-    if (in == NULL) {
-        psFree(out);
-        return NULL;
-    }
-
+    
+    
+    if ( in == NULL ) {
+            psFree( out );
+            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);
-        out->n = 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,PS_TYPE_C32, numElements);
-        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)
+    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 );
+            out->n = 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, PS_TYPE_C32, numElements );
+            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 )
 {
     psElemType type;
@@ -702,55 +709,55 @@
     unsigned int inHalfNumElements;
     unsigned int inNumElementsSquared;
-
-    if (in == NULL) {
-        psFree(out);
-        return NULL;
-    }
-
+    
+    if ( in == NULL ) {
+            psFree( out );
+            return NULL;
+        }
+        
     type = in->type.type;
     inNumElements = in->n;
-    inNumElementsSquared = inNumElements*inNumElements;
-    inHalfNumElements = inNumElements/2;
-    outNumElements = inHalfNumElements+1;
-
+    inNumElementsSquared = inNumElements * inNumElements;
+    inHalfNumElements = inNumElements / 2;
+    outNumElements = inHalfNumElements + 1;
+    
     /* 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;
-        psF32 inAbs1;
-        psF32 inAbs2;
-
-        out = psVectorRecycle(out,PS_TYPE_F32,outNumElements);
-        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;
-
-        // 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;
-        }
-
-        // 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;
-
-}
+    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, PS_TYPE_F32, outNumElements );
+            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;
+            
+            // 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;
+                }
+                
+            // 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;
+    
+}