Changeset 1983
- Timestamp:
- Oct 6, 2004, 11:31:30 AM (22 years ago)
- Location:
- trunk/psLib/src
- Files:
-
- 6 edited
-
fft/psImageFFT.c (modified) (6 diffs)
-
image/psImageConvolve.c (modified) (7 diffs)
-
image/psImageErrors.dat (modified) (1 diff)
-
image/psImageErrors.h (modified) (2 diffs)
-
image/psImageFFT.c (modified) (6 diffs)
-
imageops/psImageConvolve.c (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/fft/psImageFFT.c
r1864 r1983 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $8 * @date $Date: 2004- 09-23 18:31:49$7 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-10-06 21:31:30 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 21 21 #include "psLogMsg.h" 22 22 #include "psImageExtraction.h" 23 #include "psImageIO.h" 23 24 24 25 #include "psImageErrors.h" … … 41 42 } 42 43 44 if ( ((direction & PS_FFT_FORWARD) != 0) ) { 45 if ((direction & PS_FFT_REVERSE) != 0) { 46 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", 47 PS_ERR_BAD_PARAMETER_VALUE, true, 48 PS_ERRORTEXT_psImageFFT_FORWARD_REVERSE); 49 psFree(out); 50 return NULL; 51 } 52 if ((direction & PS_FFT_REAL_RESULT) != 0) { 53 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", 54 PS_ERR_BAD_PARAMETER_VALUE, true, 55 PS_ERRORTEXT_psImageFFT_REAL_FORWARD_NOTSUPPORTED); 56 psFree(out); 57 return NULL; 58 } 59 } else if ((direction & PS_FFT_REVERSE) == 0) { 60 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", 61 PS_ERR_BAD_PARAMETER_VALUE, true, 62 PS_ERRORTEXT_psImageFFT_NO_DIRECTION_OPTION); 63 psFree(out); 64 return NULL; 65 } 66 43 67 type = in->type.type; 44 45 if ((type != PS_TYPE_F32) && (type != PS_TYPE_C32)) {46 char* typeStr;47 PS_TYPE_NAME(typeStr,type);48 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",49 PS_ERR_BAD_PARAMETER_TYPE, true,50 PS_ERRORTEXT_psImageFFT_IMAGE_TYPE_UNSUPPORTED,51 typeStr);52 psFree(out);53 return NULL;54 }55 56 if (type != PS_TYPE_C32 && direction == PS_FFT_REVERSE) {57 char* typeStr;58 PS_TYPE_NAME(typeStr,type);59 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",60 PS_ERR_BAD_PARAMETER_TYPE, true,61 PS_ERRORTEXT_psImageFFT_REVERSE_NOT_COMPLEX,62 typeStr);63 psFree(out);64 return NULL;65 66 }67 68 if (type != PS_TYPE_F32 && direction == PS_FFT_FORWARD) {69 char* typeStr;70 PS_TYPE_NAME(typeStr,type);71 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",72 PS_ERR_BAD_PARAMETER_TYPE, true,73 PS_ERRORTEXT_psImageFFT_FORWARD_NOT_REAL,74 typeStr);75 psFree(out);76 return NULL;77 }78 68 79 69 /* make sure the system-level wisdom information is imported. */ … … 86 76 numCols = in->numCols; 87 77 88 // n.b. FFTW can perform a in-place transform at the same or faster than out-of-place. 78 // n.b. FFTW can perform a in-place transform at the same rate or faster than out-of-place. 79 int sign = ((direction & PS_FFT_FORWARD) != 0) ? FFTW_FORWARD : FFTW_BACKWARD; 89 80 out = psImageCopy(out, in, PS_TYPE_C32); 90 81 plan = fftwf_plan_dft_2d(numCols, numRows, 91 (fftwf_complex *) out->data.C32[0], 92 (fftwf_complex *) out->data.C32[0], direction, PS_FFTW_PLAN_RIGOR); 93 94 /* check if a plan exists now -- if not, it is a real problem */ 82 out->data.C32[0], 83 out->data.C32[0], 84 sign, 85 PS_FFTW_PLAN_RIGOR); 86 87 /* check if a plan exists now -- if not, it is a real problem at this point */ 95 88 if (plan == NULL) { 96 89 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", … … 105 98 106 99 fftwf_destroy_plan(plan); 100 101 if (direction & PS_FFT_REAL_RESULT) { 102 // n.b., we do this instead of using fftwf_plan_dft_c2r because that 103 // plan requires a half-image, which would require a image reordering 104 // that is not as simple as performing a normal complex transform and 105 // then taking the real part of the result. If performance here 106 // becomes an issue, the use of fftwf_plan_dft_r2c should be considered 107 // as well as fftwf_plan_dft_c2r. 108 psImage* realOut = psImageCopy(NULL,out,PS_TYPE_F32); 109 psFree(out); 110 out = realOut; 111 } 107 112 108 113 return out; … … 424 429 real = creal(inRow[col]); 425 430 imag = cimag(inRow[col]); 426 outRow[col] = real * real + imag * imag/ numElementsSquared;431 outRow[col] = (real * real + imag * imag) / numElementsSquared; 427 432 } 428 433 } -
trunk/psLib/src/image/psImageConvolve.c
r1940 r1983 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 5$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-10-0 1 20:58:32$7 * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-10-06 21:31:30 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 20 20 #include "psLogMsg.h" 21 21 #include "psError.h" 22 #include "psImageIO.h" 22 23 23 24 #include "psImageErrors.h" … … 315 316 int kRows = kernel->image->numRows; 316 317 int kCols = kernel->image->numCols; 317 int x0 = (paddedCols - kCols) / 2; 318 int y0 = (paddedRows - kRows) / 2; 319 if (x0 < 0 || y0 < 0) { 318 if (kRows >= numRows || kCols >= numCols) { 320 319 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve", 321 320 PS_ERR_BAD_PARAMETER_SIZE, true, 322 321 PS_ERRORTEXT_psImageConvolve_KERNEL_TOO_LARGE, 323 k ernel->image->numCols,kernel->image->numRows,322 kCols,kRows, 324 323 numCols, numRows); 325 324 psFree(out); … … 350 349 } 351 350 351 psImageWriteSection(paddedImage,0,0,0,NULL,0,"paddedImage.fits"); 352 352 353 // pad the kernel to the same size of paddedImage 353 354 psImage* paddedKernel = psImageAlloc(paddedCols,paddedRows,PS_TYPE_KERNEL); 354 355 memset(paddedKernel->data.U8[0],0,sizeof(psKernelType)*numCols*numRows); // zero-out image 355 for (int row=0;row<kRows;row++) { 356 psKernelType* padRow = paddedKernel->data.PS_TYPE_KERNEL_DATA[row+y0]; 357 psKernelType* kernelRow = kernel->image->data.PS_TYPE_KERNEL_DATA[row]; 358 for (int col=0;col<kCols;col++) { 359 padRow[col+x0] = kernelRow[col]; 356 int yMax = kernel->yMax; 357 int xMax = kernel->xMax; 358 for (int row = kernel->yMin; row <= yMax;row++) { 359 int padRow = row; 360 if (padRow < 0) { 361 padRow += paddedRows; 360 362 } 361 } 363 psKernelType* padData = paddedKernel->data.PS_TYPE_KERNEL_DATA[padRow]; 364 psKernelType* kernelRow = kernel->kernel[row]; 365 for (int col = kernel->xMin; col <= xMax; col++) { 366 if (col < 0) { 367 padData[col+paddedCols] = kernelRow[col]; 368 } else { 369 padData[col] = kernelRow[col]; 370 } 371 } 372 } 373 374 psImageWriteSection(paddedKernel,0,0,0,NULL,0,"paddedKernel.fits"); 362 375 363 376 psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD); … … 380 393 381 394 // convolution in fourier domain is just a pixel-wise multiplication 382 ps Image* outFourier = psImageAlloc(paddedCols,paddedRows,inFourier->type.type);383 psBinaryOp(outFourier,inFourier,"*",kernelFourier); 384 385 psImage* complexOut = psImageFFT(NULL, outFourier,PS_FFT_REVERSE);395 psBinaryOp(inFourier,inFourier,"*",kernelFourier); 396 397 psImage* complexOut = psImageFFT(NULL, inFourier, 398 PS_FFT_REVERSE); 386 399 if (complexOut == NULL) { 387 400 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve", … … 392 405 } 393 406 407 { 408 psImage* tmp = psImageReal(NULL,complexOut); 409 psImageWriteSection(tmp,0,0,0,NULL,0,"complexOut.fits"); 410 psFree(tmp); 411 } 412 394 413 // subset out the padded area now. 395 414 psImage* complexOutSansPad = psImageSubset(complexOut, … … 397 416 FOURIER_PADDING+numCols,FOURIER_PADDING+numRows); 398 417 399 // since our image was real, the result is to be real as well 400 out = psImageReal(out,complexOutSansPad); 418 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 419 float factor = 1.0f/numCols/numRows; 420 for (int row = 0; row < numRows; row++) { 421 psF32* outRow = out->data.F32[row]; 422 psC32* resultRow = complexOutSansPad->data.C32[row]; 423 for (int col = 0; col < numCols; col++) { 424 outRow[col] = crealf(resultRow[col])*factor; 425 } 426 } 401 427 402 428 psFree(complexOut); // frees complexOutSansPad, as it is a child. -
trunk/psLib/src/image/psImageErrors.dat
r1951 r1983 41 41 psImageFFT_NONREAL_NOTSUPPORTED Input psImage type, %s, is required to be either psF32 or psF64. 42 42 psImageFFT_NONCOMPLEX_NOTSUPPORTED Input psImage type, %s, is required to be either psC32 or psC64. 43 psImageFFT_REAL_FORWARD_NOTSUPPORTED The PS_FFT_FORWARD and PS_FFT_REAL_RESULT combinition is not supported. 44 psImageFFT_FORWARD_REVERSE Can not specify both PS_FFT_FORWARD and PS_FFT_REVERSE options. 45 psImageFFT_NO_DIRECTION_OPTION Must specify either PS_FFT_FORWARD or PS_FFT_REVERSE option. 43 46 # 44 47 psImageIO_FILENAME_NULL Specified filename can not be NULL. -
trunk/psLib/src/image/psImageErrors.h
r1958 r1983 7 7 * @author Robert DeSonia, MHPCC 8 8 * 9 * @version $Revision: 1. 8$ $Name: not supported by cvs2svn $10 * @date $Date: 2004-10-0 5 21:11:05$9 * @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ 10 * @date $Date: 2004-10-06 21:31:30 $ 11 11 * 12 12 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 62 62 #define PS_ERRORTEXT_psImageFFT_NONREAL_NOTSUPPORTED "Input psImage type, %s, is required to be either psF32 or psF64." 63 63 #define PS_ERRORTEXT_psImageFFT_NONCOMPLEX_NOTSUPPORTED "Input psImage type, %s, is required to be either psC32 or psC64." 64 #define PS_ERRORTEXT_psImageFFT_REAL_FORWARD_NOTSUPPORTED "The PS_FFT_FORWARD and PS_FFT_REAL_RESULT combinition is not supported." 65 #define PS_ERRORTEXT_psImageFFT_FORWARD_REVERSE "Can not specify both PS_FFT_FORWARD and PS_FFT_REVERSE options." 66 #define PS_ERRORTEXT_psImageFFT_NO_DIRECTION_OPTION "Must specify either PS_FFT_FORWARD or PS_FFT_REVERSE option." 64 67 #define PS_ERRORTEXT_psImageIO_FILENAME_NULL "Specified filename can not be NULL." 65 68 #define PS_ERRORTEXT_psImageIO_FILENAME_INVALID "Could not open file,'%s'.\nCFITSIO Error: %s" -
trunk/psLib/src/image/psImageFFT.c
r1864 r1983 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $8 * @date $Date: 2004- 09-23 18:31:49$7 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-10-06 21:31:30 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 21 21 #include "psLogMsg.h" 22 22 #include "psImageExtraction.h" 23 #include "psImageIO.h" 23 24 24 25 #include "psImageErrors.h" … … 41 42 } 42 43 44 if ( ((direction & PS_FFT_FORWARD) != 0) ) { 45 if ((direction & PS_FFT_REVERSE) != 0) { 46 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", 47 PS_ERR_BAD_PARAMETER_VALUE, true, 48 PS_ERRORTEXT_psImageFFT_FORWARD_REVERSE); 49 psFree(out); 50 return NULL; 51 } 52 if ((direction & PS_FFT_REAL_RESULT) != 0) { 53 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", 54 PS_ERR_BAD_PARAMETER_VALUE, true, 55 PS_ERRORTEXT_psImageFFT_REAL_FORWARD_NOTSUPPORTED); 56 psFree(out); 57 return NULL; 58 } 59 } else if ((direction & PS_FFT_REVERSE) == 0) { 60 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", 61 PS_ERR_BAD_PARAMETER_VALUE, true, 62 PS_ERRORTEXT_psImageFFT_NO_DIRECTION_OPTION); 63 psFree(out); 64 return NULL; 65 } 66 43 67 type = in->type.type; 44 45 if ((type != PS_TYPE_F32) && (type != PS_TYPE_C32)) {46 char* typeStr;47 PS_TYPE_NAME(typeStr,type);48 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",49 PS_ERR_BAD_PARAMETER_TYPE, true,50 PS_ERRORTEXT_psImageFFT_IMAGE_TYPE_UNSUPPORTED,51 typeStr);52 psFree(out);53 return NULL;54 }55 56 if (type != PS_TYPE_C32 && direction == PS_FFT_REVERSE) {57 char* typeStr;58 PS_TYPE_NAME(typeStr,type);59 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",60 PS_ERR_BAD_PARAMETER_TYPE, true,61 PS_ERRORTEXT_psImageFFT_REVERSE_NOT_COMPLEX,62 typeStr);63 psFree(out);64 return NULL;65 66 }67 68 if (type != PS_TYPE_F32 && direction == PS_FFT_FORWARD) {69 char* typeStr;70 PS_TYPE_NAME(typeStr,type);71 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",72 PS_ERR_BAD_PARAMETER_TYPE, true,73 PS_ERRORTEXT_psImageFFT_FORWARD_NOT_REAL,74 typeStr);75 psFree(out);76 return NULL;77 }78 68 79 69 /* make sure the system-level wisdom information is imported. */ … … 86 76 numCols = in->numCols; 87 77 88 // n.b. FFTW can perform a in-place transform at the same or faster than out-of-place. 78 // n.b. FFTW can perform a in-place transform at the same rate or faster than out-of-place. 79 int sign = ((direction & PS_FFT_FORWARD) != 0) ? FFTW_FORWARD : FFTW_BACKWARD; 89 80 out = psImageCopy(out, in, PS_TYPE_C32); 90 81 plan = fftwf_plan_dft_2d(numCols, numRows, 91 (fftwf_complex *) out->data.C32[0], 92 (fftwf_complex *) out->data.C32[0], direction, PS_FFTW_PLAN_RIGOR); 93 94 /* check if a plan exists now -- if not, it is a real problem */ 82 out->data.C32[0], 83 out->data.C32[0], 84 sign, 85 PS_FFTW_PLAN_RIGOR); 86 87 /* check if a plan exists now -- if not, it is a real problem at this point */ 95 88 if (plan == NULL) { 96 89 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT", … … 105 98 106 99 fftwf_destroy_plan(plan); 100 101 if (direction & PS_FFT_REAL_RESULT) { 102 // n.b., we do this instead of using fftwf_plan_dft_c2r because that 103 // plan requires a half-image, which would require a image reordering 104 // that is not as simple as performing a normal complex transform and 105 // then taking the real part of the result. If performance here 106 // becomes an issue, the use of fftwf_plan_dft_r2c should be considered 107 // as well as fftwf_plan_dft_c2r. 108 psImage* realOut = psImageCopy(NULL,out,PS_TYPE_F32); 109 psFree(out); 110 out = realOut; 111 } 107 112 108 113 return out; … … 424 429 real = creal(inRow[col]); 425 430 imag = cimag(inRow[col]); 426 outRow[col] = real * real + imag * imag/ numElementsSquared;431 outRow[col] = (real * real + imag * imag) / numElementsSquared; 427 432 } 428 433 } -
trunk/psLib/src/imageops/psImageConvolve.c
r1940 r1983 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 5$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-10-0 1 20:58:32$7 * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-10-06 21:31:30 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 20 20 #include "psLogMsg.h" 21 21 #include "psError.h" 22 #include "psImageIO.h" 22 23 23 24 #include "psImageErrors.h" … … 315 316 int kRows = kernel->image->numRows; 316 317 int kCols = kernel->image->numCols; 317 int x0 = (paddedCols - kCols) / 2; 318 int y0 = (paddedRows - kRows) / 2; 319 if (x0 < 0 || y0 < 0) { 318 if (kRows >= numRows || kCols >= numCols) { 320 319 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve", 321 320 PS_ERR_BAD_PARAMETER_SIZE, true, 322 321 PS_ERRORTEXT_psImageConvolve_KERNEL_TOO_LARGE, 323 k ernel->image->numCols,kernel->image->numRows,322 kCols,kRows, 324 323 numCols, numRows); 325 324 psFree(out); … … 350 349 } 351 350 351 psImageWriteSection(paddedImage,0,0,0,NULL,0,"paddedImage.fits"); 352 352 353 // pad the kernel to the same size of paddedImage 353 354 psImage* paddedKernel = psImageAlloc(paddedCols,paddedRows,PS_TYPE_KERNEL); 354 355 memset(paddedKernel->data.U8[0],0,sizeof(psKernelType)*numCols*numRows); // zero-out image 355 for (int row=0;row<kRows;row++) { 356 psKernelType* padRow = paddedKernel->data.PS_TYPE_KERNEL_DATA[row+y0]; 357 psKernelType* kernelRow = kernel->image->data.PS_TYPE_KERNEL_DATA[row]; 358 for (int col=0;col<kCols;col++) { 359 padRow[col+x0] = kernelRow[col]; 356 int yMax = kernel->yMax; 357 int xMax = kernel->xMax; 358 for (int row = kernel->yMin; row <= yMax;row++) { 359 int padRow = row; 360 if (padRow < 0) { 361 padRow += paddedRows; 360 362 } 361 } 363 psKernelType* padData = paddedKernel->data.PS_TYPE_KERNEL_DATA[padRow]; 364 psKernelType* kernelRow = kernel->kernel[row]; 365 for (int col = kernel->xMin; col <= xMax; col++) { 366 if (col < 0) { 367 padData[col+paddedCols] = kernelRow[col]; 368 } else { 369 padData[col] = kernelRow[col]; 370 } 371 } 372 } 373 374 psImageWriteSection(paddedKernel,0,0,0,NULL,0,"paddedKernel.fits"); 362 375 363 376 psImage* kernelFourier = psImageFFT(NULL, paddedKernel, PS_FFT_FORWARD); … … 380 393 381 394 // convolution in fourier domain is just a pixel-wise multiplication 382 ps Image* outFourier = psImageAlloc(paddedCols,paddedRows,inFourier->type.type);383 psBinaryOp(outFourier,inFourier,"*",kernelFourier); 384 385 psImage* complexOut = psImageFFT(NULL, outFourier,PS_FFT_REVERSE);395 psBinaryOp(inFourier,inFourier,"*",kernelFourier); 396 397 psImage* complexOut = psImageFFT(NULL, inFourier, 398 PS_FFT_REVERSE); 386 399 if (complexOut == NULL) { 387 400 psErrorMsg(PS_ERRORNAME_DOMAIN "psImageConvolve", … … 392 405 } 393 406 407 { 408 psImage* tmp = psImageReal(NULL,complexOut); 409 psImageWriteSection(tmp,0,0,0,NULL,0,"complexOut.fits"); 410 psFree(tmp); 411 } 412 394 413 // subset out the padded area now. 395 414 psImage* complexOutSansPad = psImageSubset(complexOut, … … 397 416 FOURIER_PADDING+numCols,FOURIER_PADDING+numRows); 398 417 399 // since our image was real, the result is to be real as well 400 out = psImageReal(out,complexOutSansPad); 418 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 419 float factor = 1.0f/numCols/numRows; 420 for (int row = 0; row < numRows; row++) { 421 psF32* outRow = out->data.F32[row]; 422 psC32* resultRow = complexOutSansPad->data.C32[row]; 423 for (int col = 0; col < numCols; col++) { 424 outRow[col] = crealf(resultRow[col])*factor; 425 } 426 } 401 427 402 428 psFree(complexOut); // frees complexOutSansPad, as it is a child.
Note:
See TracChangeset
for help on using the changeset viewer.
