summaryrefslogtreecommitdiffstats
path: root/cuda/2d/fft.cu
diff options
context:
space:
mode:
Diffstat (limited to 'cuda/2d/fft.cu')
-rw-r--r--cuda/2d/fft.cu383
1 files changed, 4 insertions, 379 deletions
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;
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-01 06:25:21 +0000