Index: /trunk/psLib/src/fft/psImageFFT.c =================================================================== --- /trunk/psLib/src/fft/psImageFFT.c (revision 1982) +++ /trunk/psLib/src/fft/psImageFFT.c (revision 1983) @@ -5,6 +5,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-23 18:31:49 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-06 21:31:30 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -21,4 +21,5 @@ #include "psLogMsg.h" #include "psImageExtraction.h" +#include "psImageIO.h" #include "psImageErrors.h" @@ -41,39 +42,28 @@ } + if ( ((direction & PS_FFT_FORWARD) != 0) ) { + if ((direction & PS_FFT_REVERSE) != 0) { + psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", + PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_FORWARD_REVERSE); + psFree(out); + return NULL; + } + if ((direction & PS_FFT_REAL_RESULT) != 0) { + psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", + PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_REAL_FORWARD_NOTSUPPORTED); + psFree(out); + return NULL; + } + } else if ((direction & PS_FFT_REVERSE) == 0) { + psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", + PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_NO_DIRECTION_OPTION); + psFree(out); + return NULL; + } + type = in->type.type; - - if ((type != PS_TYPE_F32) && (type != PS_TYPE_C32)) { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", - PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psImageFFT_IMAGE_TYPE_UNSUPPORTED, - typeStr); - psFree(out); - return NULL; - } - - if (type != PS_TYPE_C32 && direction == PS_FFT_REVERSE) { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", - PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psImageFFT_REVERSE_NOT_COMPLEX, - typeStr); - psFree(out); - return NULL; - - } - - if (type != PS_TYPE_F32 && direction == PS_FFT_FORWARD) { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", - PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psImageFFT_FORWARD_NOT_REAL, - typeStr); - psFree(out); - return NULL; - } /* make sure the system-level wisdom information is imported. */ @@ -86,11 +76,14 @@ numCols = in->numCols; - // n.b. FFTW can perform a in-place transform at the same or faster than out-of-place. + // n.b. FFTW can perform a in-place transform at the same rate or faster than out-of-place. + int sign = ((direction & PS_FFT_FORWARD) != 0) ? FFTW_FORWARD : FFTW_BACKWARD; out = psImageCopy(out, in, PS_TYPE_C32); plan = fftwf_plan_dft_2d(numCols, numRows, - (fftwf_complex *) out->data.C32[0], - (fftwf_complex *) out->data.C32[0], direction, PS_FFTW_PLAN_RIGOR); - - /* check if a plan exists now -- if not, it is a real problem */ + out->data.C32[0], + out->data.C32[0], + sign, + PS_FFTW_PLAN_RIGOR); + + /* check if a plan exists now -- if not, it is a real problem at this point */ if (plan == NULL) { psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", @@ -105,4 +98,16 @@ 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. + psImage* realOut = psImageCopy(NULL,out,PS_TYPE_F32); + psFree(out); + out = realOut; + } return out; @@ -424,5 +429,5 @@ real = creal(inRow[col]); imag = cimag(inRow[col]); - outRow[col] = real * real + imag * imag / numElementsSquared; + outRow[col] = (real * real + imag * imag) / numElementsSquared; } } Index: /trunk/psLib/src/image/psImageConvolve.c =================================================================== --- /trunk/psLib/src/image/psImageConvolve.c (revision 1982) +++ /trunk/psLib/src/image/psImageConvolve.c (revision 1983) @@ -5,6 +5,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-01 20:58:32 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-06 21:31:30 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -20,4 +20,5 @@ #include "psLogMsg.h" #include "psError.h" +#include "psImageIO.h" #include "psImageErrors.h" @@ -315,11 +316,9 @@ int kRows = kernel->image->numRows; int kCols = kernel->image->numCols; - int x0 = (paddedCols - kCols) / 2; - int y0 = (paddedRows - kRows) / 2; - if (x0 < 0 || y0 < 0) { + if (kRows >= numRows || kCols >= numCols) { psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve", PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_psImageConvolve_KERNEL_TOO_LARGE, - kernel->image->numCols,kernel->image->numRows, + kCols,kRows, numCols, numRows); psFree(out); @@ -350,14 +349,28 @@ } + psImageWriteSection(paddedImage,0,0,0,NULL,0,"paddedImage.fits"); + // pad the kernel to the same size of paddedImage psImage* paddedKernel = psImageAlloc(paddedCols,paddedRows,PS_TYPE_KERNEL); memset(paddedKernel->data.U8[0],0,sizeof(psKernelType)*numCols*numRows); // zero-out image - for (int row=0;rowdata.PS_TYPE_KERNEL_DATA[row+y0]; - psKernelType* kernelRow = kernel->image->data.PS_TYPE_KERNEL_DATA[row]; - for (int col=0;colyMax; + int xMax = kernel->xMax; + for (int row = kernel->yMin; row <= yMax;row++) { + int padRow = row; + if (padRow < 0) { + padRow += paddedRows; } - } + psKernelType* padData = paddedKernel->data.PS_TYPE_KERNEL_DATA[padRow]; + psKernelType* kernelRow = kernel->kernel[row]; + for (int col = kernel->xMin; col <= xMax; col++) { + if (col < 0) { + padData[col+paddedCols] = kernelRow[col]; + } else { + padData[col] = kernelRow[col]; + } + } + } + + psImageWriteSection(paddedKernel,0,0,0,NULL,0,"paddedKernel.fits"); psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD); @@ -380,8 +393,8 @@ // convolution in fourier domain is just a pixel-wise multiplication - psImage* outFourier = psImageAlloc(paddedCols,paddedRows,inFourier->type.type); - psBinaryOp(outFourier,inFourier,"*",kernelFourier); - - psImage* complexOut = psImageFFT(NULL, outFourier, PS_FFT_REVERSE); + psBinaryOp(inFourier,inFourier,"*",kernelFourier); + + psImage* complexOut = psImageFFT(NULL, inFourier, + PS_FFT_REVERSE); if (complexOut == NULL) { psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve", @@ -392,4 +405,10 @@ } + { + psImage* tmp = psImageReal(NULL,complexOut); + psImageWriteSection(tmp,0,0,0,NULL,0,"complexOut.fits"); + psFree(tmp); + } + // subset out the padded area now. psImage* complexOutSansPad = psImageSubset(complexOut, @@ -397,6 +416,13 @@ FOURIER_PADDING+numCols,FOURIER_PADDING+numRows); - // since our image was real, the result is to be real as well - out = psImageReal(out,complexOutSansPad); + out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); + float factor = 1.0f/numCols/numRows; + for (int row = 0; row < numRows; row++) { + psF32* outRow = out->data.F32[row]; + psC32* resultRow = complexOutSansPad->data.C32[row]; + for (int col = 0; col < numCols; col++) { + outRow[col] = crealf(resultRow[col])*factor; + } + } psFree(complexOut); // frees complexOutSansPad, as it is a child. Index: /trunk/psLib/src/image/psImageErrors.dat =================================================================== --- /trunk/psLib/src/image/psImageErrors.dat (revision 1982) +++ /trunk/psLib/src/image/psImageErrors.dat (revision 1983) @@ -41,4 +41,7 @@ psImageFFT_NONREAL_NOTSUPPORTED Input psImage type, %s, is required to be either psF32 or psF64. psImageFFT_NONCOMPLEX_NOTSUPPORTED Input psImage type, %s, is required to be either psC32 or psC64. +psImageFFT_REAL_FORWARD_NOTSUPPORTED The PS_FFT_FORWARD and PS_FFT_REAL_RESULT combinition is not supported. +psImageFFT_FORWARD_REVERSE Can not specify both PS_FFT_FORWARD and PS_FFT_REVERSE options. +psImageFFT_NO_DIRECTION_OPTION Must specify either PS_FFT_FORWARD or PS_FFT_REVERSE option. # psImageIO_FILENAME_NULL Specified filename can not be NULL. Index: /trunk/psLib/src/image/psImageErrors.h =================================================================== --- /trunk/psLib/src/image/psImageErrors.h (revision 1982) +++ /trunk/psLib/src/image/psImageErrors.h (revision 1983) @@ -7,6 +7,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-05 21:11:05 $ + * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-06 21:31:30 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -62,4 +62,7 @@ #define PS_ERRORTEXT_psImageFFT_NONREAL_NOTSUPPORTED "Input psImage type, %s, is required to be either psF32 or psF64." #define PS_ERRORTEXT_psImageFFT_NONCOMPLEX_NOTSUPPORTED "Input psImage type, %s, is required to be either psC32 or psC64." +#define PS_ERRORTEXT_psImageFFT_REAL_FORWARD_NOTSUPPORTED "The PS_FFT_FORWARD and PS_FFT_REAL_RESULT combinition is not supported." +#define PS_ERRORTEXT_psImageFFT_FORWARD_REVERSE "Can not specify both PS_FFT_FORWARD and PS_FFT_REVERSE options." +#define PS_ERRORTEXT_psImageFFT_NO_DIRECTION_OPTION "Must specify either PS_FFT_FORWARD or PS_FFT_REVERSE option." #define PS_ERRORTEXT_psImageIO_FILENAME_NULL "Specified filename can not be NULL." #define PS_ERRORTEXT_psImageIO_FILENAME_INVALID "Could not open file,'%s'.\nCFITSIO Error: %s" Index: /trunk/psLib/src/image/psImageFFT.c =================================================================== --- /trunk/psLib/src/image/psImageFFT.c (revision 1982) +++ /trunk/psLib/src/image/psImageFFT.c (revision 1983) @@ -5,6 +5,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-23 18:31:49 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-06 21:31:30 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -21,4 +21,5 @@ #include "psLogMsg.h" #include "psImageExtraction.h" +#include "psImageIO.h" #include "psImageErrors.h" @@ -41,39 +42,28 @@ } + if ( ((direction & PS_FFT_FORWARD) != 0) ) { + if ((direction & PS_FFT_REVERSE) != 0) { + psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", + PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_FORWARD_REVERSE); + psFree(out); + return NULL; + } + if ((direction & PS_FFT_REAL_RESULT) != 0) { + psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", + PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_REAL_FORWARD_NOTSUPPORTED); + psFree(out); + return NULL; + } + } else if ((direction & PS_FFT_REVERSE) == 0) { + psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", + PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psImageFFT_NO_DIRECTION_OPTION); + psFree(out); + return NULL; + } + type = in->type.type; - - if ((type != PS_TYPE_F32) && (type != PS_TYPE_C32)) { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", - PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psImageFFT_IMAGE_TYPE_UNSUPPORTED, - typeStr); - psFree(out); - return NULL; - } - - if (type != PS_TYPE_C32 && direction == PS_FFT_REVERSE) { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", - PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psImageFFT_REVERSE_NOT_COMPLEX, - typeStr); - psFree(out); - return NULL; - - } - - if (type != PS_TYPE_F32 && direction == PS_FFT_FORWARD) { - char* typeStr; - PS_TYPE_NAME(typeStr,type); - psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", - PS_ERR_BAD_PARAMETER_TYPE, true, - PS_ERRORTEXT_psImageFFT_FORWARD_NOT_REAL, - typeStr); - psFree(out); - return NULL; - } /* make sure the system-level wisdom information is imported. */ @@ -86,11 +76,14 @@ numCols = in->numCols; - // n.b. FFTW can perform a in-place transform at the same or faster than out-of-place. + // n.b. FFTW can perform a in-place transform at the same rate or faster than out-of-place. + int sign = ((direction & PS_FFT_FORWARD) != 0) ? FFTW_FORWARD : FFTW_BACKWARD; out = psImageCopy(out, in, PS_TYPE_C32); plan = fftwf_plan_dft_2d(numCols, numRows, - (fftwf_complex *) out->data.C32[0], - (fftwf_complex *) out->data.C32[0], direction, PS_FFTW_PLAN_RIGOR); - - /* check if a plan exists now -- if not, it is a real problem */ + out->data.C32[0], + out->data.C32[0], + sign, + PS_FFTW_PLAN_RIGOR); + + /* check if a plan exists now -- if not, it is a real problem at this point */ if (plan == NULL) { psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", @@ -105,4 +98,16 @@ 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. + psImage* realOut = psImageCopy(NULL,out,PS_TYPE_F32); + psFree(out); + out = realOut; + } return out; @@ -424,5 +429,5 @@ real = creal(inRow[col]); imag = cimag(inRow[col]); - outRow[col] = real * real + imag * imag / numElementsSquared; + outRow[col] = (real * real + imag * imag) / numElementsSquared; } } Index: /trunk/psLib/src/imageops/psImageConvolve.c =================================================================== --- /trunk/psLib/src/imageops/psImageConvolve.c (revision 1982) +++ /trunk/psLib/src/imageops/psImageConvolve.c (revision 1983) @@ -5,6 +5,6 @@ * @author Robert DeSonia, MHPCC * - * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-01 20:58:32 $ + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-06 21:31:30 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -20,4 +20,5 @@ #include "psLogMsg.h" #include "psError.h" +#include "psImageIO.h" #include "psImageErrors.h" @@ -315,11 +316,9 @@ int kRows = kernel->image->numRows; int kCols = kernel->image->numCols; - int x0 = (paddedCols - kCols) / 2; - int y0 = (paddedRows - kRows) / 2; - if (x0 < 0 || y0 < 0) { + if (kRows >= numRows || kCols >= numCols) { psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve", PS_ERR_BAD_PARAMETER_SIZE, true, PS_ERRORTEXT_psImageConvolve_KERNEL_TOO_LARGE, - kernel->image->numCols,kernel->image->numRows, + kCols,kRows, numCols, numRows); psFree(out); @@ -350,14 +349,28 @@ } + psImageWriteSection(paddedImage,0,0,0,NULL,0,"paddedImage.fits"); + // pad the kernel to the same size of paddedImage psImage* paddedKernel = psImageAlloc(paddedCols,paddedRows,PS_TYPE_KERNEL); memset(paddedKernel->data.U8[0],0,sizeof(psKernelType)*numCols*numRows); // zero-out image - for (int row=0;rowdata.PS_TYPE_KERNEL_DATA[row+y0]; - psKernelType* kernelRow = kernel->image->data.PS_TYPE_KERNEL_DATA[row]; - for (int col=0;colyMax; + int xMax = kernel->xMax; + for (int row = kernel->yMin; row <= yMax;row++) { + int padRow = row; + if (padRow < 0) { + padRow += paddedRows; } - } + psKernelType* padData = paddedKernel->data.PS_TYPE_KERNEL_DATA[padRow]; + psKernelType* kernelRow = kernel->kernel[row]; + for (int col = kernel->xMin; col <= xMax; col++) { + if (col < 0) { + padData[col+paddedCols] = kernelRow[col]; + } else { + padData[col] = kernelRow[col]; + } + } + } + + psImageWriteSection(paddedKernel,0,0,0,NULL,0,"paddedKernel.fits"); psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD); @@ -380,8 +393,8 @@ // convolution in fourier domain is just a pixel-wise multiplication - psImage* outFourier = psImageAlloc(paddedCols,paddedRows,inFourier->type.type); - psBinaryOp(outFourier,inFourier,"*",kernelFourier); - - psImage* complexOut = psImageFFT(NULL, outFourier, PS_FFT_REVERSE); + psBinaryOp(inFourier,inFourier,"*",kernelFourier); + + psImage* complexOut = psImageFFT(NULL, inFourier, + PS_FFT_REVERSE); if (complexOut == NULL) { psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve", @@ -392,4 +405,10 @@ } + { + psImage* tmp = psImageReal(NULL,complexOut); + psImageWriteSection(tmp,0,0,0,NULL,0,"complexOut.fits"); + psFree(tmp); + } + // subset out the padded area now. psImage* complexOutSansPad = psImageSubset(complexOut, @@ -397,6 +416,13 @@ FOURIER_PADDING+numCols,FOURIER_PADDING+numRows); - // since our image was real, the result is to be real as well - out = psImageReal(out,complexOutSansPad); + out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); + float factor = 1.0f/numCols/numRows; + for (int row = 0; row < numRows; row++) { + psF32* outRow = out->data.F32[row]; + psC32* resultRow = complexOutSansPad->data.C32[row]; + for (int col = 0; col < numCols; col++) { + outRow[col] = crealf(resultRow[col])*factor; + } + } psFree(complexOut); // frees complexOutSansPad, as it is a child.