Index: trunk/psLib/src/fft/psImageFFT.c =================================================================== --- trunk/psLib/src/fft/psImageFFT.c (revision 11668) +++ trunk/psLib/src/fft/psImageFFT.c (revision 11703) @@ -1,13 +1,14 @@ -/** @file psImageFFT.c - * - * @brief Contains FFT transform related functions for psImage. - * - * @author Robert DeSonia, MHPCC - * - * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-02-06 21:36:09 $ - * - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ +/// @file psImageFFT.c +/// +/// @brief Contains FFT transform related functions for psImage. +/// +/// @author Paul Price, IfA +/// @author Robert DeSonia, MHPCC +/// +/// @version $Revision: 1.21 $ $Name: not supported by cvs2svn $ +/// @date $Date: 2007-02-08 04:17:58 $ +/// +/// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +/// #ifdef HAVE_CONFIG_H @@ -20,442 +21,254 @@ #include -#include "psImageFFT.h" +#include "psAssert.h" #include "psError.h" #include "psMemory.h" #include "psLogMsg.h" +#include "psConstants.h" #include "psImageStructManip.h" - - - -#define PS_FFTW_PLAN_RIGOR FFTW_ESTIMATE - -static bool p_fftwWisdomImported = false; - -psImage* psImageFFT(psImage* out, const psImage* image, psFFTFlags direction) -{ - psU32 numCols; - psU32 numRows; - psElemType type; - fftwf_plan plan; - - /* got good image data? */ - if (image == NULL) { - psFree(out); +#include "psImageFFT.h" + + +#define FFTW_PLAN_RIGOR FFTW_ESTIMATE // How rigorous the FFTW planning is + +static psBool fftwWisdomImported = false; // Has the FFTW wisdom been imported yet? + +bool psImageForwardFFT(psImage **real, psImage **imag, const psImage *in) +{ + PS_ASSERT_IMAGE_NON_NULL(in, false); + PS_ASSERT_IMAGE_TYPE(in, PS_TYPE_F32, false); + PS_ASSERT_PTR_NON_NULL(real, false); + PS_ASSERT_PTR_NON_NULL(imag, false); + + // Make sure the system-level wisdom information is imported. + if (!fftwWisdomImported) { + fftwf_import_system_wisdom(); + fftwWisdomImported = true; + } + + int numCols = in->numCols; // Number of columns + int numRows = in->numRows; // Number of rows + + psF32 *input; // Input data, for FFTW + if (!in->parent) { + // No parent --- can just use the data buffer + input = psMemIncrRefCounter(in->p_rawDataBuffer); + } else { + // Need to copy the data + input = psAlloc(numCols * numRows); + for (int y = 0; y < numRows; y++) { + memcpy(&input[y * numRows], in->data.F32[y], numCols); + } + } + + // Do the FFT + fftwf_complex *out = fftwf_malloc((numCols/2 + 1) * numRows * sizeof(fftwf_complex)); // Output data + fftwf_plan plan = fftwf_plan_dft_r2c_2d(numRows, numCols, input, out, FFTW_PLAN_RIGOR); + fftwf_execute(plan); + fftwf_destroy_plan(plan); + psFree(input); + + // Pull the real and imaginary parts out + numCols = numCols/2 + 1; // x dimension is halved by FFTW + *real = psImageRecycle(*real, numCols, numRows, PS_TYPE_F32); + *imag = psImageRecycle(*imag, numCols, numRows, PS_TYPE_F32); + for (int y = 0, index = 0; y < numRows; y++) { + for (int x = 0; x < numCols; index++, x++) { +#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) + // C99 complex support + (*real)->data.F32[y][x] = creal(out[index]); + (*imag)->data.F32[y][x] = cimag(out[index]); +#else + // FFTW's backup complex support + (*real)->data.F32[y][x] = out[index][0]; + (*imag)->data.F32[y][x] = out[index][1]; +#endif + } + } + fftwf_free(out); + + return true; +} + +bool psImageBackwardFFT(psImage **out, const psImage *real, const psImage *imag, int origCols) +{ + PS_ASSERT_IMAGE_NON_NULL(real, false); + PS_ASSERT_IMAGE_TYPE(real, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(imag, false); + PS_ASSERT_IMAGE_TYPE(imag, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(real, imag, false); + PS_ASSERT_PTR_NON_NULL(out, false); + + int numCols = real->numCols; // Number of columns + int numRows = real->numRows; // Number of rows + + // Because of the way FFT r2c and c2r work, need the number of columns in the target to be: + // 2 * numCols = 2*(numCols/2 + 1) = origCols % 2 ? origCols + 1 : origCols + 2 + if (2 * numCols != ((origCols % 2) ? origCols + 1 : origCols + 2)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "Number of columns in FFT-ed images (%d) should be origCols (%d) / 2 + 1 = %d", + numCols, origCols, origCols/2 + 1); + return false; + } + + if (*out && (*out)->parent) { + // It has a parent, so we can't write directly into the buffer. + // It had better be the correct size and type, because we don't want to resize a child image + PS_ASSERT_IMAGE_SIZE(*out, origCols, numRows, false); + PS_ASSERT_IMAGE_TYPE(*out, PS_TYPE_F32, false); + } + + + // Make sure the system-level wisdom information is imported. + if (!fftwWisdomImported) { + fftwf_import_system_wisdom(); + fftwWisdomImported = true; + } + + // Stuff the real and imaginary parts in + psF32 *target = psAlloc(2 * numCols * numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Target for FFTW + fftwf_complex *in = fftwf_malloc(numCols * numRows * sizeof(fftwf_complex)); // Input data + for (int y = 0, index = 0; y < numRows; y++) { + for (int x = 0; x < numCols; index++, x++) { +#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) + // C99 complex support + in[index] = real->data.F32[y][x] + imag->data.F32[y][x] * I; +#else + // FFTW's backup complex support + in[index][0] = real->data.F32[y][x]; + in[index][1] = imag->data.F32[y][x]; +#endif + } + } + + // Do the FFT + fftwf_plan plan = fftwf_plan_dft_c2r_2d(numRows, origCols, in, target, FFTW_PLAN_RIGOR); + fftwf_execute(plan); + fftwf_destroy_plan(plan); + fftwf_free(in); + + if (!(*out) || !(*out)->parent) { + *out = psImageRecycle(*out, origCols, numRows, PS_TYPE_F32); + } + + // Copy the target pixels into the output + // There's a slight offset --- target has the additional padding required by FFTW, + // while the output doesn't. + for (int y = 0, index = 0; y < numRows; y++, index += origCols) { + memcpy(&(*out)->data.F32[y][0], &target[index], numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + } + psFree(target); + + return true; +} + +psImage* psImagePowerSpectrum(psImage *out, const psImage *in) +{ + PS_ASSERT_IMAGE_NON_NULL(in, NULL); + PS_ASSERT_IMAGE_TYPE(in, PS_TYPE_F32, NULL); + + psImage *real = NULL; // Real component of FFT + psImage *imag = NULL; // Imaginary component of FFT + + if (!psImageForwardFFT(&real, &imag, in)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform forward FFT."); return NULL; } - if ( ((direction & PS_FFT_FORWARD) != 0) ) { - if ((direction & PS_FFT_REVERSE) != 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Can not specify both PS_FFT_FORWARD and PS_FFT_REVERSE options.")); - psFree(out); - return NULL; - } - if ((direction & PS_FFT_REAL_RESULT) != 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("The PS_FFT_FORWARD and PS_FFT_REAL_RESULT combinition is not supported.")); - psFree(out); - return NULL; - } - } else if ((direction & PS_FFT_REVERSE) == 0) { + int numCols = real->numCols; // Number of columns + int numRows = real->numRows; // Number of rows + + float norm = 1.0 / PS_SQR(in->numCols) / PS_SQR(in->numRows); + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + // Power spectrum is the square of the complex modulus + real->data.F32[y][x] = norm * (PS_SQR(real->data.F32[y][x]) + PS_SQR(imag->data.F32[y][x])); + } + } + psFree(imag); + + return real; +} + + + +bool psImageComplexMultiply(psImage **outReal, psImage **outImag, + const psImage *in1Real, const psImage *in1Imag, + const psImage *in2Real, const psImage *in2Imag) +{ + PS_ASSERT_IMAGE_NON_NULL(in1Real, false); + PS_ASSERT_IMAGE_NON_NULL(in1Imag, false); + PS_ASSERT_IMAGE_NON_NULL(in2Real, false); + PS_ASSERT_IMAGE_NON_NULL(in2Imag, false); + PS_ASSERT_IMAGE_TYPE(in1Real, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_TYPE(in1Imag, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_TYPE(in2Real, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_TYPE(in2Imag, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(in1Imag, in1Real, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(in2Real, in1Real, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(in2Imag, in1Real, false); + PS_ASSERT_PTR_NON_NULL(outReal, false); + PS_ASSERT_PTR_NON_NULL(outImag, false); + if (*outReal) { + PS_ASSERT_IMAGE_NON_NULL(*outReal, false); + PS_ASSERT_IMAGE_NON_NULL(*outImag, false); + } else if (*outImag) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Must specify either PS_FFT_FORWARD or PS_FFT_REVERSE option.")); - psFree(out); - return NULL; - } - - type = image->type.type; - - /* make sure the system-level wisdom information is imported. */ - if (!p_fftwWisdomImported) { - fftwf_import_system_wisdom(); - p_fftwWisdomImported = true; - } - - numRows = image->numRows; - numCols = image->numCols; - - // n.b. FFTW can perform a in-place transform at the same rate or faster than out-of-place. - psS32 sign = ((direction & PS_FFT_FORWARD) != 0) ? FFTW_FORWARD : FFTW_BACKWARD; - - fftwf_complex* outBuffer = fftwf_malloc(numRows*numCols*sizeof(fftwf_complex)); - p_psImageCopyToRawBuffer(outBuffer, image, PS_TYPE_C32); - - plan = fftwf_plan_dft_2d(numRows, numCols, - outBuffer, - outBuffer, - sign, - PS_FFTW_PLAN_RIGOR); - - /* check if a plan exists now -- if not, it is a real problem at this point */ - if (plan == NULL) { - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Could not create a valid FFT plan to perform the transform.")); - psFree(out); - return NULL; - } - - /* finally, call FFTW with the plan made above */ - fftwf_execute(plan); - - fftwf_destroy_plan(plan); - - if (direction & PS_FFT_REAL_RESULT) { - // n.b., we do this instead of using fftwf_plan_dft_c2r because that - // plan requires a half-image, which would require a image reordering - // that is not as simple as performing a normal complex transform and - // then taking the real part of the result. If performance here - // becomes an issue, the use of fftwf_plan_dft_r2c should be considered - // as well as fftwf_plan_dft_c2r. - out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); - int index = 0; - for (int row=0; row < numRows; row++) { - psF32* outRow = out->data.F32[row]; - for (int col=0; col < numCols; col++) { - outRow[col] = crealf(outBuffer[index++]); // take just the real part - } - } + "If the output real part is provided, the output imaginary part must also be provided."); + return false; + } + + int numRows = in1Real->numRows; // Number of rows + int numCols = in1Real->numCols; // Number of columns + + // Need to worry if the outputs are children + psImage *targetReal, *targetImag; // Target real and imaginary parts + + if ((*outReal)->parent) { + // It had better be the correct size and type, because we don't want to resize a child image + PS_ASSERT_IMAGE_TYPE(*outReal, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_SIZE(*outReal, numCols, numRows, false); + targetReal = psImageAlloc(numCols, numRows, PS_TYPE_F32); } else { - out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); - int index = 0; - for (int row=0; row < numRows; row++) { - psC32* outRow = out->data.C32[row]; - for (int col=0; col < numCols; col++) { - outRow[col] = outBuffer[index++]; // take just the real part - } - } - // memcpy(out->rawDataBuffer, outBuffer, numRows*numCols*sizeof(fftwf_complex)); - } - - fftwf_free(outBuffer); - - return out; - -} - -psImage* psImageReal(psImage* out, const psImage* in) -{ - psElemType type; - psU32 numCols; - psU32 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 then */ - if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { - return psImageCopy(out, in, type); - } - - if (type == PS_TYPE_C32) { - psF32* outRow; - psC32* inRow; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.F32[row]; - inRow = in->data.C32[row]; - - for (psU32 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 (psU32 row = 0; row < numRows; row++) { - outRow = out->data.F64[row]; - inRow = in->data.C64[row]; - - for (psU32 col = 0; col < numCols; col++) { - outRow[col] = creal(inRow[col]); - } - } + *outReal = psImageRecycle(*outReal, numCols, numRows, PS_TYPE_F32); + targetReal = psMemIncrRefCounter(*outReal); + } + if ((*outImag)->parent) { + // It had better be the correct size and type, because we don't want to resize a child image + if ((*outReal)->numCols != numCols || (*outReal)->numRows != numRows || + (*outImag)->type.type != PS_TYPE_F32) { + // Plug potential memory leak --- need to free targetReal if there's a problem + psFree(targetReal); + } + PS_ASSERT_IMAGE_TYPE(*outReal, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_SIZE(*outReal, numCols, numRows, false); + targetImag = psImageAlloc(numCols, numRows, PS_TYPE_F32); } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psImage type, %s, is not supported."), - typeStr); - - psFree(out); - return NULL; - } - - return out; -} - -psImage* psImageImaginary(psImage* out, const psImage* in) -{ - psElemType type; - psU32 numCols; - psU32 numRows; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numCols = in->numCols; - numRows = in->numRows; - - /* if not a complex image type, this is logically just zeroed image of same size */ - if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { - 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 (psU32 row = 0; row < numRows; row++) { - outRow = out->data.F32[row]; - inRow = in->data.C32[row]; - - for (psU32 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 (psU32 row = 0; row < numRows; row++) { - outRow = out->data.F64[row]; - inRow = in->data.C64[row]; - - for (psU32 col = 0; col < numCols; col++) { - outRow[col] = cimag(inRow[col]); - } - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psImage type, %s, is not supported."), - typeStr); - psFree(out); - return NULL; - } - - return out; -} - -psImage* psImageComplex(psImage* out, const psImage* real, const psImage* imag) -{ - psElemType type; - psU32 numCols; - psU32 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) { - char* typeStrReal; - char* typeStrImag; - PS_TYPE_NAME(typeStrReal,type); - PS_TYPE_NAME(typeStrImag,imag->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Real psImage type (%s) and imaginary psImage type (%s) must be the same."), - typeStrReal,typeStrImag); - psFree(out); - return NULL; - } - - if (imag->numCols != numCols || imag->numRows != numRows) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Real psImage size (%dx%d) and imaginary psImage size (%dx%d) must be the same."), - numCols, numRows, imag->numCols, imag->numRows); - psFree(out); - return NULL; - } - - if (type == PS_TYPE_F32) { - psC32* outRow; - psF32* realRow; - psF32* imagRow; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32); - - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.C32[row]; - realRow = real->data.F32[row]; - imagRow = imag->data.F32[row]; - - for (psU32 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 (psU32 row = 0; row < numRows; row++) { - outRow = out->data.C64[row]; - realRow = real->data.F64[row]; - imagRow = imag->data.F64[row]; - - for (psU32 col = 0; col < numCols; col++) { - outRow[col] = realRow[col] + I * imagRow[col]; - } - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input psImage type, %s, is required to be either psF32 or psF64."), - typeStr); - psFree(out); - return NULL; - } - - - return out; -} - -psImage* psImageConjugate(psImage* out, const psImage* in) -{ - psElemType type; - psU32 numCols; - psU32 numRows; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numCols = in->numCols; - numRows = in->numRows; - - /* if not a complex image, this is logically just a image copy */ - if (!PS_IS_PSELEMTYPE_COMPLEX(type)) { - return psImageCopy(out, in, type); - } - - if (type == PS_TYPE_C32) { - psC32* outRow; - psC32* inRow; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_C32); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.C32[row]; - inRow = in->data.C32[row]; - - for (psU32 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 (psU32 row = 0; row < numRows; row++) { - outRow = out->data.C64[row]; - inRow = in->data.C64[row]; - - for (psU32 col = 0; col < numCols; col++) { - outRow[col] = creal(inRow[col]) - I * cimag(inRow[col]); - } - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input psImage type, %s, is required to be either psC32 or psC64."), - typeStr); - psFree(out); - return NULL; - } - - return out; -} - -psImage* psImagePowerSpectrum(psImage* out, const psImage* in) -{ - psElemType type; - psU32 numCols; - psU32 numRows; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - numCols = in->numCols; - numRows = in->numRows; - - if (type == PS_TYPE_C32) { - psF32* outRow; - psC32* inRow; - psF32 real; - psF32 imag; - psF32 fNumCols = numCols; - psF32 fNumRows = numRows; - psF32 numElementsSquared = fNumCols * fNumCols * fNumRows * fNumRows; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.F32[row]; - inRow = in->data.C32[row]; - - for (psU32 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; - psF64 numElementsSquared = numCols * numCols * numRows * numRows; - - out = psImageRecycle(out, numCols, numRows, PS_TYPE_F64); - for (psU32 row = 0; row < numRows; row++) { - outRow = out->data.F64[row]; - inRow = in->data.C64[row]; - - for (psU32 col = 0; col < numCols; col++) { - real = creal(inRow[col]); - imag = cimag(inRow[col]); - outRow[col] = (real * real + imag * imag) / numElementsSquared; - } - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input psImage type, %s, is required to be either psC32 or psC64."), - typeStr); - psFree(out); - return NULL; - } - - return out; - -} + *outImag = psImageRecycle(*outImag, numCols, numRows, PS_TYPE_F32); + targetImag = psMemIncrRefCounter(*outImag); + } + + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + // (a + bi) * (c + di) = (ac - bd) + (bc + ad)i + float real = in1Real->data.F32[y][x] * in2Real->data.F32[y][x] - + in1Imag->data.F32[y][x] * in2Imag->data.F32[y][x]; + float imag = in1Imag->data.F32[y][x] * in2Real->data.F32[y][x] + + in1Real->data.F32[y][x] * in2Imag->data.F32[y][x]; + targetReal->data.F32[y][x] = real; + targetImag->data.F32[y][x] = imag; + } + } + + if ((*outReal)->parent) { + *outReal = psImageCopy(*outReal, targetReal, PS_TYPE_F32); + } + if ((*outImag)->parent) { + *outImag = psImageCopy(*outImag, targetImag, PS_TYPE_F32); + } + + psFree(targetReal); + psFree(targetImag); + + return true; +}