diff options
Diffstat (limited to 'cuda/2d/fft.cu')
| -rw-r--r-- | cuda/2d/fft.cu | 396 |
1 files changed, 5 insertions, 391 deletions
diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu index bd8cab5..2e94b79 100644 --- a/cuda/2d/fft.cu +++ b/cuda/2d/fft.cu @@ -275,17 +275,6 @@ bool runCudaIFFT(int _iProjectionCount, const cufftComplex* _pDevSourceComplex, return true; } - -// Because the input is real, the Fourier transform is symmetric. -// CUFFT only outputs the first half (ignoring the redundant second half), -// and expects the same as input for the IFFT. -int calcFFTFourierSize(int _iFFTRealSize) -{ - int iFFTFourierSize = _iFFTRealSize / 2 + 1; - - return iFFTFourierSize; -} - void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount) { @@ -300,387 +289,13 @@ void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter, } } -void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, +void genCuFFTFilter(const SFilterConfig &_cfg, int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, - int _iFFTFourierDetectorCount, float _fParameter /* = -1.0f */) + int _iFFTFourierDetectorCount) { - float * pfFilt = new float[_iFFTFourierDetectorCount]; - float * pfW = new float[_iFFTFourierDetectorCount]; - - // 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; - } - } - - pfData[0] = 0.25f; - - cdft(2*_iFFTRealDetectorCount, -1, pfData, ip, w); - delete[] ip; - delete[] w; - - iFilterCacheSize = _iFFTRealDetectorCount; - } - - for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount; - - pfFilt[iDetectorIndex] = 2.0f * pfData[2*iDetectorIndex]; - pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex; - } - - switch(_eFilter) - { - case FILTER_RAMLAK: - { - // do nothing - break; - } - case FILTER_SHEPPLOGAN: - { - // filt(2:end) = filt(2:end) .* (sin(w(2:end)/(2*d))./(w(2:end)/(2*d))) - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (sinf(pfW[iDetectorIndex] / 2.0f / _fD) / (pfW[iDetectorIndex] / 2.0f / _fD)); - } - break; - } - case FILTER_COSINE: - { - // filt(2:end) = filt(2:end) .* cos(w(2:end)/(2*d)) - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * cosf(pfW[iDetectorIndex] / 2.0f / _fD); - } - break; - } - case FILTER_HAMMING: - { - // filt(2:end) = filt(2:end) .* (.54 + .46 * cos(w(2:end)/d)) - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * ( 0.54f + 0.46f * cosf(pfW[iDetectorIndex] / _fD)); - } - break; - } - case FILTER_HANN: - { - // filt(2:end) = filt(2:end) .*(1+cos(w(2:end)./d)) / 2 - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (1.0f + cosf(pfW[iDetectorIndex] / _fD)) / 2.0f; - } - break; - } - case FILTER_TUKEY: - { - float fAlpha = _fParameter; - if(_fParameter < 0.0f) fAlpha = 0.5f; - float fN = (float)_iFFTFourierDetectorCount; - float fHalfN = fN / 2.0f; - float fEnumTerm = fAlpha * fHalfN; - float fDenom = (1.0f - fAlpha) * fHalfN; - float fBlockStart = fHalfN - fEnumTerm; - float fBlockEnd = fHalfN + fEnumTerm; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fAbsSmallN = fabs((float)iDetectorIndex); - float fStoredValue = 0.0f; - - if((fBlockStart <= fAbsSmallN) && (fAbsSmallN <= fBlockEnd)) - { - fStoredValue = 1.0f; - } - else - { - float fEnum = fAbsSmallN - fEnumTerm; - float fCosInput = M_PI * fEnum / fDenom; - fStoredValue = 0.5f * (1.0f + cosf(fCosInput)); - } - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_LANCZOS: - { - float fDenum = (float)(_iFFTFourierDetectorCount - 1); - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fX = 2.0f * fSmallN / fDenum - 1.0f; - float fSinInput = M_PI * fX; - float fStoredValue = 0.0f; - - if(fabsf(fSinInput) > 0.001f) - { - fStoredValue = sin(fSinInput)/fSinInput; - } - else - { - fStoredValue = 1.0f; - } - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_TRIANGULAR: - { - float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1); - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fAbsInput = fSmallN - fNMinusOne / 2.0f; - float fParenInput = fNMinusOne / 2.0f - fabsf(fAbsInput); - float fStoredValue = 2.0f / fNMinusOne * fParenInput; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_GAUSSIAN: - { - float fSigma = _fParameter; - if(_fParameter < 0.0f) fSigma = 0.4f; - float fN = (float)_iFFTFourierDetectorCount; - float fQuotient = (fN - 1.0f) / 2.0f; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fEnum = fSmallN - fQuotient; - float fDenom = fSigma * fQuotient; - float fPower = -0.5f * (fEnum / fDenom) * (fEnum / fDenom); - float fStoredValue = expf(fPower); - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_BARTLETTHANN: - { - const float fA0 = 0.62f; - const float fA1 = 0.48f; - const float fA2 = 0.38f; - float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fAbsInput = fSmallN / fNMinusOne - 0.5f; - float fFirstTerm = fA1 * fabsf(fAbsInput); - float fCosInput = 2.0f * M_PI * fSmallN / fNMinusOne; - float fSecondTerm = fA2 * cosf(fCosInput); - float fStoredValue = fA0 - fFirstTerm - fSecondTerm; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_BLACKMAN: - { - float fAlpha = _fParameter; - if(_fParameter < 0.0f) fAlpha = 0.16f; - float fA0 = (1.0f - fAlpha) / 2.0f; - float fA1 = 0.5f; - float fA2 = fAlpha / 2.0f; - float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1); - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fCosInput1 = 2.0f * M_PI * 0.5f * fSmallN / fNMinusOne; - float fCosInput2 = 4.0f * M_PI * 0.5f * fSmallN / fNMinusOne; - float fStoredValue = fA0 - fA1 * cosf(fCosInput1) + fA2 * cosf(fCosInput2); - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_NUTTALL: - { - const float fA0 = 0.355768f; - const float fA1 = 0.487396f; - const float fA2 = 0.144232f; - const float fA3 = 0.012604f; - float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fBaseCosInput = M_PI * fSmallN / fNMinusOne; - float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput); - float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput); - float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput); - float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_BLACKMANHARRIS: - { - const float fA0 = 0.35875f; - const float fA1 = 0.48829f; - const float fA2 = 0.14128f; - const float fA3 = 0.01168f; - float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fBaseCosInput = M_PI * fSmallN / fNMinusOne; - float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput); - float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput); - float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput); - float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_BLACKMANNUTTALL: - { - const float fA0 = 0.3635819f; - const float fA1 = 0.4891775f; - const float fA2 = 0.1365995f; - const float fA3 = 0.0106411f; - float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fBaseCosInput = M_PI * fSmallN / fNMinusOne; - float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput); - float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput); - float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput); - float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_FLATTOP: - { - const float fA0 = 1.0f; - const float fA1 = 1.93f; - const float fA2 = 1.29f; - const float fA3 = 0.388f; - const float fA4 = 0.032f; - float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f; - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fBaseCosInput = M_PI * fSmallN / fNMinusOne; - float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput); - float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput); - float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput); - float fFourthTerm = fA4 * cosf(8.0f * fBaseCosInput); - float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm + fFourthTerm; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_KAISER: - { - float fAlpha = _fParameter; - if(_fParameter < 0.0f) fAlpha = 3.0f; - float fPiTimesAlpha = M_PI * fAlpha; - float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1); - float fDenom = (float)j0((double)fPiTimesAlpha); - - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fSquareInput = 2.0f * fSmallN / fNMinusOne - 1; - float fSqrtInput = 1.0f - fSquareInput * fSquareInput; - float fBesselInput = fPiTimesAlpha * sqrt(fSqrtInput); - float fEnum = (float)j0((double)fBesselInput); - float fStoredValue = fEnum / fDenom; - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - case FILTER_PARZEN: - { - for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fSmallN = (float)iDetectorIndex; - float fQ = fSmallN / (float)(_iFFTFourierDetectorCount - 1); - float fStoredValue = 0.0f; - - if(fQ <= 0.5f) - { - fStoredValue = 1.0f - 6.0f * fQ * fQ * (1.0f - fQ); - } - else - { - float fCubedValue = 1.0f - fQ; - fStoredValue = 2.0f * fCubedValue * fCubedValue * fCubedValue; - } - - pfFilt[iDetectorIndex] *= fStoredValue; - } - - break; - } - default: - { - ASTRA_ERROR("Cannot serve requested filter"); - } - } - - // filt(w>pi*d) = 0; - float fPiTimesD = M_PI * _fD; - for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) - { - float fWValue = pfW[iDetectorIndex]; - - if(fWValue > fPiTimesD) - { - pfFilt[iDetectorIndex] = 0.0f; - } - } + float * pfFilt = astra::genFilter(_cfg, + _iFFTRealDetectorCount, + _iFFTFourierDetectorCount); for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { @@ -695,7 +310,6 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, } delete[] pfFilt; - delete[] pfW; } |