Index: trunk/psLib/src/fft/psImageFFT.c =================================================================== --- trunk/psLib/src/fft/psImageFFT.c (revision 11701) +++ 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; +} Index: trunk/psLib/src/fft/psImageFFT.h =================================================================== --- trunk/psLib/src/fft/psImageFFT.h (revision 11701) +++ trunk/psLib/src/fft/psImageFFT.h (revision 11703) @@ -1,11 +1,12 @@ -/* @file psImageFFT.h - * @brief Contains FFT transform related functions for psImage - * - * @author Robert DeSonia, MHPCC - * - * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-01-23 22:47:22 $ - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ +/// @file psImageFFT.h +/// @brief Contains FFT transform related functions for psImage +/// +/// @author Paul Price, IfA +/// @author Robert DeSonia, MHPCC +/// +/// @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ +/// @date $Date: 2007-02-08 04:17:58 $ +/// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +/// #ifndef PS_IMAGE_FFT_H @@ -13,69 +14,50 @@ #include "psImage.h" -#include "psVectorFFT.h" // for psFFTFlags /// @addtogroup MathOps Mathematical Operations /// @{ -/** 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* image, ///< the psImage to apply transform to - psFFTFlags direction ///< the direction of the transform -); +/// Forward FFT of an image +/// +/// Applies a forward FFT (exponent -1), with the result returned in both real and imaginary parts. The FFT +/// is not normalised (a forward followed by a reverse is the original scaled by the total number of pixels). +/// The FFT takes advantage of the fact that the input is purely real; hence the output number of columns is +/// numCols/2 + 1 (with division rounding down). Only implemented for F32 input. +bool psImageForwardFFT(psImage **real, ///< Real part of FFT + psImage **imag, ///< Imaginary part of FFT + const psImage *in///< Input image (F32) + ); -/** 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 -); +/// Backward FFT of an image +/// +/// Applies a backward FFT (exponent +1) from the real and imaginary parts with the (purely) real result +/// returned. The FFT is not normalised (a forward followed by a reverse is the original scaled by the total +/// number of pixels). The FFT takes advantage of the fact that the output will be purely real; hence the +/// input number of columns is numCols/2 + 1 (with division rounding down); to manage the redundancy (is the +/// original number of columns even or odd?), we need the original number of columns (the number of columns of +/// the image that was input to psImageForwardFFT) to be provided. Only implemented for F32 input. +bool psImageBackwardFFT(psImage **out,///< Output image + const psImage *real, ///< Real input (F32) + const psImage *imag, ///< Imaginary input (F32) + int origCols ///< Original number of columns + ); -/** 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 -); +/// Power spectrum of an image +/// +/// Generates the power spectrum of an image. Only implemented for F32 input. +psImage* psImagePowerSpectrum(psImage *out, const psImage *in); -/** 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 -); +/// Multiply complex images +/// +/// The input images are the real and imaginary parts of each of two images. The real and imaginary parts +/// of the output image are returned. Only implemented for F32 input. +bool psImageComplexMultiply(psImage **outReal, ///< Real part of output + psimage **outImag, ///< Imaginary part of output + const psImage *in1Real, ///< Real part of input 1 + const psImage *in1Imag, ///< Imaginary part of input 1 + const psImage *in2Real, ///< Real part of input 2 + const psImage *in2Imag ///< Imaginary part of input 2 + ); /// @} Index: trunk/psLib/src/fft/psVectorFFT.c =================================================================== --- trunk/psLib/src/fft/psVectorFFT.c (revision 11701) +++ trunk/psLib/src/fft/psVectorFFT.c (revision 11703) @@ -3,8 +3,9 @@ * @brief Contains FFT transform related functions for psVector * + * @author Paul Price, IfA * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.36 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-02-06 21:36:09 $ + * @version $Revision: 1.37 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-02-08 04:17:58 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -12,5 +13,5 @@ #ifdef HAVE_CONFIG_H -# include "config.h" +#include "config.h" #endif @@ -21,401 +22,161 @@ #include -#include "psVectorFFT.h" +#include "psAssert.h" #include "psError.h" #include "psMemory.h" #include "psLogMsg.h" +#include "psConstants.h" +#include "psVectorFFT.h" + +#define FFTW_PLAN_RIGOR FFTW_ESTIMATE // How rigorous the FFTW planning is + +static psBool fftwWisdomImported = false; // Has the system wisdom been imported? + +bool psVectorForwardFFT(psVector **real, psVector **imag, const psVector *in) +{ + PS_ASSERT_VECTOR_NON_NULL(in, false); + PS_ASSERT_VECTOR_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; + } + + long num = in->n; // Number of elements + + // Do the FFT + fftwf_complex *out = fftwf_malloc(num * sizeof(fftwf_complex)); // Output data + fftwf_plan plan = fftwf_plan_dft_r2c_1d(num, in->data.F32, out, FFTW_PLAN_RIGOR); + fftwf_execute(plan); + fftwf_destroy_plan(plan); + + // Pull the real and imaginary parts out + *real = psVectorRecycle(*real, num/2 + 1, PS_TYPE_F32); + *imag = psVectorRecycle(*imag, num/2 + 1, PS_TYPE_F32); + for (int i = 0; i < num; i++) { +#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) + // C99 complex support + (*real)->data.F32[i] = creal(out[i]); + (*imag)->data.F32[i] = cimag(out[i]); +#else + // FFTW's backup complex support + (*real)->data.F32[i] = out[i][0]; + (*imag)->data.F32[i] = out[i][1]; +#endif + } + fftwf_free(out); + + return true; +} + +bool psVectorBackwardFFT(psVector **out, const psVector *real, const psVector *imag) +{ + PS_ASSERT_VECTOR_NON_NULL(real, false); + PS_ASSERT_VECTOR_TYPE(real, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_NON_NULL(imag, false); + PS_ASSERT_VECTOR_TYPE(imag, PS_TYPE_F32, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(real, imag, false); + PS_ASSERT_PTR_NON_NULL(out, false); + + // Make sure the system-level wisdom information is imported. + if (!fftwWisdomImported) { + fftwf_import_system_wisdom(); + fftwWisdomImported = true; + } + + long num = real->n; // Number of elements + + // Stuff the real and imaginary parts in + fftwf_complex *in = fftwf_malloc(num * sizeof(fftwf_complex)); // Input data + for (int i = 0; i < num; i++) { +#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) + // C99 complex support + in[i] = real->data.F32[i] + imag->data.F32[i] * I; +#else + // FFTW's backup complex support + in[i][0] = real->data.F32[i]; + in[i][1] = imag->data.F32[i]; +#endif + } + // Do the FFT + *out = psVectorRecycle(*out, 2 * num, PS_TYPE_F32); + fftwf_plan plan = fftwf_plan_dft_c2r_1d(2 * num, in, (*out)->data.F32, FFTW_PLAN_RIGOR); + fftwf_execute(plan); + fftwf_destroy_plan(plan); -#define P_FFTW_PLAN_RIGOR FFTW_ESTIMATE + fftwf_free(in); -static bool p_fftwWisdomImported = false; + return true; +} -psVector* psVectorFFT(psVector* out, const psVector* in, psFFTFlags direction) +psVector *psVectorPowerSpectrum(psVector* out, const psVector* in) { - psU32 numElements; - psElemType type; - fftwf_plan plan; + PS_ASSERT_VECTOR_NON_NULL(in, NULL); + PS_ASSERT_VECTOR_TYPE(in, PS_TYPE_F32, NULL); - /* got good image data? */ - if (in == NULL) { - psFree(out); + psVector *real = NULL; // Real component of FFT + psVector *imag = NULL; // Imaginary component of FFT + + if (!psVectorForwardFFT(&real, &imag, in)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform forward FFT."); return NULL; } - type = in->type.type; + int num = real->n; // Number of elements - /* make sure the system-level wisdom information is imported. */ - if (!p_fftwWisdomImported) { - fftwf_import_system_wisdom(); - p_fftwWisdomImported = true; + float norm = 1.0 / PS_SQR(in->n); + for (int i = 0; i < num; i++) { + // Power spectrum is the square of the complex modulus + real->data.F32[i] = norm * (PS_SQR(real->data.F32[i]) + PS_SQR(imag->data.F32[i])); + } + psFree(imag); + + return real; +} + +bool psVectorComplexMultiply(psVector **outReal, psVector **outImag, + const psVector *in1Real, const psVector *in1Imag, + const psVector *in2Real, const psVector *in2Imag) +{ + PS_ASSERT_VECTOR_NON_NULL(in1Real, false); + PS_ASSERT_VECTOR_NON_NULL(in1Imag, false); + PS_ASSERT_VECTOR_NON_NULL(in2Real, false); + PS_ASSERT_VECTOR_NON_NULL(in2Imag, false); + PS_ASSERT_VECTOR_TYPE(in1Real, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_TYPE(in1Imag, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_TYPE(in2Real, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_TYPE(in2Imag, PS_TYPE_F32, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(in1Imag, in1Real, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(in2Real, in1Real, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(in2Imag, in1Real, false); + PS_ASSERT_PTR_NON_NULL(outReal, false); + PS_ASSERT_PTR_NON_NULL(outImag, false); + if (*outReal) { + PS_ASSERT_VECTOR_NON_NULL(*outReal, false); + PS_ASSERT_VECTOR_NON_NULL(*outImag, false); + } else if (*outImag) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "If the output real part is provided, the output imaginary part must also be provided."); + return false; } - numElements = in->n; + int num = in1Real->n; // Number of elements - out = psVectorCopy(out, in, PS_TYPE_C32); - out->n = numElements; + *outReal = psVectorRecycle(*outReal, num, PS_TYPE_F32); + *outImag = psVectorRecycle(*outImag, num, PS_TYPE_F32); - if ((direction & PS_FFT_FORWARD) != 0) { - plan = fftwf_plan_dft_1d(numElements, - (fftwf_complex *) out->data.C32, - (fftwf_complex *) out->data.C32, FFTW_FORWARD, P_FFTW_PLAN_RIGOR); - } else if ((direction & PS_FFT_REVERSE) != 0) { - plan = fftwf_plan_dft_1d(numElements, - (fftwf_complex *) out->data.C32, - (fftwf_complex *) out->data.C32, FFTW_BACKWARD, P_FFTW_PLAN_RIGOR); - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Must specify the direction as either PS_FFT_FORWARD or PS_FFT_REVERSE.")); - psFree(out); - return NULL; + for (int i = 0; i < num; i++) { + // (a + bi) * (c + di) = (ac - bd) + (bc + ad)i + (*outReal)->data.F32[i] = in1Real->data.F32[i] * in2Real->data.F32[i] - + in1Imag->data.F32[i] * in2Imag->data.F32[i]; + (*outImag)->data.F32[i] = in1Imag->data.F32[i] * in2Real->data.F32[i] + + in1Real->data.F32[i] * in2Imag->data.F32[i]; } - /* check if a plan exists now */ - 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) != 0) { - for (psS32 i = 0; i < numElements; i++) { - out->data.F32[i] = out->data.C32[i]; - } - out->type.type = PS_TYPE_F32; - out->data.U8 = psRealloc(out->data.U8,PSELEMTYPE_SIZEOF(PS_TYPE_F32)*out->nalloc); - } - - return out; + return true; } - -psVector* psVectorReal(psVector* out, const psVector* in) -{ - psElemType type; - psU32 numElements; - - 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, numElements, type); - out->n = numElements; - memcpy(out->data.U8, in->data.U8, numElements * PSELEMTYPE_SIZEOF(type)); - return out; - } - - if (type == PS_TYPE_C32) { - psF32* outVec; - psC32* inVec = in->data.C32; - - out = psVectorRecycle(out, numElements, PS_TYPE_F32); - out->n = numElements; - outVec = out->data.F32; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = crealf(inVec[i]); - } - } else if (type == PS_TYPE_C64) { - psF64* outVec; - psC64* inVec = in->data.C64; - - out = psVectorRecycle(out, numElements, PS_TYPE_F64); - out->n = numElements; - outVec = out->data.F64; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = creal(inVec[i]); - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psVector type, %s, is not supported."), - typeStr); - psFree(out); - return NULL; - } - - return out; -} - -psVector* psVectorImaginary(psVector* out, const psVector* in) -{ - psElemType type; - psU32 numElements; - - 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, numElements, type); - out->n = numElements; - memset(out->data.U8, 0, PSELEMTYPE_SIZEOF(type) * numElements); - return out; - } - - if (type == PS_TYPE_C32) { - psF32* outVec; - psC32* inVec = in->data.C32; - - out = psVectorRecycle(out, numElements, PS_TYPE_F32); - out->n = numElements; - outVec = out->data.F32; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = cimagf(inVec[i]); - } - } else if (type == PS_TYPE_C64) { - psF64* outVec; - psC64* inVec = in->data.C64; - - out = psVectorRecycle(out, numElements, PS_TYPE_F64); - out->n = numElements; - outVec = out->data.F64; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = cimag(inVec[i]); - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psVector type, %s, is not supported."), - typeStr); - psFree(out); - return NULL; - } - - return out; -} - -psVector* psVectorComplex(psVector* out, const psVector* real, const psVector* imag) -{ - psElemType type; - psU32 numElements; - - 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) { - char* typeStrReal; - char* typeStrImag; - PS_TYPE_NAME(typeStrReal,type); - PS_TYPE_NAME(typeStrImag,imag->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Real psVector type, %s, and imaginary psVector type, %s, must be the same."), - typeStrReal,typeStrImag); - psFree(out); - return NULL; - } - - if (type == PS_TYPE_F32) { - psC32* outVec; - psF32* realVec = real->data.F32; - psF32* imagVec = imag->data.F32; - - out = psVectorRecycle(out, numElements, PS_TYPE_C32); - out->n = numElements; - outVec = out->data.C32; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = realVec[i] + I * imagVec[i]; - } - } else if (type == PS_TYPE_F64) { - psC64* outVec; - psF64* realVec = real->data.F64; - psF64* imagVec = imag->data.F64; - - out = psVectorRecycle(out, numElements, PS_TYPE_C64); - out->n = numElements; - outVec = out->data.C64; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = realVec[i] + I * imagVec[i]; - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input psVector type, %s, is required to be either psF32 or psF64."), - typeStr); - psFree(out); - return NULL; - } - - return out; -} - -psVector* psVectorConjugate(psVector* out, const psVector* in) -{ - psElemType type; - psU32 numElements; - - 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, numElements, type); - out->n = numElements; - memcpy(out->data.U8, in->data.U8, PSELEMTYPE_SIZEOF(type) * numElements); - return out; - } - - if (type == PS_TYPE_C32) { - psC32* outVec; - psC32* inVec = in->data.C32; - - out = psVectorRecycle(out, numElements, PS_TYPE_C32); - out->n = numElements; - outVec = out->data.C32; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = crealf(inVec[i]) - I * cimagf(inVec[i]); - } - } else if (type == PS_TYPE_C64) { - psC64* outVec; - psC64* inVec = in->data.C64; - - out = psVectorRecycle(out, numElements, PS_TYPE_C64); - out->n = numElements; - outVec = out->data.C64; - - for (psU32 i = 0; i < numElements; i++) { - outVec[i] = creal(inVec[i]) - I * cimag(inVec[i]); - } - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input psVector type, %s, is required to be either psC32 or psC64."), - typeStr); - psFree(out); - return NULL; - } - - return out; -} - -psVector* psVectorPowerSpectrum(psVector* out, const psVector* in) -{ - psElemType type; - psU32 outNumElements; - psU32 inNumElements; - psU32 inHalfNumElements; - psU32 inNumElementsSquared; - - if (in == NULL) { - psFree(out); - return NULL; - } - - type = in->type.type; - inNumElements = in->n; - inNumElementsSquared = inNumElements * inNumElements; - inHalfNumElements = inNumElements / 2; - outNumElements = inHalfNumElements + 1; - - if (type == PS_TYPE_C32) { - psF32* outVec; - psC32* inVec = in->data.C32; - psF32 inAbs1; - psF32 inAbs2; - - out = psVectorRecycle(out, outNumElements, PS_TYPE_F32); - out->n = outNumElements; - outVec = out->data.F32; - - // from ADD: P_0 = |C_0|^2/N^2 - inAbs1 = cabsf(inVec[0]); - outVec[0] = inAbs1 * inAbs1 / inNumElementsSquared; - - // from ADD: P_j = (|C_j|^2+|C_N-j|^2)/N^2, where j = 1,2,...,(N/2-1) - for (psU32 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 if (type == PS_TYPE_C64) { - psF64* outVec; - psC64* inVec = in->data.C64; - psF64 inAbs1; - psF64 inAbs2; - - out = psVectorRecycle(out, outNumElements, PS_TYPE_F64); - out->n = outNumElements; - outVec = out->data.F64; - - // from ADD: P_0 = |C_0|^2/N^2 - inAbs1 = cabs(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 (psU32 i = 1; i < inHalfNumElements; i++) { - inAbs1 = cabs(inVec[i]); - inAbs2 = cabs(inVec[inNumElements - i]); - outVec[i] = (inAbs1 * inAbs1 + inAbs2 * inAbs2) / inNumElementsSquared; - } - - // from ADD: P_N/2 = |C_N/2|^2/N^2 - inAbs1 = cabs(inVec[inHalfNumElements]); - outVec[inHalfNumElements] = inAbs1 * inAbs1 / inNumElementsSquared; - } else { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input psVector type, %s, is required to be either psC32 or psC64."), - typeStr); - psFree(out); - return NULL; - } - - return out; - -} Index: trunk/psLib/src/fft/psVectorFFT.h =================================================================== --- trunk/psLib/src/fft/psVectorFFT.h (revision 11701) +++ trunk/psLib/src/fft/psVectorFFT.h (revision 11703) @@ -1,11 +1,12 @@ -/* @file psVectorFFT.h - * @brief Contains FFT transform related functions for psVector - * - * @author Robert DeSonia, MHPCC - * - * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-01-23 22:47:22 $ - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ +/// @file psVectorFFT.h +/// @brief Contains FFT transform related functions for psVector +/// +/// @author Paul Price, IfA +/// @author Robert DeSonia, MHPCC +/// +/// @version $Revision: 1.20 $ $Name: not supported by cvs2svn $ +/// @date $Date: 2007-02-08 04:17:58 $ +/// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +/// #ifndef PS_VECTOR_FFT_H @@ -17,77 +18,47 @@ /// @{ -/** Specify direction of FFT */ -typedef enum { - /// psImageFFT/psVectorFFT should perform a forward FFT. - PS_FFT_FORWARD = 1, +/// Forward FFT of a vector +/// +/// Applies a forward FFT (exponent -1), with the result returned in both real and imaginary parts. The FFT +/// is not normalised (a forward followed by a reverse is the original scaled by the number of elements). The +/// FFT takes advantage of the fact that the input is purely real; hence the output size is N/2 + 1 (with +/// division rounding down). Only implemented for F32 input. +bool psVectorForwardFFT(psVector **real,///< Real part of FFT + psVector **imag,///< Imaginary part of FFT + const psVector *in ///< Input vector (F32) + ); - /// psImageFFT/psVectorFFT should perform a reverse FFT. - PS_FFT_REVERSE = 2, +/// Backward FFT of a vector +/// +/// Applies a backward FFT (exponent +1) from the real and imaginary parts with the (purely) real result +/// returned. The FFT is not normalised (a forward followed by a reverse is the original scaled by the number +/// of elements). The FFT takes advantage of the fact that the output will be purely real; hence the input +/// size is N/2 + 1 (with division rounding down); to manage the redundancy (is the original size even or +/// odd?), we need the original size (the size of the array that was input to psVectorForwardFFT) to be +/// provided. Only implemented for F32 input. +bool psVectorBackwardFFT(psVector **out,///< Output vector + const psVector *real, ///< Real input (F32) + const psVector *imag, ///< Imaginary input (F32) + int origN ///< Original number of elements + ); - /// psImageFFT/psVectorFFT should perform a reverse FFT with a real result. - PS_FFT_REAL_RESULT = 4 -} psFFTFlags; +/// Power spectrum of a vector +/// +/// Generates the power spectrum of a vector. Only implemented for F32 input. +psVector *psVectorPowerSpectrum(psVector *out, ///< Output power spectrum, or NULL + const psVector* in ///< Input vector (F32) + ); - -/** Forward and reverse FFT calculations. - * - * This takes as input the vector of interest (in) and the direction - * (direction), which is specified by an enumerated type psFftDirection. - * The input vector may be of type psF32 or psC32, the result is always - * psC32. If the input vector is psF32, the direction must be forward. - * - * @return psVector* the FFT transformation result - */ -psVector* psVectorFFT( - psVector* out, ///< a psVector to recycle. If NULL, a new psVector is made. - const psVector* in, ///< the vector to apply transform to - psFFTFlags direction ///< the direction of the transform -); - -/** extract the real portion of a complex vector - * - * @return psVector* real portion of the input vector. - */ -psVector* psVectorReal( - psVector* out, ///< a psVector to recycle. If NULL, a new psVector is made. - const psVector* in ///< the psVector to extract real portion from -); - -/** extract the imaginary portion of a complex vector - * - * @return psVector* imaginary portion of the input vector. - */ -psVector* psVectorImaginary( - psVector* out, ///< a psVector to recycle. If NULL, a new psVector is made. - const psVector* in ///< the psVector to extract imaginary portion from -); - -/** creates a complex vector from separate real and imaginary vectors - * - * @return psVector* resulting complex vector - */ -psVector* psVectorComplex( - psVector* out, ///< a psVector to recycle. If NULL, a new psVector is made. - const psVector* real, ///< the real vector - const psVector* imag ///< the imaginary vector -); - -/** computes the complex conjugate of a vector - * - * @return psVector* the complex conjugate of the 'in' vector - */ -psVector* psVectorConjugate( - psVector* out, ///< a psVector to recycle. If NULL, a new psVector is made. - const psVector* in ///< the psVector to compute conjugate of -); - -/** computes the power spectrum of a vector - * - * @return psVector* the power spectrum of the 'in' vector - */ -psVector* psVectorPowerSpectrum( - psVector* out, ///< a psVector to recycle. If NULL, a new psVector is made. - const psVector* in ///< the psVector to power spectrum of -); +/// Multiply complex vectors +/// +/// The input vectors are the real and imaginary parts of each of two vectors. The real and imaginary parts +/// of the output vector are returned. Only implemented for F32 input. +bool psVectorComplexMultiply(psVector **outReal, ///< Real part of output + psVector **outImag, ///< Imaginary part of output + const psVector *in1Real, ///< Real part of input 1 + const psVector *in1Imag, ///< Imaginary part of input 1 + const psVector *in2Real, ///< Real part of input 2 + const psVector *in2Imag ///< Imaginary part of input 2 + ); /// @} Index: trunk/psLib/src/imageops/psImageConvolve.c =================================================================== --- trunk/psLib/src/imageops/psImageConvolve.c (revision 11701) +++ trunk/psLib/src/imageops/psImageConvolve.c (revision 11703) @@ -1,42 +1,47 @@ -/* @file psImageConvolve.c - * - * @brief Contains FFT transform related functions for psImage. - * - * @author Robert DeSonia, MHPCC - * - * @version $Revision: 1.42 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-01-19 04:30:33 $ - * - * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - */ +/// @file psImageConvolve.c +/// +/// @brief Contains FFT transform related functions for psImage. +/// +/// @author Robert DeSonia, MHPCC +/// @author Paul Price, IfA +/// @author Eugene Magnier, IfA +/// +/// @version $Revision: 1.43 $ $Name: not supported by cvs2svn $ +/// @date $Date: 2007-02-08 04:17:58 $ +/// +/// Copyright 2004-2007 Institute for Astronomy, University of Hawaii +/// #ifdef HAVE_CONFIG_H -# include "config.h" +#include "config.h" #endif #include #include -#include "psImageConvolve.h" -#include "psImageFFT.h" -#include "psImageStructManip.h" -#include "psBinaryOp.h" +#include "psAbort.h" #include "psMemory.h" #include "psLogMsg.h" #include "psError.h" #include "psAssert.h" - - - -#define FOURIER_PADDING 32 /* padding amount in every side of the image for fourier convolution */ - -static void freeKernel(psKernel* ptr); - -psKernel* psKernelAlloc(int xMin, int xMax, int yMin, int yMax) -{ - psKernel* result; - psS32 numRows; - psS32 numCols; - - // following is explicitly spelled out in the SDRS as a requirement +#include "psScalar.h" +#include "psBinaryOp.h" +#include "psImageFFT.h" +#include "psImageStructManip.h" +#include "psImagePixelManip.h" + +#include "psImageConvolve.h" + +static void kernelFree(psKernel *kernel) +{ + if (kernel) { + psFree(kernel->image); + psFree(kernel->p_kernelRows); + } + return; +} + +psKernel *psKernelAlloc(int xMin, int xMax, int yMin, int yMax) +{ + // Check the inputs to make sure max > min; if not, switch. if (yMin > yMax) { psLogMsg(__func__, PS_LOG_WARN, @@ -44,5 +49,5 @@ yMin, yMax); - psS32 temp = yMin; + int temp = yMin; yMin = yMax; yMax = temp; @@ -55,40 +60,33 @@ xMin, xMax); - psS32 temp = xMin; + int temp = xMin; xMin = xMax; xMax = temp; } - numRows = yMax - yMin + 1; - numCols = xMax - xMin + 1; - - result = psAlloc(sizeof(psKernel)); - result->xMin = xMin; - result->xMax = xMax; - result->yMin = yMin; - result->yMax = yMax; - result->image = psImageAlloc(numCols,numRows,PS_TYPE_KERNEL); - memset((result->image->p_rawDataBuffer),0,numCols*numRows*PSELEMTYPE_SIZEOF(PS_TYPE_KERNEL)); - result->p_kernelRows = psAlloc(sizeof(float*)*numRows); - - float** kernelRows = result->p_kernelRows; - float** imageRows = result->image->data.PS_TYPE_KERNEL_DATA; - for (psS32 i = 0; i < numRows; i++) { - kernelRows[i] = imageRows[i] - xMin; - } - result->kernel = kernelRows - yMin; - - psMemSetDeallocator(result,(psFreeFunc)freeKernel); - - return result; -} - -void freeKernel(psKernel* ptr) -{ - if (ptr != NULL) { - psFree(ptr->image); - psFree(ptr->p_kernelRows); - } -} + int numRows = yMax - yMin + 1; // Number of rows for kernel image + int numCols = xMax - xMin + 1; // Number of columns for kernel image + + psKernel *kernel = psAlloc(sizeof(psKernel)); // The kernel, to be returned + psMemSetDeallocator(kernel,(psFreeFunc)kernelFree); + + kernel->xMin = xMin; + kernel->xMax = xMax; + kernel->yMin = yMin; + kernel->yMax = yMax; + kernel->image = psImageAlloc(numCols, numRows, PS_TYPE_KERNEL); + psImageInit(kernel->image, 0.0); + + // Set up indirections, so we can refer to kernel->kernel[-1][-3] for the (-1,-1) element, instead of + // kernel->image[kernel->yMax - kernel->yMin + 1][kernel->yMax - kernel->yMin + 3] (yuk!). + kernel->p_kernelRows = psAlloc(sizeof(float*)*numRows); + for (int i = 0; i < numRows; i++) { + kernel->p_kernelRows[i] = kernel->image->data.PS_TYPE_KERNEL_DATA[i] - xMin; + } + kernel->kernel = kernel->p_kernelRows - yMin; + + return kernel; +} + @@ -96,397 +94,264 @@ { PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc)freeKernel ); -} - - -psKernel* psKernelGenerate(const psVector* tShifts, - const psVector* xShifts, - const psVector* yShifts, - bool relative) -{ - psS32 lastX; - psS32 lastY; - psS32 lastT; - psS32 x; - psS32 y; - psS32 t; - psS32 xMin = 0; - psS32 xMax = 0; - psS32 yMin = 0; - psS32 yMax = 0; - psS32 length = 0; - float normalizeTime = 1.0; // fraction of total time for each shift clock - psKernel* result = NULL; - float** kernel = NULL; - - // got non-NULL vectors? - if (tShifts == NULL || xShifts == NULL || yShifts == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("Specified shift vectors can not be NULL.")); + return ( psMemGetDeallocator(ptr) == (psFreeFunc)kernelFree ); +} + + +psKernel *psKernelGenerate(const psVector *tShifts, + const psVector *xShifts, + const psVector *yShifts, + bool tRelative, + bool xyRelative) +{ + PS_ASSERT_VECTOR_NON_NULL(tShifts, NULL); + PS_ASSERT_VECTOR_NON_NULL(xShifts, NULL); + PS_ASSERT_VECTOR_NON_NULL(yShifts, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(tShifts, xShifts, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(tShifts, yShifts, NULL); + PS_ASSERT_VECTOR_TYPE(tShifts, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE(xShifts, PS_TYPE_S32, NULL); + PS_ASSERT_VECTOR_TYPE(yShifts, PS_TYPE_S32, NULL); + + // If there are no shifts, the kernel is just a 1 at 0,0 + long num = tShifts->n; // Number of shifts + if (num == 0) { + psKernel *kernel = psKernelAlloc(0,0,0,0); + kernel->kernel[0][0] = 1; + return kernel; + } + + // Get dimensions and scaling + int xMin, xMax, yMin, yMax; // Range of values for kernel + int xLast, yLast; // Last location, for relative shifts + xLast = xMin = xMax = xShifts->data.S32[0]; + yLast = yMin = yMax = yShifts->data.S32[0]; + float tSum = tShifts->data.F32[0]; // Sum of the times + for (long i = 1; i < num; i++) { + int x = xShifts->data.S32[i]; // x position in kernel + int y = yShifts->data.S32[i]; // y position in kernel + if (xyRelative) { + x += xLast; + y += yLast; + xLast = x; + yLast = y; + } + if (x < xMin) { + xMin = x; + } + if (x > xMax) { + xMax = x; + } + if (y < yMin) { + yMin = y; + } + if (y > yMax) { + yMax = y; + } + + if (tRelative) { + tSum += tShifts->data.F32[i]; + } + } + + if (!tRelative) { + // Then the total time is simply the final value + // NB: We assume the counter starts at zero! + tSum = tShifts->data.F32[tShifts->n - 1]; + } + + // One more pass through to set the kernel + psKernel *kernel = psKernelAlloc(xMin, xMax, yMin, yMax); // The kernel + xLast = xShifts->data.S32[0]; + yLast = yShifts->data.S32[0]; + float tLast = 0.0; // Last value for t + for (int i = 0; i < num; i++) { + int x = xShifts->data.S32[i]; // x position in kernel + int y = yShifts->data.S32[i]; // y position in kernel + if (xyRelative) { + x += xLast; + y += yLast; + xLast = x; + yLast = y; + } + float t = tShifts->data.F32[i]; + if (tRelative) { + t -= tLast; + tLast = tShifts->data.F32[i]; + } + + kernel->kernel[y][x] += t; + } + + // Normalise the kernel by the total time (kernel sum should be unity) + psBinaryOp(kernel->image, kernel->image, "*", psScalarAlloc(1.0 / tSum, PS_TYPE_F32)); + + return kernel; +} + +psImage *psImageConvolveDirect(psImage *out, + const psImage *in, + const psKernel *kernel) +{ + PS_ASSERT_IMAGE_NON_NULL(in, NULL); + PS_ASSERT_IMAGE_TYPE_F32_OR_F64(in, NULL); + PS_ASSERT_PTR_NON_NULL(kernel, NULL); + PS_ASSERT_PTR_NON_NULL(kernel->kernel, NULL); + + // Pull out kernel parameters, for convenience + int xMin = kernel->xMin; + int xMax = kernel->xMax; + int yMin = kernel->yMin; + int yMax = kernel->yMax; + float **kernelData = kernel->kernel; + + int numRows = in->numRows; // Number of rows + int numCols = in->numCols; // Number of columns + +#define SPATIAL_CONVOLVE_CASE(TYPE) \ + case PS_TYPE_##TYPE: { \ + ps##TYPE **inData = in->data.TYPE; \ + out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \ + for (int row = 0; row < numRows; row++) { \ + ps##TYPE *outRow = out->data.TYPE[row]; \ + for (int col = 0; col < numCols; col++) { \ + ps##TYPE pixel = 0.0; \ + for (int kRow = PS_MAX(yMin, -row); kRow <= PS_MIN(yMax, numRows - row - 1); kRow++) { \ + for (int kCol = PS_MAX(xMin, -col); kCol <= PS_MIN(xMax, numCols - col - 1); kCol++) { \ + pixel += kernelData[kRow][kCol] * inData[row + kRow][col + kCol]; \ + } \ + } \ + outRow[col] = pixel; \ + } \ + } \ + } \ + break; + + switch (in->type.type) { + SPATIAL_CONVOLVE_CASE(F32); + SPATIAL_CONVOLVE_CASE(F64); + default: + psAbort(PS_FILE_LINE, "Should never get here: bad type that was asserted on previously."); + } + + return out; +} + +psImage *psImageConvolveFFT(psImage *out, + const psImage *in, + const psKernel *kernel, + float pad) +{ + PS_ASSERT_IMAGE_NON_NULL(in, NULL); + PS_ASSERT_IMAGE_TYPE(in, PS_TYPE_F32, NULL); + PS_ASSERT_PTR_NON_NULL(kernel, NULL); + PS_ASSERT_IMAGE_NON_NULL(kernel->image, NULL); + + // Pull out kernel parameters, for convenience + int xMin = kernel->xMin; + int xMax = kernel->xMax; + int yMin = kernel->yMin; + int yMax = kernel->yMax; + + int numRows = in->numRows; // Number of rows in input image + int numCols = in->numCols; // Number of columns in input image + + // Need to pad the input image to protect from wrap-around effects + if (xMax - xMin > numCols || yMax - yMin > numRows) { + // Cannot pad the image if the kernel is larger. + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + _("Kernel cannot extend further than input image size (%dx%d vs %dx%d)."), + xMax, yMax, numCols, numRows); return NULL; } - - // types match? - if (xShifts->type.type != yShifts->type.type || - tShifts->type.type != xShifts->type.type) { - char* typeXStr; - char* typeYStr; - char* typeTStr; - PS_TYPE_NAME(typeXStr,xShifts->type.type); - PS_TYPE_NAME(typeYStr,yShifts->type.type); - PS_TYPE_NAME(typeTStr,tShifts->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Input t-, x-, and y-shift vector types (%s/%s/%s) must match."), - typeTStr, typeXStr, typeYStr); + int paddedCols = numCols + PS_MAX(-xMin, xMax); // Number of columns in padded image + int paddedRows = numRows + PS_MAX(-yMin, yMax); // Number of rows in padded image + + // Generate padded image + psImage *paddedImage = psImageAlloc(paddedCols,paddedRows,in->type.type); // Padded input image + psImageOverlaySection(paddedImage, in, 0, 0, "="); + for (int y = 0; y < numRows; y++) { + for (int x = numCols; x < paddedCols; x++) { + paddedImage->data.F32[y][x] = pad; + } + } + for (int y = numRows; y < paddedRows; y++) { + for (int x = 0; x < paddedCols; x++) { + paddedImage->data.F32[y][x] = pad; + } + } + + // Result of FFT + psImage *inRealFFT = NULL, *inImagFFT = NULL; + if (!psImageForwardFFT(&inRealFFT, &inImagFFT, paddedImage)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to fourier transform input image.")); + psFree(paddedImage); return NULL; } - - // sizes match? - length = xShifts->n; - if (length != yShifts->n || - length != tShifts->n) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - "Shift vectors can not be of different sizes."); + psFree(paddedImage); + + // Generate padded kernel image + psImage *paddedKernel = psImageAlloc(paddedCols, paddedRows, PS_TYPE_F32); + psImageInit(paddedKernel, 0.0); + for (int y = PS_MIN(-1, yMin); y <= PS_MIN(-1, yMax); y++) { + // y is negative + for (int x = PS_MIN(-1, xMin); x <= PS_MIN(-1, xMax); x++) { + // x is negative + paddedKernel->data.F32[paddedRows + y][paddedCols + x] = kernel->kernel[y][x]; + } + for (int x = PS_MAX(0, xMin); x <= PS_MAX(0, xMax); x++) { + // x is positive + paddedKernel->data.F32[paddedRows + y][x] = kernel->kernel[y][x]; + } + } + for (int y = PS_MAX(0, yMin); y <= PS_MAX(0, yMax); y++) { + // y is positive + for (int x = PS_MIN(-1, xMin); x <= PS_MIN(-1, xMax); x++) { + // x is negative + paddedKernel->data.F32[y][paddedCols + x] = kernel->kernel[y][x]; + } + for (int x = PS_MAX(0, xMin); x <= PS_MAX(0, xMax); x++) { + // x is positive + paddedKernel->data.F32[y][x] = kernel->kernel[y][x]; + } + } + + psImage *kernelRealFFT = NULL, *kernelImagFFT = NULL; + if (!psImageForwardFFT(&kernelRealFFT, &kernelImagFFT, paddedKernel)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to fourier transform kernel.")); + psFree(inRealFFT); + psFree(inImagFFT); + psFree(paddedKernel); return NULL; } - - // if no shifts, the kernel is just a 1 at 0,0 - if (length < 1) { - result = psKernelAlloc(0,0,0,0); - result->kernel[0][0] = 1; - return result; - } - - #define KERNEL_GENERATE_CASE(TYPE) \ -case PS_TYPE_##TYPE: { \ - ps##TYPE *tShiftData = tShifts->data.TYPE; \ - ps##TYPE *xShiftData = xShifts->data.TYPE; \ - ps##TYPE *yShiftData = yShifts->data.TYPE; \ - lastX = xShiftData[length-1]; \ - lastY = yShiftData[length-1]; \ - lastT = tShiftData[length-1]; \ - \ - for (int lcv = 0; lcv < length; lcv++) { \ - x = lastX - xShiftData[lcv]; \ - y = lastY - yShiftData[lcv]; \ - \ - if (x < xMin) { \ - xMin = x; \ - } else if (x > xMax) { \ - xMax = x; \ - } \ - if (y < yMin) { \ - yMin = y; \ - } else if (y > yMax) { \ - yMax = y; \ - } \ - } \ - \ - normalizeTime = 1.0 / (float)(tShiftData[length-1]); \ - result = psKernelAlloc(xMin,xMax,yMin,yMax); \ - kernel = result->kernel; \ - \ - psS32 prevT = 0; \ - for (int i = 0; i < length; i++) { \ - t = tShiftData[i] - prevT; \ - x = lastX - xShiftData[i]; \ - y = lastY - yShiftData[i]; \ - \ - kernel[y][x] += (float)t / (float)lastT; \ - prevT = tShiftData[i]; \ - } \ - break; \ - } - - #define RELATIVE_KERNEL_GENERATE_CASE(TYPE) \ -case PS_TYPE_##TYPE: { \ - ps##TYPE *tShiftData = tShifts->data.TYPE; \ - ps##TYPE *xShiftData = xShifts->data.TYPE; \ - ps##TYPE *yShiftData = yShifts->data.TYPE; \ - \ - x = 0; \ - y = 0; \ - t = 0; \ - \ - for (int lcv = length-1; lcv >= 0; lcv--) { \ - t += tShiftData[lcv]; \ - \ - if (x < xMin) { \ - xMin = x; \ - } else if (x > xMax) { \ - xMax = x; \ - } \ - if (y < yMin) { \ - yMin = y; \ - } else if (y > yMax) { \ - yMax = y; \ - } \ - x -= xShiftData[lcv]; \ - y -= yShiftData[lcv]; \ - \ - } \ - result = psKernelAlloc(xMin,xMax,yMin,yMax); \ - kernel = result->kernel; \ - \ - normalizeTime = 1.0 / (float)t; \ - x = 0; \ - y = 0; \ - for (psS32 i = length-1; i >= 0; i--) { \ - kernel[y][x] += (float)(tShiftData[i]) * normalizeTime; \ - x -= xShiftData[i]; \ - y -= yShiftData[i]; \ - \ - } \ - break; \ - } - - if (relative) { - switch (xShifts->type.type) { - RELATIVE_KERNEL_GENERATE_CASE(U8); - RELATIVE_KERNEL_GENERATE_CASE(U16); - RELATIVE_KERNEL_GENERATE_CASE(U32); - RELATIVE_KERNEL_GENERATE_CASE(U64); - RELATIVE_KERNEL_GENERATE_CASE(S8); - RELATIVE_KERNEL_GENERATE_CASE(S16); - RELATIVE_KERNEL_GENERATE_CASE(S32); - RELATIVE_KERNEL_GENERATE_CASE(S64); - RELATIVE_KERNEL_GENERATE_CASE(F32); - RELATIVE_KERNEL_GENERATE_CASE(F64); - RELATIVE_KERNEL_GENERATE_CASE(C32); - RELATIVE_KERNEL_GENERATE_CASE(C64); - - default: { - char* typeStr; - PS_TYPE_NAME(typeStr,xShifts->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psImage type, %s, is not supported."), - typeStr); - } - } - } else { - switch (xShifts->type.type) { - KERNEL_GENERATE_CASE(U8); - KERNEL_GENERATE_CASE(U16); - KERNEL_GENERATE_CASE(U32); - KERNEL_GENERATE_CASE(U64); - KERNEL_GENERATE_CASE(S8); - KERNEL_GENERATE_CASE(S16); - KERNEL_GENERATE_CASE(S32); - KERNEL_GENERATE_CASE(S64); - KERNEL_GENERATE_CASE(F32); - KERNEL_GENERATE_CASE(F64); - KERNEL_GENERATE_CASE(C32); - KERNEL_GENERATE_CASE(C64); - - default: { - char* typeStr; - PS_TYPE_NAME(typeStr,xShifts->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psImage type, %s, is not supported."), - typeStr); - } - } - } - - return result; -} - -psImage* psImageConvolve(psImage* out, - const psImage* in, - const psKernel* kernel, - bool direct) -{ - if (in == NULL) { - psFree(out); + psFree(paddedKernel); + + // Convolution in fourier domain is just a pixel-wise multiplication + if (!psImageComplexMultiply(&inRealFFT, &inImagFFT, inRealFFT, inImagFFT, kernelRealFFT, kernelImagFFT)) { + psError(PS_ERR_UNKNOWN, false, _("Unable to multiply fourier transformts.")); + psFree(inRealFFT); + psFree(inImagFFT); + psFree(kernelRealFFT); + psFree(kernelImagFFT); return NULL; } - - if (kernel == NULL) { - psError(PS_ERR_BAD_PARAMETER_NULL, true, - _("Specified psKernel can not be NULL.")); - psFree(out); + psFree(kernelRealFFT); + psFree(kernelImagFFT); + + psImage *paddedConvolved = NULL; // Padded convolved image + if (!psImageBackwardFFT(&paddedConvolved, inRealFFT, inImagFFT, paddedCols)) { + psError(PS_ERR_UNKNOWN, false, _("Failed to invert fourier transform of convolution image.")); + psFree(inRealFFT); + psFree(inImagFFT); return NULL; } - psS32 xMin = kernel->xMin; - psS32 xMax = kernel->xMax; - psS32 yMin = kernel->yMin; - psS32 yMax = kernel->yMax; - float** kData = kernel->kernel; - - // make the output image to the proper size and type - psS32 numRows = in->numRows; - psS32 numCols = in->numCols; - - - - if (direct) { - // spatial convolution - - #define SPATIAL_CONVOLVE_CASE(TYPE) \ - case PS_TYPE_##TYPE: { \ - ps##TYPE** inData = in->data.TYPE; \ - out = psImageRecycle(out, numCols, numRows, PS_TYPE_##TYPE); \ - for (psS32 row=0;rowdata.TYPE[row]; \ - for (psS32 col=0;col= 0 && row-kRow < numRows) { \ - for (psS32 kCol = xMin; kCol < xMax; kCol++) { \ - if (col-kCol >= 0 && col-kCol < numCols) { \ - pixel += kData[kRow][kCol] * inData[row-kRow][col-kCol]; \ - } \ - } \ - } \ - } \ - outRow[col] = pixel; \ - } \ - } \ - } \ - break; - - switch (in->type.type) { - SPATIAL_CONVOLVE_CASE(U8) - SPATIAL_CONVOLVE_CASE(U16) - SPATIAL_CONVOLVE_CASE(U32) - SPATIAL_CONVOLVE_CASE(U64) - SPATIAL_CONVOLVE_CASE(S8) - SPATIAL_CONVOLVE_CASE(S16) - SPATIAL_CONVOLVE_CASE(S32) - SPATIAL_CONVOLVE_CASE(S64) - SPATIAL_CONVOLVE_CASE(F32) - SPATIAL_CONVOLVE_CASE(F64) - SPATIAL_CONVOLVE_CASE(C32) - SPATIAL_CONVOLVE_CASE(C64) - - default: { - char* typeStr; - PS_TYPE_NAME(typeStr,in->type.type); - psError(PS_ERR_BAD_PARAMETER_TYPE, true, - _("Specified psImage type, %s, is not supported."), - typeStr); - psFree(out); - return NULL; - - } - } - - - } else { - // fourier convolution - psS32 paddedCols = numCols+2*FOURIER_PADDING; - psS32 paddedRows = numRows+2*FOURIER_PADDING; - - // check to see if kernel is smaller, otherwise padding it up will fail. - psS32 kRows = kernel->image->numRows; - psS32 kCols = kernel->image->numCols; - if (kRows >= numRows || kCols >= numCols) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - _("Specified psKernel size, %dx%d, can not be larger than input psImage size, %dx%d."), - kCols,kRows, - numCols, numRows); - psFree(out); - return NULL; - } - - // pad the image - psImage* paddedImage = psImageAlloc(paddedCols,paddedRows,in->type.type); - psS32 elementSize = PSELEMTYPE_SIZEOF(in->type.type); - - // zero out padded area on top and bottom - memset(paddedImage->data.U8[0],0,FOURIER_PADDING*paddedCols*elementSize); - memset(paddedImage->data.U8[FOURIER_PADDING+numRows-1],0,FOURIER_PADDING*paddedCols*elementSize); - - // fill in the image-containing rows. - psS32 sidePaddingSize = FOURIER_PADDING*elementSize; - psS32 imageRowSize = numCols*elementSize; - psU8* paddedData = paddedImage->data.U8[FOURIER_PADDING]; - for (psS32 row=0;rowdata.U8[row],imageRowSize); - paddedData += imageRowSize; - // zero out padded area on right edge. - memset(paddedData,0,sidePaddingSize); - paddedData += sidePaddingSize; - } - - // pad the kernel to the same size of paddedImage - psImage* paddedKernel = psImageAlloc(paddedCols,paddedRows,PS_TYPE_KERNEL); - memset(paddedKernel->data.U8[0],0,sizeof(float)*numCols*numRows); // zero-out image - psS32 yMax = kernel->yMax; - psS32 xMax = kernel->xMax; - for (psS32 row = kernel->yMin; row <= yMax;row++) { - psS32 padRow = row; - if (padRow < 0) { - padRow += paddedRows; - } - float* padData = paddedKernel->data.PS_TYPE_KERNEL_DATA[padRow]; - float* kernelRow = kernel->kernel[row]; - for (psS32 col = kernel->xMin; col <= xMax; col++) { - if (col < 0) { - padData[col+paddedCols] = kernelRow[col]; - } else { - padData[col] = kernelRow[col]; - } - } - } - - psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD); - if (kernelFourier == NULL) { - psError(PS_ERR_UNKNOWN, false, - _("Failed to perform a fourier transform of kernel.")); - psFree(out); - return NULL; - } - - psImage* inFourier = psImageFFT(NULL, paddedImage, PS_FFT_FORWARD); - if (inFourier == NULL) { - psError(PS_ERR_UNKNOWN, false, - _("Failed to perform a fourier transform of input image.")); - psFree(out); - return NULL; - } - - // convolution in fourier domain is just a pixel-wise multiplication - for (int row = 0; row < paddedRows; row++) { - psC32* inRow = inFourier->data.C32[row]; - psC32* kRow = kernelFourier->data.C32[row]; - for (int col = 0; col < paddedCols; col++) { - inRow[col] *= kRow[col]; - } - } - - psImage* complexOut = psImageFFT(NULL, inFourier, - PS_FFT_REVERSE); - if (complexOut == NULL) { - psError(PS_ERR_UNKNOWN, false, - _("Failed to perform a fourier transform of input image.")); - psFree(out); - return NULL; - } - - // subset out the padded area now. - psImage* complexOutSansPad = psImageSubset(complexOut, - psRegionSet(FOURIER_PADDING, FOURIER_PADDING+numCols,FOURIER_PADDING,FOURIER_PADDING+numRows)); - - out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); - float factor = 1.0f/(float)paddedCols/(float)paddedRows; - for (psS32 row = 0; row < numRows; row++) { - psF32* outRow = out->data.F32[row]; - psC32* resultRow = complexOutSansPad->data.C32[row]; - for (psS32 col = 0; col < numCols; col++) { - outRow[col] = crealf(resultRow[col])*factor; - } - } - - psFree(complexOut); // frees complexOutSansPad, as it is a child. - psFree(kernelFourier); - psFree(inFourier); - psFree(paddedImage); - psFree(paddedKernel); - - } + psFree(inRealFFT); + psFree(inImagFFT); + + // Trim off the padding, then renormalise (which also does a copy, so there's no parent for the output) + psImage *convolved = psImageSubset(paddedConvolved, psRegionSet(0, numCols, 0, numRows)); + out = (psImage*)psBinaryOp(out, convolved, "*", + psScalarAlloc(1.0 / paddedCols / paddedRows, PS_TYPE_F32)); + psFree(convolved); + psFree(paddedConvolved); return out; @@ -661,5 +526,5 @@ IMAGESMOOTH_CASE(F64); default: { - char* typeStr; + char *typeStr; PS_TYPE_NAME(typeStr,image->type.type); psError(PS_ERR_BAD_PARAMETER_TYPE, true, Index: trunk/psLib/src/imageops/psImageConvolve.h =================================================================== --- trunk/psLib/src/imageops/psImageConvolve.h (revision 11701) +++ trunk/psLib/src/imageops/psImageConvolve.h (revision 11703) @@ -5,6 +5,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.15 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-01-23 22:47:23 $ + * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-02-08 04:17:58 $ * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii */ @@ -20,55 +20,48 @@ #include "psType.h" -#define PS_TYPE_KERNEL PS_TYPE_F32 /**< the data member to use for kernel image */ -#define PS_TYPE_KERNEL_DATA F32 /**< the data member to use for kernel image */ -#define PS_TYPE_KERNEL_NAME "psF32" /**< the data type for kernel as a string */ +#define PS_TYPE_KERNEL PS_TYPE_F32 ///< the data member to use for kernel image */ +#define PS_TYPE_KERNEL_DATA F32 ///< the data member to use for kernel image */ +#define PS_TYPE_KERNEL_NAME "psF32" ///< the data type for kernel as a string */ -/** Kernel Type - * - * A floating-point data type used for storing kernel data. - * - */ -//typedef float psKernelType; - -/** A convolution kernel */ +/// A convolution kernel typedef struct { - psImage* image; ///< Kernel data, in the form of an image + psImage *image; ///< Kernel data, in the form of an image int xMin; ///< Most negative x index int yMin; ///< Most negative y index int xMax; ///< Most positive x index int yMax; ///< Most positive y index - float** kernel; ///< Pointer to the kernel data - float** p_kernelRows; ///< Pointer to the rows of the kernel data; not intended for user use. + float **kernel; ///< Pointer to the kernel data + float **p_kernelRows; ///< Pointer to the rows of the kernel data; not intended for user use. } psKernel; -/** Allocates a convolution kernel of the given range - * - * In order to perform a convolution, we need to define the convolution - * kernel. We need a more general object than a psImage so that we can - * incorporate the offset from the (0, 0) pixel to the (0, 0) value of the - * kernel. It might be convenient to allow both positive and negative - * indices to convey the positive and negative shifts. One might consider - * setting the x0 and y0 members of a psImage to the appropriate offsets, - * but this is not the purpose of these members, and doing so may affect the - * behavior of other psImage operations. - * - * This construction allows the kernel member to use negative indices, while - * preserving the location of psMemBlocks relative to allocated memory. - * - * The maximum extent of the kernel shifts shall be defined by the xMin, - * xMax, yMin and yMax members. Note that xMin and yMin, under normal - * circumstances, should be negative numbers. That is, - * myKernel->kernel[-3][-2] may be defined if yMin and xMin are equal to or - * more negative than -3 and -2, respectively. - * - * In the event that one of the minimum values is greater than the - * corresponding maximum value, the function shall generate a warning, and - * the offending values shall be exchanged. - * - * @return psKernel* A new kernel object - */ -psKernel* psKernelAlloc( +/// Allocates a convolution kernel of the given range +/// +/// In order to perform a convolution, we need to define the convolution +/// kernel. We need a more general object than a psImage so that we can +/// incorporate the offset from the (0, 0) pixel to the (0, 0) value of the +/// kernel. It might be convenient to allow both positive and negative +/// indices to convey the positive and negative shifts. One might consider +/// setting the x0 and y0 members of a psImage to the appropriate offsets, +/// but this is not the purpose of these members, and doing so may affect the +/// behavior of other psImage operations. +/// +/// This construction allows the kernel member to use negative indices, while +/// preserving the location of psMemBlocks relative to allocated memory. +/// +/// The maximum extent of the kernel shifts shall be defined by the xMin, +/// xMax, yMin and yMax members. Note that xMin and yMin, under normal +/// circumstances, should be negative numbers. That is, +/// myKernel->kernel[-3][-2] may be defined if yMin and xMin are equal to or +/// more negative than -3 and -2, respectively. +/// +/// In the event that one of the minimum values is greater than the +/// corresponding maximum value, the function shall generate a warning, and +/// the offending values shall be exchanged. +/// +/// @return psKernel* A new kernel object +/// +psKernel *psKernelAlloc( int xMin, ///< Most negative x index int xMax, ///< Most positive x index @@ -77,10 +70,10 @@ ); -/** Checks the type of a particular pointer. - * - * Uses the appropriate deallocation function in psMemBlock to check the ptr datatype. - * - * @return bool: True if the pointer matches a psKernel structure, false otherwise. - */ +/// Checks the type of a particular pointer. +/// +/// Uses the appropriate deallocation function in psMemBlock to check the ptr datatype. +/// +/// @return bool: True if the pointer matches a psKernel structure, false otherwise. +/// bool psMemCheckKernel( psPtr ptr ///< the pointer whose type to check @@ -88,60 +81,58 @@ -/** Generates a kernel given a list of shift values - * - * Given a list of values (e.g., shifts made in the course of OT guiding), - * psKernelGenerate shall return the appropriate kernel. The vectors xShifts - * and yShifts, which are a list of shifts relative to some starting point, - * will be supplied by the user. The elements of the vectors should be of an - * integer type; otherwise the values shall be truncated to integers. The - * output kernel shall be normalized such that the sum over the kernel is - * unity. - * - * If the vectors are not of the same number of elements, then the function - * shall generate a warning shall be generated, following which, the longer - * vector trimmed to the length of the shorter, and the function shall continue. - * - * @return psKernel* new Kernel object - */ -psKernel* psKernelGenerate( - const psVector* tShifts, ///< list of time shifts - const psVector* xShifts, ///< list of x-axis shifts - const psVector* yShifts, ///< list of y-axis shifts - bool relative - /**< specifies the starting point for the shifts; true=relative to previous shift - * false = relative to some other starting point. */ +/// Generates a kernel given a list of shift values +/// +/// Given a list of values (e.g., shifts made in the course of OT guiding), +/// psKernelGenerate shall return the appropriate kernel. The vectors xShifts +/// and yShifts, which are a list of shifts relative to some starting point, +/// will be supplied by the user. The elements of the vectors should be of an +/// integer type; otherwise the values shall be truncated to integers. The +/// output kernel shall be normalized such that the sum over the kernel is +/// unity. +/// +/// If the vectors are not of the same number of elements, then the function +/// shall generate a warning shall be generated, following which, the longer +/// vector trimmed to the length of the shorter, and the function shall continue. +/// +/// @return psKernel* new Kernel object +/// +psKernel *psKernelGenerate( + const psVector *tShifts, ///< list of time shifts (F32) + const psVector *xShifts, ///< list of x-axis shifts (S32) + const psVector *yShifts, ///< list of y-axis shifts (S32) + bool tRelative, ///< Are times relative (durations) or absolute? + bool xyRelative ///< Are x,y positions relative (shifts) or absolute? ); -/** convolve an image with a kernel - * - * Given an input image and the convolution kernel, psImageConvolve shall - * convolve the input image, in, with the kernel, kernel and return the - * convolved image, out. - * - * Two methods shall be available for the convolution: if direct is true, - * then the convolution shall be performed in real space (appropriate for - * small kernels); otherwise, the convolution shall be performed using Fast - * Fourier Transforms (FFTs; appropriate for larger kernels). The latter - * option involves padding the input image, copying the kernel into an image - * of the same size as the padded input image, performing an FFT on each, - * multiplying the FFTs, and performing an inverse FFT before trimming the - * image back to the original size. - * - * @return psImage* resulting image - */ -psImage* psImageConvolve( - psImage* out, ///< a psImage to recycle. If NULL, a new psImage is made. - const psImage* in, ///< the psImage to convolve - const psKernel* kernel, ///< kernel to colvolve with - bool direct ///< specifies method, true=direct convolution, false=fourier +/// Convolve an image with a kernel, using a direct convolution +/// +/// This is appropriate for small kernels, where there is no time saving to use FFT method. +/// +/// @return psImage* resulting image +/// +psImage *psImageConvolveDirect( + psImage *out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage *in, ///< the psImage to convolve + const psKernel *kernel ///< kernel to colvolve with ); -/** Smooths an image by parts using 1D Gaussian independently in x and y. - * - * Applies a circularly symmetric Gaussian smoothing first in x and then in y - * directions with just a vector. This process is 2N faster than 2D convolutions (in general). - * - * @return bool TRUE if successful, otherwise FALSE - */ +/// Convolve an image with a kernel, using the FFT +/// +/// This is appropriate for larger kernels, where the direct convolution is slow. The input image and kernel +/// are suitably padded to avoid wrap-around effects. +psImage *psImageConvolveFFT( + psImage *out, ///< a psImage to recycle. If NULL, a new psImage is made. + const psImage *in, ///< the psImage to convolve + const psKernel *kernel, ///< kernel to colvolve with + float pad ///< Value to use to pad the input image +); + +/// Smooths an image by parts using 1D Gaussian independently in x and y. +/// +/// Applies a circularly symmetric Gaussian smoothing first in x and then in y +/// directions with just a vector. This process is 2N faster than 2D convolutions (in general). +/// +/// @return bool TRUE if successful, otherwise FALSE +/// bool psImageSmooth( psImage *image, ///< the image to be smoothed Index: trunk/psLib/test/imageops/.cvsignore =================================================================== --- trunk/psLib/test/imageops/.cvsignore (revision 11701) +++ trunk/psLib/test/imageops/.cvsignore (revision 11703) @@ -23,2 +23,4 @@ tap_psImageShift tap_psImageShiftKernel +tap_psImageConvolve +tap_psImageConvolve2 Index: trunk/psLib/test/imageops/Makefile.am =================================================================== --- trunk/psLib/test/imageops/Makefile.am (revision 11701) +++ trunk/psLib/test/imageops/Makefile.am (revision 11703) @@ -21,4 +21,5 @@ tap_psImageStructManip \ tap_psImageConvolve \ + tap_psImageConvolve2 \ tap_psImagePixelExtract Index: trunk/psLib/test/imageops/tap_psImageConvolve2.c =================================================================== --- trunk/psLib/test/imageops/tap_psImageConvolve2.c (revision 11703) +++ trunk/psLib/test/imageops/tap_psImageConvolve2.c (revision 11703) @@ -0,0 +1,149 @@ +#include +#include + +#include "tap.h" +#include "pstap.h" + +#define IMAGE_SIZE 21 +#define KERNEL_SIZE 3 +#define TOL 1.0e-6 + + +// Generate image with single high pixel +static psImage *generateImage(void) +{ + psImage *image = psImageAlloc(IMAGE_SIZE, IMAGE_SIZE, PS_TYPE_F32); + psImageInit(image, 0.0); + image->data.F32[IMAGE_SIZE/2][IMAGE_SIZE/2] = 1.0; + return image; +} + +// Generate Gaussian kernel +static psKernel *generateKernel(void) +{ + psKernel *kernel = psKernelAlloc(- KERNEL_SIZE, KERNEL_SIZE, - KERNEL_SIZE, KERNEL_SIZE); + for (int y = - KERNEL_SIZE; y <= KERNEL_SIZE; y++) { + for (int x = - KERNEL_SIZE; x <= KERNEL_SIZE; x++) { + kernel->kernel[y][x] = exp(- x*x - y*y); + } + } + return kernel; +} + +// Check the convolved image matches the kernel +static bool checkConvolved(const psImage *image, const psKernel *kernel) +{ + bool correct = true; + for (int y = 0; y < IMAGE_SIZE; y++) { + for (int x = 0; x < IMAGE_SIZE; x++) { + if (x < IMAGE_SIZE/2 - KERNEL_SIZE || x > IMAGE_SIZE/2 + KERNEL_SIZE || + y < IMAGE_SIZE/2 - KERNEL_SIZE || y > IMAGE_SIZE/2 + KERNEL_SIZE) { + if (fabs(image->data.F32[y][x]) > TOL) { + diag("%d,%d --> %f", x, y, fabs(image->data.F32[y][x])); + correct = false; + } + } else { + // Position relative to the centre + int kx = x - IMAGE_SIZE/2; + int ky = y - IMAGE_SIZE/2; + if (fabs(image->data.F32[y][x] - kernel->kernel[ky][kx]) > TOL) { + diag("%d,%d --> %f", x, y, fabs(image->data.F32[y][x] - kernel->kernel[ky][kx])); + correct = false; + } + } + } + } + return correct; +} + + +int main(int argc, char *argv[]) +{ + plan_tests(16); + + diag("psImageConvolve tests"); + + { + // Direct convolution + psMemId id = psMemGetId(); + + psImage *image = generateImage(); + psKernel *kernel = generateKernel(); + + psImage *convolved = psImageConvolveDirect(NULL, image, kernel); + ok(convolved, "convolution result"); + skip_start(!convolved, 3, "convolution failed"); + ok(convolved->type.type == PS_TYPE_F32, "output type"); + ok(convolved->numCols == IMAGE_SIZE && convolved->numRows == IMAGE_SIZE, "output size %dx%d", + convolved->numCols, convolved->numRows); + ok(checkConvolved(convolved, kernel), "convolved image correct"); + psFree(convolved); + skip_end(); + + psFree(kernel); + psFree(image); + + ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); + } + + { + // FFT forward, then back --- do I get what I started with? + psMemId id = psMemGetId(); + + psImage *old = generateImage(); + psImage *fftReal = NULL, *fftImag = NULL; + bool result = psImageForwardFFT(&fftReal, &fftImag, old); + ok(result, "forward fft result"); + skip_start(!result, 3, "forward fft failed"); + ok(fftReal->type.type == PS_TYPE_F32 && fftImag->type.type == PS_TYPE_F32, "forward fft types"); + psImage *new = NULL; + result = psImageBackwardFFT(&new, fftReal, fftImag, old->numCols); + ok(result, "backward fft result"); + skip_start(!result, 2, "backward fft failed"); + ok(new->type.type == PS_TYPE_F32, "backward fft type"); + bool correct = true; + for (int y = 0; y < old->numRows; y++) { + for (int x = 0; x < old->numCols; x++) { + if (fabs(new->data.F32[y][x] / IMAGE_SIZE / IMAGE_SIZE - old->data.F32[y][x]) > TOL) { + correct = false; + } + } + } + ok(correct, "new matches old"); + psFree(new); + skip_end(); + skip_end(); + + psFree(fftReal); + psFree(fftImag); + psFree(old); + ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); + + } + + { + // FFT convolution + psMemId id = psMemGetId(); + + psImage *image = generateImage(); + psKernel *kernel = generateKernel(); + + psImage *convolved = psImageConvolveFFT(NULL, image, kernel, 0.0); + ok(convolved, "convolution result"); + skip_start(!convolved, 3, "convolution failed"); + ok(convolved->type.type == PS_TYPE_F32, "output type"); + ok(convolved->numCols == IMAGE_SIZE && convolved->numRows == IMAGE_SIZE, "output size %dx%d", + convolved->numCols, convolved->numRows); + ok(checkConvolved(convolved, kernel), "convolved image correct"); + psFree(convolved); + skip_end(); + + psFree(kernel); + psFree(image); + + ok(!psMemCheckLeaks(id, NULL, NULL, false), "no memory leaks"); + } + + + return exit_status(); +}