Changeset 19058
- Timestamp:
- Aug 13, 2008, 5:23:13 PM (18 years ago)
- Location:
- trunk/psLib/src/fft
- Files:
-
- 4 edited
-
psFFT.c (modified) (1 diff)
-
psFFT.h (modified) (1 diff)
-
psImageFFT.c (modified) (16 diffs)
-
psVectorFFT.c (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/fft/psFFT.c
r19004 r19058 4 4 5 5 #include <fftw3.h> 6 #include <pthread.h> 6 7 7 8 #include "psAssert.h" 8 9 #include "psLogMsg.h" 9 10 #include "psFFT.h" 11 12 static pthread_mutex_t fftwLock = PTHREAD_MUTEX_INITIALIZER; // Lock for FFTW 13 14 void psFFTLock(void) 15 { 16 pthread_mutex_lock(&fftwLock); 17 } 18 19 void psFFTUnlock(void) 20 { 21 pthread_mutex_unlock(&fftwLock); 22 } 10 23 11 24 bool psFFTThreads(int threads) -
trunk/psLib/src/fft/psFFT.h
r18956 r19058 2 2 #define PS_FFT_H 3 3 4 /// Set number of threads to use with all FFTs 4 // Only fftw_execute is thread-safe: all other calls have to be wrapped in a mutex 5 6 /// Lock FFT mutex 7 void psFFTLock(void); 8 9 /// Unlock FFT mutex 10 void psFFTUnlock(void); 11 12 /// Set number of threads for FFTW to use 5 13 /// 6 /// Pass zero to turn off threading 14 /// Pass zero to turn off threading. Note that the FFTW threads are *independent* of any other threads set up 15 /// throuhg, e.g., psThread. 7 16 bool psFFTThreads(int threads // Number of threads 8 17 ); -
trunk/psLib/src/fft/psImageFFT.c
r17377 r19058 6 6 /// @author Robert DeSonia, MHPCC 7 7 /// 8 /// @version $Revision: 1.2 6$ $Name: not supported by cvs2svn $9 /// @date $Date: 2008-0 4-08 01:25:42$8 /// @version $Revision: 1.27 $ $Name: not supported by cvs2svn $ 9 /// @date $Date: 2008-08-14 03:23:13 $ 10 10 /// 11 11 /// Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 20 20 #include <complex.h> 21 21 #include <fftw3.h> 22 #include <pthread.h> 22 23 23 24 #include "psAssert.h" … … 28 29 #include "psImageStructManip.h" 29 30 #include "psImageConvolve.h" 31 #include "psThread.h" 32 #include "psFFT.h" 30 33 #include "psImageFFT.h" 31 34 35 // Lock FFTW access 36 #define FFTW_LOCK \ 37 if (threaded) { \ 38 psFFTLock(); \ 39 } 40 // Unlock FFTW access 41 #define FFTW_UNLOCK \ 42 if (threaded) { \ 43 psFFTUnlock(); \ 44 } 32 45 33 46 #define FFTW_PLAN_RIGOR FFTW_ESTIMATE // How rigorous the FFTW planning is … … 42 55 PS_ASSERT_PTR_NON_NULL(imag, false); 43 56 57 bool threaded = psThreadPoolSize(); // Are we running threaded? 58 44 59 // Make sure the system-level wisdom information is imported. 60 FFTW_LOCK; 45 61 if (!fftwWisdomImported) { 46 62 fftwf_import_system_wisdom(); 47 63 fftwWisdomImported = true; 48 64 } 65 FFTW_UNLOCK; 49 66 50 67 int numCols = in->numCols; // Number of columns … … 64 81 65 82 // Do the FFT 83 84 FFTW_LOCK; 66 85 fftwf_complex *out = fftwf_malloc((numCols/2 + 1) * numRows * sizeof(fftwf_complex)); // Output data 67 86 fftwf_plan plan = fftwf_plan_dft_r2c_2d(numRows, numCols, input, out, FFTW_PLAN_RIGOR); 87 FFTW_UNLOCK; 88 68 89 fftwf_execute(plan); 90 91 FFTW_LOCK; 69 92 fftwf_destroy_plan(plan); 93 FFTW_UNLOCK; 94 70 95 psFree(input); 71 96 … … 87 112 } 88 113 } 114 115 FFTW_LOCK; 89 116 fftwf_free(out); 117 FFTW_UNLOCK; 90 118 91 119 return true; … … 101 129 PS_ASSERT_PTR_NON_NULL(out, false); 102 130 131 bool threaded = psThreadPoolSize(); // Are we running threaded? 132 103 133 int numCols = real->numCols; // Number of columns 104 134 int numRows = real->numRows; // Number of rows … … 120 150 } 121 151 122 123 152 // Make sure the system-level wisdom information is imported. 153 FFTW_LOCK; 124 154 if (!fftwWisdomImported) { 125 155 fftwf_import_system_wisdom(); 126 156 fftwWisdomImported = true; 127 157 } 158 FFTW_UNLOCK; 128 159 129 160 // Stuff the real and imaginary parts in 130 161 psF32 *target = psAlloc(2 * numCols * numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Target for FFTW 162 FFTW_LOCK; 131 163 fftwf_complex *in = fftwf_malloc(numCols * numRows * sizeof(fftwf_complex)); // Input data 164 FFTW_UNLOCK; 132 165 133 166 for (int y = 0, index = 0; y < numRows; y++) { … … 145 178 146 179 // Do the FFT 180 FFTW_LOCK; 147 181 fftwf_plan plan = fftwf_plan_dft_c2r_2d(numRows, origCols, in, target, FFTW_PLAN_RIGOR); 182 FFTW_UNLOCK; 183 148 184 fftwf_execute(plan); 185 186 FFTW_LOCK; 149 187 fftwf_destroy_plan(plan); 150 188 fftwf_free(in); 189 FFTW_UNLOCK; 151 190 152 191 // Copy the target pixels into the output … … 288 327 } 289 328 329 bool threaded = psThreadPoolSize(); // Are we running threaded? 330 290 331 int numCols = in->numCols, numRows = in->numRows; // Size of image 291 332 int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Kernel sizes … … 304 345 305 346 // Create data array containing the padded image and padded kernel 347 FFTW_LOCK; 306 348 psF32 *data = fftwf_malloc(2 * numPadded * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Data for FFTW 349 FFTW_UNLOCK; 307 350 psF32 *dataPtr = data; // Pointer into FFTW data 308 351 psF32 **imageData = in->data.F32; // Pointer into image data … … 317 360 memset(dataPtr, 0, (paddedRows - numRows) * paddedCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); 318 361 319 #if 1362 #if 0 320 363 { 321 364 // Use this for inspecting the result of copying the image … … 379 422 } 380 423 381 #if 1424 #if 0 382 425 { 383 426 // Use this for inspecting the result of copying the kernel … … 409 452 // Do the forward FFT 410 453 // Note that the FFT images have different size from the input 454 FFTW_LOCK; 411 455 fftwf_complex *fft = fftwf_malloc(2 * (paddedCols/2 + 1) * paddedRows * sizeof(fftwf_complex)); // FFT 456 FFTW_UNLOCK; 412 457 int size[] = { paddedRows, paddedCols }; // Size of transforms 413 458 int fftCols = paddedCols/2 + 1, fftRows = paddedRows; // Size of FFT images 414 459 int fftPixels = fftCols * fftRows; // Number of pixels in FFT image 460 461 FFTW_LOCK; 415 462 fftwf_plan forward = fftwf_plan_many_dft_r2c(2, size, 2, data, NULL, 1, paddedCols * paddedRows, 416 463 fft, NULL, 1, fftPixels, FFTW_PLAN_RIGOR); 464 FFTW_UNLOCK; 465 417 466 fftwf_execute(forward); 467 468 FFTW_LOCK; 418 469 fftwf_destroy_plan(forward); 470 FFTW_UNLOCK; 419 471 420 472 // Multiply the two transforms … … 434 486 435 487 // Do the backward FFT 488 FFTW_LOCK; 436 489 fftwf_plan backward = fftwf_plan_dft_c2r_2d(paddedRows, paddedCols, fft, data, FFTW_PLAN_RIGOR); 490 FFTW_UNLOCK; 491 437 492 fftwf_execute(backward); 493 494 FFTW_LOCK; 438 495 fftwf_destroy_plan(backward); 439 496 fftwf_free(fft); 497 FFTW_UNLOCK; 440 498 441 499 // Copy into the target, without the padding … … 446 504 memcpy(*outData, dataPtr, goodBytes); 447 505 } 506 507 FFTW_LOCK; 448 508 fftwf_free(data); 509 FFTW_UNLOCK; 449 510 450 511 return out; -
trunk/psLib/src/fft/psVectorFFT.c
r11719 r19058 6 6 * @author Robert DeSonia, MHPCC 7 7 * 8 * @version $Revision: 1. 39$ $Name: not supported by cvs2svn $9 * @date $Date: 200 7-02-09 00:45:51$8 * @version $Revision: 1.40 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2008-08-14 03:23:13 $ 10 10 * 11 11 * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii … … 27 27 #include "psLogMsg.h" 28 28 #include "psConstants.h" 29 #include "psThread.h" 30 #include "psFFT.h" 29 31 #include "psVectorFFT.h" 32 33 // Lock FFTW access 34 #define FFTW_LOCK \ 35 if (threaded) { \ 36 psFFTLock(); \ 37 } 38 // Unlock FFTW access 39 #define FFTW_UNLOCK \ 40 if (threaded) { \ 41 psFFTUnlock(); \ 42 } 30 43 31 44 #define FFTW_PLAN_RIGOR FFTW_ESTIMATE // How rigorous the FFTW planning is … … 40 53 PS_ASSERT_PTR_NON_NULL(imag, false); 41 54 55 bool threaded = psThreadPoolSize(); // Are we running threaded? 56 42 57 // Make sure the system-level wisdom information is imported. 58 FFTW_LOCK; 43 59 if (!fftwWisdomImported) { 44 60 fftwf_import_system_wisdom(); 45 61 fftwWisdomImported = true; 46 62 } 63 FFTW_UNLOCK; 47 64 48 65 long num = in->n; // Number of elements 49 66 50 67 // Do the FFT 68 FFTW_LOCK; 51 69 fftwf_complex *out = fftwf_malloc((num/2 + 1) * sizeof(fftwf_complex)); // Output data 52 70 fftwf_plan plan = fftwf_plan_dft_r2c_1d(num, in->data.F32, out, FFTW_PLAN_RIGOR); 71 FFTW_UNLOCK; 72 53 73 fftwf_execute(plan); 74 75 FFTW_LOCK; 54 76 fftwf_destroy_plan(plan); 77 FFTW_UNLOCK; 55 78 56 79 // Pull the real and imaginary parts out … … 69 92 #endif 70 93 } 94 95 FFTW_LOCK; 71 96 fftwf_free(out); 97 FFTW_UNLOCK; 72 98 73 99 return true; … … 83 109 PS_ASSERT_PTR_NON_NULL(out, false); 84 110 111 bool threaded = psThreadPoolSize(); // Are we running threaded? 112 85 113 // Make sure the system-level wisdom information is imported. 114 FFTW_LOCK; 86 115 if (!fftwWisdomImported) { 87 116 fftwf_import_system_wisdom(); 88 117 fftwWisdomImported = true; 89 118 } 119 FFTW_UNLOCK; 90 120 91 121 long num = real->n; // Number of elements 92 122 93 123 // Stuff the real and imaginary parts in 124 FFTW_LOCK; 94 125 fftwf_complex *in = fftwf_malloc(num * sizeof(fftwf_complex)); // Input data 126 FFTW_UNLOCK; 95 127 for (int i = 0; i < num; i++) { 96 128 #if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I) … … 107 139 // Do the FFT 108 140 *out = psVectorRecycle(*out, origNum, PS_TYPE_F32); 141 FFTW_LOCK; 109 142 fftwf_plan plan = fftwf_plan_dft_c2r_1d(origNum, in, (*out)->data.F32, FFTW_PLAN_RIGOR); 143 FFTW_UNLOCK; 144 110 145 fftwf_execute(plan); 146 147 FFTW_LOCK; 111 148 fftwf_destroy_plan(plan); 112 113 149 fftwf_free(in); 150 FFTW_UNLOCK; 114 151 115 152 return true;
Note:
See TracChangeset
for help on using the changeset viewer.
