From d58f6b51493c7931cbc02feb3ffeb0eca6ea3a4e Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 5 Jul 2018 16:13:37 +0200 Subject: Refactor a few filter-related functions out of cuda code --- cuda/2d/fbp.cu | 2 +- cuda/2d/fft.cu | 383 +-------------------------------------------------------- cuda/3d/fdk.cu | 2 +- 3 files changed, 6 insertions(+), 381 deletions(-) (limited to 'cuda') diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu index 48fb7dc..7a8d2e9 100644 --- a/cuda/2d/fbp.cu +++ b/cuda/2d/fbp.cu @@ -130,7 +130,7 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* case astra::FILTER_FLATTOP: case astra::FILTER_PARZEN: { - genFilter(_eFilter, _fD, dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount, _fFilterParameter); + genCuFFTFilter(_eFilter, _fD, dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount, _fFilterParameter); uploadComplexArrayToDevice(dims.iProjAngles, iFreqBinCount, pHostFilter, (cufftComplex*)D_filter); break; diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu index bd8cab5..4454746 100644 --- a/cuda/2d/fft.cu +++ b/cuda/2d/fft.cu @@ -300,387 +300,13 @@ void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter, } } -void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, +void genCuFFTFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, cufftComplex * _pFilter, int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount, float _fParameter /* = -1.0f */) { - 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(_eFilter, _fD, _iProjectionCount, + _iFFTRealDetectorCount, + _iFFTFourierDetectorCount, _fParameter); for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++) { @@ -695,7 +321,6 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount, } delete[] pfFilt; - delete[] pfW; } diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu index 46c07e7..194d2fb 100644 --- a/cuda/3d/fdk.cu +++ b/cuda/3d/fdk.cu @@ -253,7 +253,7 @@ bool FDK_Filter(cudaPitchedPtr D_projData, memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iHalfFFTSize); if (pfFilter == 0){ - astraCUDA::genFilter(astra::FILTER_RAMLAK, 1.0f, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize); + astraCUDA::genCuFFTFilter(astra::FILTER_RAMLAK, 1.0f, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize); } else { for (int i = 0; i < dims.iProjAngles * iHalfFFTSize; i++) { pHostFilter[i].x = pfFilter[i]; -- cgit v1.2.3