diff options
Diffstat (limited to 'cuda/2d')
| -rw-r--r-- | cuda/2d/fft.cu | 70 |
1 files changed, 31 insertions, 39 deletions
diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu index aab755a..3dd1c22 100644 --- a/cuda/2d/fft.cu +++ b/cuda/2d/fft.cu @@ -304,44 +304,42 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, float * pfFilt = new float[_iFFTFourierDetectorCount]; float * pfW = new float[_iFFTFourierDetectorCount]; -#if 1 + // We cache one Fourier transform for repeated FBP's of the same size + static float *pfData = 0; + static int iFilterCacheSize = 0; + + if (!pfData || iFilterCacheSize != _iFFTRealDetectorCount) { + // Compute filter in spatial domain + + delete[] pfData; + pfData = new float[2*_iFFTRealDetectorCount]; + int *ip = new int[int(2+sqrt(_iFFTRealDetectorCount)+1)]; + ip[0] = 0; + float32 *w = new float32[_iFFTRealDetectorCount/2]; + + for (int i = 0; i < _iFFTRealDetectorCount; ++i) { + pfData[2*i+1] = 0.0f; + + if (i & 1) { + int j = i; + if (2*j > _iFFTRealDetectorCount) + j = _iFFTRealDetectorCount - j; + float f = M_PI * j; + pfData[2*i] = -1 / (f*f); + } else { + pfData[2*i] = 0.0f; + } + } - for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount; + pfData[0] = 0.25f; - // filt = 2*( 0:(order/2) )./order; - pfFilt[iDetectorIndex] = 2.0f * fRelIndex; - //pfFilt[iDetectorIndex] = 1.0f; + cdft(2*_iFFTRealDetectorCount, -1, pfData, ip, w); + delete[] ip; + delete[] w; - // w = 2*pi*(0:size(filt,2)-1)/order - pfW[iDetectorIndex] = 3.1415f * 2.0f * fRelIndex; - } -#else - - float *pfData = new float[2*_iFFTRealDetectorCount]; - int *ip = new int[int(2+sqrt(_iFFTRealDetectorCount)+1)]; - ip[0]=0; - float32 *w = new float32[_iFFTRealDetectorCount/2]; - - for (int i = 0; i < _iFFTRealDetectorCount; ++i) { - pfData[2*i+1] = 0.0f; - - if (i & 1) { - int j = i; - if (2*j > _iFFTRealDetectorCount) - j = _iFFTRealDetectorCount - j; - float f = M_PI * j; - pfData[2*i] = -1 / (f*f); - } else { - pfData[2*i] = 0.0f; - } + iFilterCacheSize = _iFFTRealDetectorCount; } - pfData[0] = 0.25f; - - cdft(2*_iFFTRealDetectorCount, -1, pfData, ip, w); - for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount; @@ -350,11 +348,6 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex; } - delete[] pfData; - delete[] ip; - delete[] w; -#endif - switch(_eFilter) { case FILTER_RAMLAK: @@ -906,5 +899,4 @@ void downloadDebugFilterReal(float * _pfHostSinogram, int _iProjectionCount, free(pfHostFilter); } - #endif |