Index: /trunk/psLib/src/dataManip/psFFT.c =================================================================== --- /trunk/psLib/src/dataManip/psFFT.c (revision 1624) +++ /trunk/psLib/src/dataManip/psFFT.c (revision 1625) @@ -1,10 +1,10 @@ -/** @file psFFT.c +/** @file psVectorFFT.c * - * @brief Contains FFT transforms functions + * @brief Contains FFT transform related functions for psVector * * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-25 20:51:57 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-25 21:10:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -26,394 +26,4 @@ static bool p_fftwWisdomImported = false; - -psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction) -{ - unsigned int numCols; - unsigned int numRows; - 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 (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; - } - - 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); - - /* 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) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - - 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) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - - 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, const psImage* real, const psImage* imag) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - - 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) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - - 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) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - int numElementsSquared; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numCols = in->numCols; - numRows = in->numRows; - numElementsSquared = numCols * numCols * numRows * numRows; - - /* if not a complex number, this is not implemented */ - if (!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) Index: /trunk/psLib/src/dataManip/psFFT.h =================================================================== --- /trunk/psLib/src/dataManip/psFFT.h (revision 1624) +++ /trunk/psLib/src/dataManip/psFFT.h (revision 1625) @@ -1,6 +1,5 @@ - -/** @file psFFT.h +/** @file psVectorFFT.h * - * @brief Contains FFT transforms functions + * @brief Contains FFT transform related functions for psVector * * @ingroup Transform @@ -8,14 +7,13 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-10 18:54:38 $ + * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-25 21:10:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ -#ifndef PS_FFT_H -#define PS_FFT_H +#ifndef PS_VECTOR_FFT_H +#define PS_VECTOR_FFT_H -#include "psImage.h" #include "psVector.h" @@ -32,64 +30,4 @@ } psFftDirection; -/** Forward and reverse FFT calculations. - * - * This takes as input the image of interest (in) and the direction - * (direction), which is specified by an enumerated type psFftDirection. - * The input image may be of type psF32 or psC32, the result is always - * psC32. If the input vector is psF32, the direction must be forward. - * - * @return psImage* the FFT transformation result - */ -psImage* psImageFFT( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in, ///< the psImage to apply transform to - psFftDirection direction ///< the direction of the transform -); - -/** extract the real portion of a complex image - * - * @return psImage* real portion of the input image. - */ -psImage* psImageReal( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to extract real portion from -); - -/** extract the imaginary portion of a complex image - * - * @return psImage* imaginary portion of the input image. - */ -psImage* psImageImaginary( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to extract imaginary portion from -); - -/** creates a complex image from separate real and imaginary plane images - * - * @return psImage* resulting complex image - */ -psImage* psImageComplex( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* real, ///< the real plane image - const psImage* imag ///< the imaginary plane image -); - -/** computes the complex conjugate of an image - * - * @return psImage* the complex conjugate of the 'in' image - */ -psImage* psImageConjugate( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to compute conjugate of -); - -/** computes the power spectrum of an image - * - * @return psImage* the power spectrum of the 'in' image - */ -psImage* psImagePowerSpectrum( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to power spectrum of -); /** Forward and reverse FFT calculations. Index: /trunk/psLib/src/dataManip/psVectorFFT.c =================================================================== --- /trunk/psLib/src/dataManip/psVectorFFT.c (revision 1624) +++ /trunk/psLib/src/dataManip/psVectorFFT.c (revision 1625) @@ -1,10 +1,10 @@ -/** @file psFFT.c +/** @file psVectorFFT.c * - * @brief Contains FFT transforms functions + * @brief Contains FFT transform related functions for psVector * * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-25 20:51:57 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-25 21:10:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -26,394 +26,4 @@ static bool p_fftwWisdomImported = false; - -psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction) -{ - unsigned int numCols; - unsigned int numRows; - 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 (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; - } - - 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); - - /* 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) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - - 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) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - - 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, const psImage* real, const psImage* imag) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - - 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) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - - 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) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - int numElementsSquared; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numCols = in->numCols; - numRows = in->numRows; - numElementsSquared = numCols * numCols * numRows * numRows; - - /* if not a complex number, this is not implemented */ - if (!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) Index: /trunk/psLib/src/dataManip/psVectorFFT.h =================================================================== --- /trunk/psLib/src/dataManip/psVectorFFT.h (revision 1624) +++ /trunk/psLib/src/dataManip/psVectorFFT.h (revision 1625) @@ -1,6 +1,5 @@ - -/** @file psFFT.h +/** @file psVectorFFT.h * - * @brief Contains FFT transforms functions + * @brief Contains FFT transform related functions for psVector * * @ingroup Transform @@ -8,14 +7,13 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-10 18:54:38 $ + * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-25 21:10:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ -#ifndef PS_FFT_H -#define PS_FFT_H +#ifndef PS_VECTOR_FFT_H +#define PS_VECTOR_FFT_H -#include "psImage.h" #include "psVector.h" @@ -32,64 +30,4 @@ } psFftDirection; -/** Forward and reverse FFT calculations. - * - * This takes as input the image of interest (in) and the direction - * (direction), which is specified by an enumerated type psFftDirection. - * The input image may be of type psF32 or psC32, the result is always - * psC32. If the input vector is psF32, the direction must be forward. - * - * @return psImage* the FFT transformation result - */ -psImage* psImageFFT( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in, ///< the psImage to apply transform to - psFftDirection direction ///< the direction of the transform -); - -/** extract the real portion of a complex image - * - * @return psImage* real portion of the input image. - */ -psImage* psImageReal( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to extract real portion from -); - -/** extract the imaginary portion of a complex image - * - * @return psImage* imaginary portion of the input image. - */ -psImage* psImageImaginary( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to extract imaginary portion from -); - -/** creates a complex image from separate real and imaginary plane images - * - * @return psImage* resulting complex image - */ -psImage* psImageComplex( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* real, ///< the real plane image - const psImage* imag ///< the imaginary plane image -); - -/** computes the complex conjugate of an image - * - * @return psImage* the complex conjugate of the 'in' image - */ -psImage* psImageConjugate( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to compute conjugate of -); - -/** computes the power spectrum of an image - * - * @return psImage* the power spectrum of the 'in' image - */ -psImage* psImagePowerSpectrum( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to power spectrum of -); /** Forward and reverse FFT calculations. Index: /trunk/psLib/src/fft/psImageFFT.c =================================================================== --- /trunk/psLib/src/fft/psImageFFT.c (revision 1625) +++ /trunk/psLib/src/fft/psImageFFT.c (revision 1625) @@ -0,0 +1,414 @@ +/** @file psImageFFT.c + * + * @brief Contains FFT transform related functions for psImage. + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-25 21:10:09 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ +#include +#include +#include +#include +#include + +#include "psImageFFT.h" +#include "psError.h" +#include "psMemory.h" +#include "psLogMsg.h" +#include "psImageExtraction.h" + +#define P_FFTW_PLAN_RIGOR FFTW_ESTIMATE + +static bool p_fftwWisdomImported = false; + +psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction) +{ + unsigned int numCols; + unsigned int numRows; + 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 (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; + } + + 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); + + /* 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) +{ + psElemType type; + unsigned int numCols; + unsigned int numRows; + + 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) +{ + psElemType type; + unsigned int numCols; + unsigned int numRows; + + 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, const psImage* real, const psImage* imag) +{ + psElemType type; + unsigned int numCols; + unsigned int numRows; + + 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) +{ + psElemType type; + unsigned int numCols; + unsigned int numRows; + + 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) +{ + psElemType type; + unsigned int numCols; + unsigned int numRows; + int numElementsSquared; + + if (in == NULL) { + psFree(out); + return NULL; + } + + type = in->type.type; + numCols = in->numCols; + numRows = in->numRows; + numElementsSquared = numCols * numCols * numRows * numRows; + + /* if not a complex number, this is not implemented */ + if (!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; + +} Index: /trunk/psLib/src/fft/psImageFFT.h =================================================================== --- /trunk/psLib/src/fft/psImageFFT.h (revision 1625) +++ /trunk/psLib/src/fft/psImageFFT.h (revision 1625) @@ -0,0 +1,87 @@ +/** @file psImageFFT.h + * + * @brief Contains FFT transform related functions for psImage + * + * @ingroup Transform + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-25 21:10:09 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PS_IMAGE_FFT_H +#define PS_IMAGE_FFT_H + +#include "psImage.h" +#include "psVectorFFT.h" // for psFftDirection + +/// @addtogroup Transform +/// @{ + +/** Forward and reverse FFT calculations. + * + * This takes as input the image of interest (in) and the direction + * (direction), which is specified by an enumerated type psFftDirection. + * The input image may be of type psF32 or psC32, the result is always + * psC32. If the input vector is psF32, the direction must be forward. + * + * @return psImage* the FFT transformation result + */ +psImage* psImageFFT( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in, ///< the psImage to apply transform to + psFftDirection direction ///< the direction of the transform +); + +/** extract the real portion of a complex image + * + * @return psImage* real portion of the input image. + */ +psImage* psImageReal( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to extract real portion from +); + +/** extract the imaginary portion of a complex image + * + * @return psImage* imaginary portion of the input image. + */ +psImage* psImageImaginary( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to extract imaginary portion from +); + +/** creates a complex image from separate real and imaginary plane images + * + * @return psImage* resulting complex image + */ +psImage* psImageComplex( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* real, ///< the real plane image + const psImage* imag ///< the imaginary plane image +); + +/** computes the complex conjugate of an image + * + * @return psImage* the complex conjugate of the 'in' image + */ +psImage* psImageConjugate( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to compute conjugate of +); + +/** computes the power spectrum of an image + * + * @return psImage* the power spectrum of the 'in' image + */ +psImage* psImagePowerSpectrum( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to power spectrum of +); + +/// @} + +#endif Index: /trunk/psLib/src/fft/psVectorFFT.c =================================================================== --- /trunk/psLib/src/fft/psVectorFFT.c (revision 1624) +++ /trunk/psLib/src/fft/psVectorFFT.c (revision 1625) @@ -1,10 +1,10 @@ -/** @file psFFT.c +/** @file psVectorFFT.c * - * @brief Contains FFT transforms functions + * @brief Contains FFT transform related functions for psVector * * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.23 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-25 20:51:57 $ + * @version $Revision: 1.24 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-25 21:10:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -26,394 +26,4 @@ static bool p_fftwWisdomImported = false; - -psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction) -{ - unsigned int numCols; - unsigned int numRows; - 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 (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; - } - - 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); - - /* 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) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - - 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) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - - 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, const psImage* real, const psImage* imag) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - - 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) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - - 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) -{ - psElemType type; - unsigned int numCols; - unsigned int numRows; - int numElementsSquared; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numCols = in->numCols; - numRows = in->numRows; - numElementsSquared = numCols * numCols * numRows * numRows; - - /* if not a complex number, this is not implemented */ - if (!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) Index: /trunk/psLib/src/fft/psVectorFFT.h =================================================================== --- /trunk/psLib/src/fft/psVectorFFT.h (revision 1624) +++ /trunk/psLib/src/fft/psVectorFFT.h (revision 1625) @@ -1,6 +1,5 @@ - -/** @file psFFT.h +/** @file psVectorFFT.h * - * @brief Contains FFT transforms functions + * @brief Contains FFT transform related functions for psVector * * @ingroup Transform @@ -8,14 +7,13 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-10 18:54:38 $ + * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-25 21:10:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ -#ifndef PS_FFT_H -#define PS_FFT_H +#ifndef PS_VECTOR_FFT_H +#define PS_VECTOR_FFT_H -#include "psImage.h" #include "psVector.h" @@ -32,64 +30,4 @@ } psFftDirection; -/** Forward and reverse FFT calculations. - * - * This takes as input the image of interest (in) and the direction - * (direction), which is specified by an enumerated type psFftDirection. - * The input image may be of type psF32 or psC32, the result is always - * psC32. If the input vector is psF32, the direction must be forward. - * - * @return psImage* the FFT transformation result - */ -psImage* psImageFFT( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in, ///< the psImage to apply transform to - psFftDirection direction ///< the direction of the transform -); - -/** extract the real portion of a complex image - * - * @return psImage* real portion of the input image. - */ -psImage* psImageReal( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to extract real portion from -); - -/** extract the imaginary portion of a complex image - * - * @return psImage* imaginary portion of the input image. - */ -psImage* psImageImaginary( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to extract imaginary portion from -); - -/** creates a complex image from separate real and imaginary plane images - * - * @return psImage* resulting complex image - */ -psImage* psImageComplex( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* real, ///< the real plane image - const psImage* imag ///< the imaginary plane image -); - -/** computes the complex conjugate of an image - * - * @return psImage* the complex conjugate of the 'in' image - */ -psImage* psImageConjugate( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to compute conjugate of -); - -/** computes the power spectrum of an image - * - * @return psImage* the power spectrum of the 'in' image - */ -psImage* psImagePowerSpectrum( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in ///< the psImage to power spectrum of -); /** Forward and reverse FFT calculations. Index: /trunk/psLib/src/image/Makefile =================================================================== --- /trunk/psLib/src/image/Makefile (revision 1624) +++ /trunk/psLib/src/image/Makefile (revision 1625) @@ -3,6 +3,6 @@ ## Makefile: collections ## -## $Revision: 1.2 $ $Name: not supported by cvs2svn $ -## $Date: 2004-07-29 21:43:11 $ +## $Revision: 1.3 $ $Name: not supported by cvs2svn $ +## $Date: 2004-08-25 21:10:08 $ ## ## Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -34,5 +34,6 @@ psImageIO.o \ psImageManip.o \ - psImageStats.o + psImageStats.o \ + psImageFFT.o OBJS = $(addprefix makedir/,$(SRC_OBJS)) Index: /trunk/psLib/src/image/psImageFFT.c =================================================================== --- /trunk/psLib/src/image/psImageFFT.c (revision 1625) +++ /trunk/psLib/src/image/psImageFFT.c (revision 1625) @@ -0,0 +1,414 @@ +/** @file psImageFFT.c + * + * @brief Contains FFT transform related functions for psImage. + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-25 21:10:09 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ +#include +#include +#include +#include +#include + +#include "psImageFFT.h" +#include "psError.h" +#include "psMemory.h" +#include "psLogMsg.h" +#include "psImageExtraction.h" + +#define P_FFTW_PLAN_RIGOR FFTW_ESTIMATE + +static bool p_fftwWisdomImported = false; + +psImage* psImageFFT(psImage* out, const psImage* in, psFftDirection direction) +{ + unsigned int numCols; + unsigned int numRows; + 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 (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; + } + + 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); + + /* 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) +{ + psElemType type; + unsigned int numCols; + unsigned int numRows; + + 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) +{ + psElemType type; + unsigned int numCols; + unsigned int numRows; + + 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, const psImage* real, const psImage* imag) +{ + psElemType type; + unsigned int numCols; + unsigned int numRows; + + 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) +{ + psElemType type; + unsigned int numCols; + unsigned int numRows; + + 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) +{ + psElemType type; + unsigned int numCols; + unsigned int numRows; + int numElementsSquared; + + if (in == NULL) { + psFree(out); + return NULL; + } + + type = in->type.type; + numCols = in->numCols; + numRows = in->numRows; + numElementsSquared = numCols * numCols * numRows * numRows; + + /* if not a complex number, this is not implemented */ + if (!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; + +} Index: /trunk/psLib/src/image/psImageFFT.h =================================================================== --- /trunk/psLib/src/image/psImageFFT.h (revision 1625) +++ /trunk/psLib/src/image/psImageFFT.h (revision 1625) @@ -0,0 +1,87 @@ +/** @file psImageFFT.h + * + * @brief Contains FFT transform related functions for psImage + * + * @ingroup Transform + * + * @author Robert DeSonia, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-25 21:10:09 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PS_IMAGE_FFT_H +#define PS_IMAGE_FFT_H + +#include "psImage.h" +#include "psVectorFFT.h" // for psFftDirection + +/// @addtogroup Transform +/// @{ + +/** Forward and reverse FFT calculations. + * + * This takes as input the image of interest (in) and the direction + * (direction), which is specified by an enumerated type psFftDirection. + * The input image may be of type psF32 or psC32, the result is always + * psC32. If the input vector is psF32, the direction must be forward. + * + * @return psImage* the FFT transformation result + */ +psImage* psImageFFT( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in, ///< the psImage to apply transform to + psFftDirection direction ///< the direction of the transform +); + +/** extract the real portion of a complex image + * + * @return psImage* real portion of the input image. + */ +psImage* psImageReal( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to extract real portion from +); + +/** extract the imaginary portion of a complex image + * + * @return psImage* imaginary portion of the input image. + */ +psImage* psImageImaginary( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to extract imaginary portion from +); + +/** creates a complex image from separate real and imaginary plane images + * + * @return psImage* resulting complex image + */ +psImage* psImageComplex( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* real, ///< the real plane image + const psImage* imag ///< the imaginary plane image +); + +/** computes the complex conjugate of an image + * + * @return psImage* the complex conjugate of the 'in' image + */ +psImage* psImageConjugate( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to compute conjugate of +); + +/** computes the power spectrum of an image + * + * @return psImage* the power spectrum of the 'in' image + */ +psImage* psImagePowerSpectrum( + psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage* in ///< the psImage to power spectrum of +); + +/// @} + +#endif Index: /trunk/psLib/src/pslib.h =================================================================== --- /trunk/psLib/src/pslib.h (revision 1624) +++ /trunk/psLib/src/pslib.h (revision 1625) @@ -9,6 +9,6 @@ * @author Eric Van Alst, MHPCC * -* @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-08-25 20:51:57 $ +* @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-08-25 21:10:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -128,4 +128,5 @@ /// @ingroup DataManip #include "psVectorFFT.h" +#include "psImageFFT.h" #include "psFunctions.h"