summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--astra_vc14.vcxproj3
-rw-r--r--astra_vc14.vcxproj.filters9
-rw-r--r--build/linux/Makefile.in1
-rw-r--r--build/msvc/gen.py3
-rw-r--r--cuda/2d/fbp.cu55
-rw-r--r--cuda/2d/fft.cu396
-rw-r--r--cuda/3d/fdk.cu6
-rw-r--r--include/astra/CudaFilteredBackProjectionAlgorithm.h9
-rw-r--r--include/astra/FilteredBackProjectionAlgorithm.h5
-rw-r--r--include/astra/Filters.h (renamed from include/astra/cuda/2d/fbp_filters.h)43
-rw-r--r--include/astra/cuda/2d/astra.h1
-rw-r--r--include/astra/cuda/2d/fbp.h6
-rw-r--r--include/astra/cuda/2d/fft.h10
-rw-r--r--samples/matlab/s023_FBP_filters.m96
-rw-r--r--samples/python/s023_FBP_filters.py116
-rw-r--r--src/CudaFDKAlgorithm3D.cpp4
-rw-r--r--src/CudaFilteredBackProjectionAlgorithm.cpp199
-rw-r--r--src/FilteredBackProjectionAlgorithm.cpp137
-rw-r--r--src/Filters.cpp608
19 files changed, 1050 insertions, 657 deletions
diff --git a/astra_vc14.vcxproj b/astra_vc14.vcxproj
index 64c08ee..893305f 100644
--- a/astra_vc14.vcxproj
+++ b/astra_vc14.vcxproj
@@ -509,6 +509,7 @@
<ClCompile Include="src\FanFlatProjectionGeometry2D.cpp" />
<ClCompile Include="src\FanFlatVecProjectionGeometry2D.cpp" />
<ClCompile Include="src\FilteredBackProjectionAlgorithm.cpp" />
+ <ClCompile Include="src\Filters.cpp" />
<ClCompile Include="src\Float32Data.cpp" />
<ClCompile Include="src\Float32Data2D.cpp" />
<ClCompile Include="src\Float32Data3D.cpp" />
@@ -596,6 +597,7 @@
<ClInclude Include="include\astra\FanFlatProjectionGeometry2D.h" />
<ClInclude Include="include\astra\FanFlatVecProjectionGeometry2D.h" />
<ClInclude Include="include\astra\FilteredBackProjectionAlgorithm.h" />
+ <ClInclude Include="include\astra\Filters.h" />
<ClInclude Include="include\astra\Float32Data.h" />
<ClInclude Include="include\astra\Float32Data2D.h" />
<ClInclude Include="include\astra\Float32Data3D.h" />
@@ -656,7 +658,6 @@
<ClInclude Include="include\astra\cuda\2d\fan_bp.h" />
<ClInclude Include="include\astra\cuda\2d\fan_fp.h" />
<ClInclude Include="include\astra\cuda\2d\fbp.h" />
- <ClInclude Include="include\astra\cuda\2d\fbp_filters.h" />
<ClInclude Include="include\astra\cuda\2d\fft.h" />
<ClInclude Include="include\astra\cuda\2d\par_bp.h" />
<ClInclude Include="include\astra\cuda\2d\par_fp.h" />
diff --git a/astra_vc14.vcxproj.filters b/astra_vc14.vcxproj.filters
index 3e9ff9a..fef6a8a 100644
--- a/astra_vc14.vcxproj.filters
+++ b/astra_vc14.vcxproj.filters
@@ -168,6 +168,9 @@
<ClCompile Include="src\Config.cpp">
<Filter>Global &amp; Other\source</Filter>
</ClCompile>
+ <ClCompile Include="src\Filters.cpp">
+ <Filter>Global &amp; Other\source</Filter>
+ </ClCompile>
<ClCompile Include="src\Fourier.cpp">
<Filter>Global &amp; Other\source</Filter>
</ClCompile>
@@ -434,6 +437,9 @@
<ClInclude Include="include\astra\Config.h">
<Filter>Global &amp; Other\headers</Filter>
</ClInclude>
+ <ClInclude Include="include\astra\Filters.h">
+ <Filter>Global &amp; Other\headers</Filter>
+ </ClInclude>
<ClInclude Include="include\astra\Fourier.h">
<Filter>Global &amp; Other\headers</Filter>
</ClInclude>
@@ -638,9 +644,6 @@
<ClInclude Include="include\astra\cuda\2d\fan_fp.h">
<Filter>CUDA\cuda headers</Filter>
</ClInclude>
- <ClInclude Include="include\astra\cuda\2d\fbp_filters.h">
- <Filter>CUDA\cuda headers</Filter>
- </ClInclude>
<ClInclude Include="include\astra\cuda\2d\fbp.h">
<Filter>CUDA\cuda headers</Filter>
</ClInclude>
diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in
index 1627a2e..0b90bd9 100644
--- a/build/linux/Makefile.in
+++ b/build/linux/Makefile.in
@@ -142,6 +142,7 @@ BASE_OBJECTS=\
src/FanFlatProjectionGeometry2D.lo \
src/FanFlatVecProjectionGeometry2D.lo \
src/FilteredBackProjectionAlgorithm.lo \
+ src/Filters.lo \
src/Float32Data2D.lo \
src/Float32Data3D.lo \
src/Float32Data3DMemory.lo \
diff --git a/build/msvc/gen.py b/build/msvc/gen.py
index fcc12d2..47e7440 100644
--- a/build/msvc/gen.py
+++ b/build/msvc/gen.py
@@ -214,6 +214,7 @@ P_astra["filters"]["Global &amp; Other\\source"] = [
"src\\AstraObjectManager.cpp",
"src\\CompositeGeometryManager.cpp",
"src\\Config.cpp",
+"src\\Filters.cpp",
"src\\Fourier.cpp",
"src\\Globals.cpp",
"src\\Logging.cpp",
@@ -292,7 +293,6 @@ P_astra["filters"]["CUDA\\cuda headers"] = [
"include\\astra\\cuda\\2d\\em.h",
"include\\astra\\cuda\\2d\\fan_bp.h",
"include\\astra\\cuda\\2d\\fan_fp.h",
-"include\\astra\\cuda\\2d\\fbp_filters.h",
"include\\astra\\cuda\\2d\\fbp.h",
"include\\astra\\cuda\\2d\\fft.h",
"include\\astra\\cuda\\2d\\par_bp.h",
@@ -354,6 +354,7 @@ P_astra["filters"]["Global &amp; Other\\headers"] = [
"include\\astra\\clog.h",
"include\\astra\\CompositeGeometryManager.h",
"include\\astra\\Config.h",
+"include\\astra\\Filters.h",
"include\\astra\\Fourier.h",
"include\\astra\\Globals.h",
"include\\astra\\Logging.h",
diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu
index 48fb7dc..a5b8a7a 100644
--- a/cuda/2d/fbp.cu
+++ b/cuda/2d/fbp.cu
@@ -35,27 +35,18 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#include "astra/cuda/3d/fdk.h"
#include "astra/Logging.h"
+#include "astra/Filters.h"
#include <cuda.h>
namespace astraCUDA {
-
-static int calcNextPowerOfTwo(int n)
-{
- int x = 1;
- while (x < n)
- x *= 2;
-
- return x;
-}
-
// static
int FBP::calcFourierFilterSize(int _iDetectorCount)
{
- int iFFTRealDetCount = calcNextPowerOfTwo(2 * _iDetectorCount);
- int iFreqBinCount = calcFFTFourierSize(iFFTRealDetCount);
+ int iFFTRealDetCount = astra::calcNextPowerOfTwo(2 * _iDetectorCount);
+ int iFreqBinCount = astra::calcFFTFourierSize(iFFTRealDetCount);
// CHECKME: Matlab makes this at least 64. Do we also need to?
return iFreqBinCount;
@@ -88,7 +79,7 @@ bool FBP::init()
return true;
}
-bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /* = NULL */, int _iFilterWidth /* = 0 */, float _fD /* = 1.0f */, float _fFilterParameter /* = -1.0f */)
+bool FBP::setFilter(const astra::SFilterConfig &_cfg)
{
if (D_filter)
{
@@ -96,19 +87,19 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
D_filter = 0;
}
- if (_eFilter == astra::FILTER_NONE)
+ if (_cfg.m_eType == astra::FILTER_NONE)
return true; // leave D_filter set to 0
- int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets);
- int iFreqBinCount = calcFFTFourierSize(iFFTRealDetCount);
+ int iFFTRealDetCount = astra::calcNextPowerOfTwo(2 * dims.iProjDets);
+ int iFreqBinCount = astra::calcFFTFourierSize(iFFTRealDetCount);
cufftComplex * pHostFilter = new cufftComplex[dims.iProjAngles * iFreqBinCount];
memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iFreqBinCount);
allocateComplexOnDevice(dims.iProjAngles, iFreqBinCount, (cufftComplex**)&D_filter);
- switch(_eFilter)
+ switch(_cfg.m_eType)
{
case astra::FILTER_NONE:
// handled above
@@ -130,7 +121,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(_cfg, dims.iProjAngles, pHostFilter, iFFTRealDetCount, iFreqBinCount);
uploadComplexArrayToDevice(dims.iProjAngles, iFreqBinCount, pHostFilter, (cufftComplex*)D_filter);
break;
@@ -138,11 +129,11 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
case astra::FILTER_PROJECTION:
{
// make sure the offered filter has the correct size
- assert(_iFilterWidth == iFreqBinCount);
+ assert(_cfg.m_iCustomFilterWidth == iFreqBinCount);
for(int iFreqBinIndex = 0; iFreqBinIndex < iFreqBinCount; iFreqBinIndex++)
{
- float fValue = _pfHostFilter[iFreqBinIndex];
+ float fValue = _cfg.m_pfCustomFilter[iFreqBinIndex];
for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++)
{
@@ -156,13 +147,13 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
case astra::FILTER_SINOGRAM:
{
// make sure the offered filter has the correct size
- assert(_iFilterWidth == iFreqBinCount);
+ assert(_cfg.m_iCustomFilterWidth == iFreqBinCount);
for(int iFreqBinIndex = 0; iFreqBinIndex < iFreqBinCount; iFreqBinIndex++)
{
for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++)
{
- float fValue = _pfHostFilter[iFreqBinIndex + iProjectionIndex * _iFilterWidth];
+ float fValue = _cfg.m_pfCustomFilter[iFreqBinIndex + iProjectionIndex * _cfg.m_iCustomFilterWidth];
pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].x = fValue;
pHostFilter[iFreqBinIndex + iProjectionIndex * iFreqBinCount].y = 0.0f;
@@ -178,16 +169,16 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
float * pfHostRealFilter = new float[iRealFilterElementCount];
memset(pfHostRealFilter, 0, sizeof(float) * iRealFilterElementCount);
- int iUsedFilterWidth = min(_iFilterWidth, iFFTRealDetCount);
- int iStartFilterIndex = (_iFilterWidth - iUsedFilterWidth) / 2;
+ int iUsedFilterWidth = min(_cfg.m_iCustomFilterWidth, iFFTRealDetCount);
+ int iStartFilterIndex = (_cfg.m_iCustomFilterWidth - iUsedFilterWidth) / 2;
int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth;
- int iFilterShiftSize = _iFilterWidth / 2;
+ int iFilterShiftSize = _cfg.m_iCustomFilterWidth / 2;
for(int iDetectorIndex = iStartFilterIndex; iDetectorIndex < iMaxFilterIndex; iDetectorIndex++)
{
int iFFTInFilterIndex = (iDetectorIndex + iFFTRealDetCount - iFilterShiftSize) % iFFTRealDetCount;
- float fValue = _pfHostFilter[iDetectorIndex];
+ float fValue = _cfg.m_pfCustomFilter[iDetectorIndex];
for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++)
{
@@ -213,11 +204,11 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
float* pfHostRealFilter = new float[iRealFilterElementCount];
memset(pfHostRealFilter, 0, sizeof(float) * iRealFilterElementCount);
- int iUsedFilterWidth = min(_iFilterWidth, iFFTRealDetCount);
- int iStartFilterIndex = (_iFilterWidth - iUsedFilterWidth) / 2;
+ int iUsedFilterWidth = min(_cfg.m_iCustomFilterWidth, iFFTRealDetCount);
+ int iStartFilterIndex = (_cfg.m_iCustomFilterWidth - iUsedFilterWidth) / 2;
int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth;
- int iFilterShiftSize = _iFilterWidth / 2;
+ int iFilterShiftSize = _cfg.m_iCustomFilterWidth / 2;
for(int iDetectorIndex = iStartFilterIndex; iDetectorIndex < iMaxFilterIndex; iDetectorIndex++)
{
@@ -225,7 +216,7 @@ bool FBP::setFilter(astra::E_FBPFILTER _eFilter, const float * _pfHostFilter /*
for(int iProjectionIndex = 0; iProjectionIndex < (int)dims.iProjAngles; iProjectionIndex++)
{
- float fValue = _pfHostFilter[iDetectorIndex + iProjectionIndex * _iFilterWidth];
+ float fValue = _cfg.m_pfCustomFilter[iDetectorIndex + iProjectionIndex * _cfg.m_iCustomFilterWidth];
pfHostRealFilter[iFFTInFilterIndex + iProjectionIndex * iFFTRealDetCount] = fValue;
}
}
@@ -310,8 +301,8 @@ bool FBP::iterate(unsigned int iterations)
if (D_filter) {
- int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets);
- int iFFTFourDetCount = calcFFTFourierSize(iFFTRealDetCount);
+ int iFFTRealDetCount = astra::calcNextPowerOfTwo(2 * dims.iProjDets);
+ int iFFTFourDetCount = astra::calcFFTFourierSize(iFFTRealDetCount);
cufftComplex * pDevComplexSinogram = NULL;
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;
}
diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu
index 46c07e7..1294721 100644
--- a/cuda/3d/fdk.cu
+++ b/cuda/3d/fdk.cu
@@ -246,14 +246,16 @@ bool FDK_Filter(cudaPitchedPtr D_projData,
// Generate filter
// TODO: Check errors
int iPaddedDetCount = calcNextPowerOfTwo(2 * dims.iProjU);
- int iHalfFFTSize = astraCUDA::calcFFTFourierSize(iPaddedDetCount);
+ int iHalfFFTSize = astra::calcFFTFourierSize(iPaddedDetCount);
cufftComplex *pHostFilter = new cufftComplex[dims.iProjAngles * iHalfFFTSize];
memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iHalfFFTSize);
if (pfFilter == 0){
- astraCUDA::genFilter(astra::FILTER_RAMLAK, 1.0f, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize);
+ astra::SFilterConfig filter;
+ filter.m_eType = astra::FILTER_RAMLAK;
+ astraCUDA::genCuFFTFilter(filter, dims.iProjAngles, pHostFilter, iPaddedDetCount, iHalfFFTSize);
} else {
for (int i = 0; i < dims.iProjAngles * iHalfFFTSize; i++) {
pHostFilter[i].x = pfFilter[i];
diff --git a/include/astra/CudaFilteredBackProjectionAlgorithm.h b/include/astra/CudaFilteredBackProjectionAlgorithm.h
index 1280e9a..8ef5a57 100644
--- a/include/astra/CudaFilteredBackProjectionAlgorithm.h
+++ b/include/astra/CudaFilteredBackProjectionAlgorithm.h
@@ -33,6 +33,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#include "Float32ProjectionData2D.h"
#include "Float32VolumeData2D.h"
#include "CudaReconstructionAlgorithm2D.h"
+#include "Filters.h"
#include "cuda/2d/astra.h"
@@ -45,15 +46,9 @@ public:
static std::string type;
private:
- E_FBPFILTER m_eFilter;
- float * m_pfFilter;
- int m_iFilterWidth; // number of elements per projection direction in filter
- float m_fFilterParameter; // some filters allow for parameterization (value < 0.0f -> no parameter)
- float m_fFilterD; // frequency cut-off
+ SFilterConfig m_filterConfig;
bool m_bShortScan; // short-scan mode for fan beam
- static E_FBPFILTER _convertStringToFilter(const char * _filterType);
-
public:
CCudaFilteredBackProjectionAlgorithm();
virtual ~CCudaFilteredBackProjectionAlgorithm();
diff --git a/include/astra/FilteredBackProjectionAlgorithm.h b/include/astra/FilteredBackProjectionAlgorithm.h
index 1cd4296..a234845 100644
--- a/include/astra/FilteredBackProjectionAlgorithm.h
+++ b/include/astra/FilteredBackProjectionAlgorithm.h
@@ -35,6 +35,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#include "Projector2D.h"
#include "Float32ProjectionData2D.h"
#include "Float32VolumeData2D.h"
+#include "Filters.h"
namespace astra {
@@ -144,6 +145,10 @@ public:
*/
virtual std::string description() const;
+protected:
+
+ SFilterConfig m_filterConfig;
+
};
// inline functions
diff --git a/include/astra/cuda/2d/fbp_filters.h b/include/astra/Filters.h
index 7c1121a..a1dec97 100644
--- a/include/astra/cuda/2d/fbp_filters.h
+++ b/include/astra/Filters.h
@@ -25,13 +25,18 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
-----------------------------------------------------------------------
*/
-#ifndef FBP_FILTERS_H
-#define FBP_FILTERS_H
+#ifndef _INC_ASTRA_FILTERS_H
+#define _INC_ASTRA_FILTERS_H
namespace astra {
+struct Config;
+class CAlgorithm;
+class CProjectionGeometry2D;
+
enum E_FBPFILTER
{
+ FILTER_ERROR, //< not a valid filter
FILTER_NONE, //< no filter (regular BP)
FILTER_RAMLAK, //< default FBP filter
FILTER_SHEPPLOGAN, //< Shepp-Logan
@@ -54,8 +59,40 @@ enum E_FBPFILTER
FILTER_SINOGRAM, //< every projection direction has its own filter
FILTER_RPROJECTION, //< projection filter in real space (as opposed to fourier space)
FILTER_RSINOGRAM, //< sinogram filter in real space
+
};
+struct SFilterConfig {
+ E_FBPFILTER m_eType;
+ float m_fD;
+ float m_fParameter;
+
+ float *m_pfCustomFilter;
+ int m_iCustomFilterWidth;
+ int m_iCustomFilterHeight;
+
+ SFilterConfig() : m_eType(FILTER_ERROR), m_fD(1.0f), m_fParameter(-1.0f),
+ m_pfCustomFilter(0), m_iCustomFilterWidth(0),
+ m_iCustomFilterHeight(0) { };
+};
+
+// Generate filter of given size and parameters. Returns newly allocated array.
+float *genFilter(const SFilterConfig &_cfg,
+ int _iFFTRealDetectorCount,
+ int _iFFTFourierDetectorCount);
+
+// Convert string to filter type. Returns FILTER_ERROR if unrecognized.
+E_FBPFILTER convertStringToFilter(const char * _filterType);
+
+SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg);
+
+bool checkCustomFilterSize(const SFilterConfig &_cfg, const CProjectionGeometry2D &_geom);
+
+
+int calcNextPowerOfTwo(int _iValue);
+int calcFFTFourierSize(int _iFFTRealSize);
+
+
}
-#endif /* FBP_FILTERS_H */
+#endif
diff --git a/include/astra/cuda/2d/astra.h b/include/astra/cuda/2d/astra.h
index 6f0e2f0..a2f2219 100644
--- a/include/astra/cuda/2d/astra.h
+++ b/include/astra/cuda/2d/astra.h
@@ -28,7 +28,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#ifndef _CUDA_ASTRA_H
#define _CUDA_ASTRA_H
-#include "fbp_filters.h"
#include "dims.h"
#include "algo.h"
diff --git a/include/astra/cuda/2d/fbp.h b/include/astra/cuda/2d/fbp.h
index 8666646..1adf3b1 100644
--- a/include/astra/cuda/2d/fbp.h
+++ b/include/astra/cuda/2d/fbp.h
@@ -26,7 +26,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
*/
#include "algo.h"
-#include "fbp_filters.h"
+#include "astra/Filters.h"
namespace astraCUDA {
@@ -75,9 +75,7 @@ public:
// FILTER_COSINE, FILTER_HAMMING, and FILTER_HANN)
// have a D variable, which gives the cutoff point in the frequency domain.
// Setting this value to 1.0 will include the whole filter
- bool setFilter(astra::E_FBPFILTER _eFilter,
- const float * _pfHostFilter = NULL,
- int _iFilterWidth = 0, float _fD = 1.0f, float _fFilterParameter = -1.0f);
+ bool setFilter(const astra::SFilterConfig &_cfg);
bool setShortScan(bool ss) { m_bShortScan = ss; return true; }
diff --git a/include/astra/cuda/2d/fft.h b/include/astra/cuda/2d/fft.h
index d36cae2..f77cbde 100644
--- a/include/astra/cuda/2d/fft.h
+++ b/include/astra/cuda/2d/fft.h
@@ -31,7 +31,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#include <cufft.h>
#include <cuda.h>
-#include "fbp_filters.h"
+#include "astra/Filters.h"
namespace astraCUDA {
@@ -58,11 +58,9 @@ bool runCudaIFFT(int _iProjectionCount, const cufftComplex* _pDevSourceComplex,
void applyFilter(int _iProjectionCount, int _iFreqBinCount,
cufftComplex * _pSinogram, cufftComplex * _pFilter);
-int calcFFTFourierSize(int _iFFTRealSize);
-
-void genFilter(astra::E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
- cufftComplex * _pFilter, int _iFFTRealDetectorCount,
- int _iFFTFourierDetectorCount, float _fParameter = -1.0f);
+void genCuFFTFilter(const astra::SFilterConfig &_cfg, int _iProjectionCount,
+ cufftComplex * _pFilter, int _iFFTRealDetectorCount,
+ int _iFFTFourierDetectorCount);
void genIdenFilter(int _iProjectionCount, cufftComplex * _pFilter,
int _iFFTRealDetectorCount, int _iFFTFourierDetectorCount);
diff --git a/samples/matlab/s023_FBP_filters.m b/samples/matlab/s023_FBP_filters.m
new file mode 100644
index 0000000..4abec7e
--- /dev/null
+++ b/samples/matlab/s023_FBP_filters.m
@@ -0,0 +1,96 @@
+% -----------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+%
+% Copyright: 2010-2018, imec Vision Lab, University of Antwerp
+% 2014-2018, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Contact: astra@astra-toolbox.com
+% Website: http://www.astra-toolbox.com/
+% -----------------------------------------------------------------------
+
+
+% This sample script illustrates three ways of passing filters to FBP.
+% They work with both the FBP (CPU) and the FBP_CUDA (GPU) algorithms.
+
+N = 256;
+
+vol_geom = astra_create_vol_geom(N, N);
+proj_geom = astra_create_proj_geom('parallel', 1.0, N, linspace2(0,pi,180));
+
+proj_id = astra_create_projector('strip', proj_geom, vol_geom);
+
+P = phantom(256);
+
+[sinogram_id, sinogram] = astra_create_sino(P, proj_id);
+
+rec_id = astra_mex_data2d('create', '-vol', vol_geom);
+
+cfg = astra_struct('FBP');
+cfg.ReconstructionDataId = rec_id;
+cfg.ProjectionDataId = sinogram_id;
+cfg.ProjectorId = proj_id;
+
+
+% 1. Use a standard Ram-Lak filter
+cfg.FilterType = 'ram-lak';
+
+alg_id = astra_mex_algorithm('create', cfg);
+astra_mex_algorithm('run', alg_id);
+rec_RL = astra_mex_data2d('get', rec_id);
+astra_mex_algorithm('delete', alg_id);
+
+
+% 2. Define a filter in Fourier space
+% This is assumed to be symmetric, and ASTRA therefore expects only half
+
+% The full filter size should be the smallest power of two that is at least
+% twice the number of detector pixels.
+fullFilterSize = 2*N;
+kernel = [linspace2(0, 1, floor(fullFilterSize / 2)) linspace2(1, 0, ceil(fullFilterSize / 2))];
+halfFilterSize = floor(fullFilterSize / 2) + 1;
+filter = kernel(1:halfFilterSize);
+
+filter_geom = astra_create_proj_geom('parallel', 1.0, halfFilterSize, [0]);
+filter_id = astra_mex_data2d('create', '-sino', filter_geom, filter);
+
+cfg.FilterType = 'projection';
+cfg.FilterSinogramId = filter_id;
+
+alg_id = astra_mex_algorithm('create', cfg);
+astra_mex_algorithm('run', alg_id);
+rec_filter = astra_mex_data2d('get', rec_id);
+astra_mex_algorithm('delete', alg_id);
+
+% 3. Define a (spatial) convolution kernel directly
+% For a kernel of odd size 2*k+1, the central component is at kernel(k+1)
+% For a kernel of even size 2*k, the central component is at kernel(k+1)
+
+kernel = zeros(1, N);
+for i = 0:floor(N/4)-1
+ f = pi * (2*i + 1);
+ val = -2.0 / (f * f);
+ kernel(floor(N/2) + 1 + (2*i+1)) = val;
+ kernel(floor(N/2) + 1 - (2*i+1)) = val;
+end
+kernel(floor(N/2)+1) = 0.5;
+
+kernel_geom = astra_create_proj_geom('parallel', 1.0, N, [0]);
+kernel_id = astra_mex_data2d('create', '-sino', kernel_geom, kernel);
+
+cfg.FilterType = 'rprojection';
+cfg.FilterSinogramId = kernel_id;
+
+alg_id = astra_mex_algorithm('create', cfg);
+astra_mex_algorithm('run', alg_id);
+rec_kernel = astra_mex_data2d('get', rec_id);
+astra_mex_algorithm('delete', alg_id);
+
+figure(1); imshow(P, []);
+figure(2); imshow(rec_RL, []);
+figure(3); imshow(rec_filter, []);
+figure(4); imshow(rec_kernel, []);
+
+
+astra_mex_data2d('delete', rec_id);
+astra_mex_data2d('delete', sinogram_id);
+astra_mex_projector('delete', proj_id);
diff --git a/samples/python/s023_FBP_filters.py b/samples/python/s023_FBP_filters.py
new file mode 100644
index 0000000..11518ac
--- /dev/null
+++ b/samples/python/s023_FBP_filters.py
@@ -0,0 +1,116 @@
+# -----------------------------------------------------------------------
+# Copyright: 2010-2018, imec Vision Lab, University of Antwerp
+# 2013-2018, CWI, Amsterdam
+#
+# Contact: astra@astra-toolbox.com
+# Website: http://www.astra-toolbox.com/
+#
+# This file is part of the ASTRA Toolbox.
+#
+#
+# The ASTRA Toolbox is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# The ASTRA Toolbox is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+#
+# -----------------------------------------------------------------------
+
+import astra
+import numpy as np
+import scipy.io
+
+
+# This sample script illustrates three ways of passing filters to FBP.
+# They work with both the FBP (CPU) and the FBP_CUDA (GPU) algorithms.
+
+
+N = 256
+
+vol_geom = astra.create_vol_geom(N, N)
+proj_geom = astra.create_proj_geom('parallel', 1.0, N, np.linspace(0,np.pi,180,False))
+
+P = scipy.io.loadmat('phantom.mat')['phantom256']
+
+proj_id = astra.create_projector('strip',proj_geom,vol_geom)
+sinogram_id, sinogram = astra.create_sino(P, proj_id)
+
+rec_id = astra.data2d.create('-vol', vol_geom)
+cfg = astra.astra_dict('FBP')
+cfg['ReconstructionDataId'] = rec_id
+cfg['ProjectionDataId'] = sinogram_id
+cfg['ProjectorId'] = proj_id
+
+
+
+# 1. Use a standard Ram-Lak filter
+cfg['FilterType'] = 'ram-lak'
+
+alg_id = astra.algorithm.create(cfg)
+astra.algorithm.run(alg_id)
+rec_RL = astra.data2d.get(rec_id)
+astra.algorithm.delete(alg_id)
+
+# 2. Define a filter in Fourier space
+# This is assumed to be symmetric, and ASTRA therefore expects only half
+
+# The full filter size should be the smallest power of two that is at least
+# twice the number of detector pixels.
+fullFilterSize = 2*N
+kernel = np.append( np.linspace(0, 1, fullFilterSize//2, endpoint=False), np.linspace(1, 0, fullFilterSize//2, endpoint=False) )
+halfFilterSize = fullFilterSize // 2 + 1
+filter = np.reshape(kernel[0:halfFilterSize], (1, halfFilterSize))
+
+filter_geom = astra.create_proj_geom('parallel', 1.0, halfFilterSize, [0]);
+filter_id = astra.data2d.create('-sino', filter_geom, filter);
+
+cfg['FilterType'] = 'projection'
+cfg['FilterSinogramId'] = filter_id
+alg_id = astra.algorithm.create(cfg)
+astra.algorithm.run(alg_id)
+rec_filter = astra.data2d.get(rec_id)
+astra.algorithm.delete(alg_id)
+
+
+# 3. Define a (spatial) convolution kernel directly
+# For a kernel of odd size 2*k+1, the central component is at kernel[k]
+# For a kernel of even size 2*k, the central component is at kernel[k]
+kernel = np.zeros((1, N))
+for i in range(0,N//4):
+ f = np.pi * (2*i + 1)
+ val = -2.0 / (f * f)
+ kernel[0, N//2 + (2*i+1)] = val
+ kernel[0, N//2 - (2*i+1)] = val
+kernel[0, N//2] = 0.5
+kernel_geom = astra.create_proj_geom('parallel', 1.0, N, [0]);
+kernel_id = astra.data2d.create('-sino', kernel_geom, kernel);
+
+cfg['FilterType'] = 'rprojection'
+cfg['FilterSinogramId'] = kernel_id
+alg_id = astra.algorithm.create(cfg)
+astra.algorithm.run(alg_id)
+rec_kernel = astra.data2d.get(rec_id)
+astra.algorithm.delete(alg_id)
+
+import pylab
+pylab.figure()
+pylab.imshow(P)
+pylab.figure()
+pylab.imshow(rec_RL)
+pylab.figure()
+pylab.imshow(rec_filter)
+pylab.figure()
+pylab.imshow(rec_kernel)
+pylab.show()
+
+astra.data2d.delete(rec_id)
+astra.data2d.delete(sinogram_id)
+astra.projector.delete(proj_id)
+
diff --git a/src/CudaFDKAlgorithm3D.cpp b/src/CudaFDKAlgorithm3D.cpp
index 4d8fc5a..24ed04f 100644
--- a/src/CudaFDKAlgorithm3D.cpp
+++ b/src/CudaFDKAlgorithm3D.cpp
@@ -35,9 +35,9 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#include "astra/CompositeGeometryManager.h"
#include "astra/Logging.h"
+#include "astra/Filters.h"
#include "astra/cuda/3d/astra3d.h"
-#include "astra/cuda/2d/fft.h"
#include "astra/cuda/3d/util3d.h"
using namespace std;
@@ -156,7 +156,7 @@ bool CCudaFDKAlgorithm3D::initialize(const Config& _cfg)
const CProjectionGeometry3D* projgeom = m_pSinogram->getGeometry();
const CProjectionGeometry2D* filtgeom = pFilterData->getGeometry();
int iPaddedDetCount = calcNextPowerOfTwo(2 * projgeom->getDetectorColCount());
- int iHalfFFTSize = astraCUDA::calcFFTFourierSize(iPaddedDetCount);
+ int iHalfFFTSize = calcFFTFourierSize(iPaddedDetCount);
if(filtgeom->getDetectorCount()!=iHalfFFTSize || filtgeom->getProjectionAngleCount()!=projgeom->getProjectionCount()){
ASTRA_ERROR("Filter size does not match required size (%i angles, %i detectors)",projgeom->getProjectionCount(),iHalfFFTSize);
return false;
diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp
index 944f165..88e235b 100644
--- a/src/CudaFilteredBackProjectionAlgorithm.cpp
+++ b/src/CudaFilteredBackProjectionAlgorithm.cpp
@@ -27,10 +27,10 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#include <astra/CudaFilteredBackProjectionAlgorithm.h>
#include <astra/FanFlatProjectionGeometry2D.h>
-#include <cstring>
#include "astra/AstraObjectManager.h"
#include "astra/CudaProjector2D.h"
+#include "astra/Filters.h"
#include "astra/cuda/2d/astra.h"
#include "astra/cuda/2d/fbp.h"
@@ -45,15 +45,12 @@ CCudaFilteredBackProjectionAlgorithm::CCudaFilteredBackProjectionAlgorithm()
{
m_bIsInitialized = false;
CCudaReconstructionAlgorithm2D::_clear();
- m_pfFilter = NULL;
- m_fFilterParameter = -1.0f;
- m_fFilterD = 1.0f;
}
CCudaFilteredBackProjectionAlgorithm::~CCudaFilteredBackProjectionAlgorithm()
{
- delete[] m_pfFilter;
- m_pfFilter = NULL;
+ delete[] m_filterConfig.m_pfCustomFilter;
+ m_filterConfig.m_pfCustomFilter = NULL;
}
bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
@@ -71,59 +68,7 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
if (!m_bIsInitialized)
return false;
-
- // filter type
- XMLNode node = _cfg.self.getSingleNode("FilterType");
- if (node)
- m_eFilter = _convertStringToFilter(node.getContent().c_str());
- else
- m_eFilter = FILTER_RAMLAK;
- CC.markNodeParsed("FilterType");
-
- // filter
- node = _cfg.self.getSingleNode("FilterSinogramId");
- if (node)
- {
- int id = node.getContentInt();
- const CFloat32ProjectionData2D * pFilterData = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id));
- m_iFilterWidth = pFilterData->getGeometry()->getDetectorCount();
- int iFilterProjectionCount = pFilterData->getGeometry()->getProjectionAngleCount();
-
- m_pfFilter = new float[m_iFilterWidth * iFilterProjectionCount];
- memcpy(m_pfFilter, pFilterData->getDataConst(), sizeof(float) * m_iFilterWidth * iFilterProjectionCount);
- }
- else
- {
- m_iFilterWidth = 0;
- m_pfFilter = NULL;
- }
- CC.markNodeParsed("FilterSinogramId"); // TODO: Only for some types!
-
- // filter parameter
- node = _cfg.self.getSingleNode("FilterParameter");
- if (node)
- {
- float fParameter = node.getContentNumerical();
- m_fFilterParameter = fParameter;
- }
- else
- {
- m_fFilterParameter = -1.0f;
- }
- CC.markNodeParsed("FilterParameter"); // TODO: Only for some types!
-
- // D value
- node = _cfg.self.getSingleNode("FilterD");
- if (node)
- {
- float fD = node.getContentNumerical();
- m_fFilterD = fD;
- }
- else
- {
- m_fFilterD = 1.0f;
- }
- CC.markNodeParsed("FilterD"); // TODO: Only for some types!
+ m_filterConfig = getFilterConfigForAlgorithm(_cfg, this);
// Fan beam short scan mode
if (m_pSinogram && dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pSinogram->getGeometry())) {
@@ -153,8 +98,8 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D *
m_pReconstruction = _pReconstruction;
m_iGPUIndex = _iGPUIndex;
- m_eFilter = _eFilter;
- m_iFilterWidth = _iFilterWidth;
+ m_filterConfig.m_eType = _eFilter;
+ m_filterConfig.m_iCustomFilterWidth = _iFilterWidth;
m_bShortScan = false;
// success
@@ -167,7 +112,7 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D *
{
int iFilterElementCount = 0;
- if((_eFilter != FILTER_SINOGRAM) && (_eFilter != FILTER_RSINOGRAM))
+ if((m_filterConfig.m_eType != FILTER_SINOGRAM) && (m_filterConfig.m_eType != FILTER_RSINOGRAM))
{
iFilterElementCount = _iFilterWidth;
}
@@ -176,15 +121,15 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D *
iFilterElementCount = m_pSinogram->getAngleCount();
}
- m_pfFilter = new float[iFilterElementCount];
- memcpy(m_pfFilter, _pfFilter, iFilterElementCount * sizeof(float));
+ m_filterConfig.m_pfCustomFilter = new float[iFilterElementCount];
+ memcpy(m_filterConfig.m_pfCustomFilter, _pfFilter, iFilterElementCount * sizeof(float));
}
else
{
- m_pfFilter = NULL;
+ m_filterConfig.m_pfCustomFilter = NULL;
}
- m_fFilterParameter = _fFilterParameter;
+ m_filterConfig.m_fParameter = _fFilterParameter;
return check();
}
@@ -196,7 +141,7 @@ void CCudaFilteredBackProjectionAlgorithm::initCUDAAlgorithm()
astraCUDA::FBP* pFBP = dynamic_cast<astraCUDA::FBP*>(m_pAlgo);
- bool ok = pFBP->setFilter(m_eFilter, m_pfFilter, m_iFilterWidth, m_fFilterD, m_fFilterParameter);
+ bool ok = pFBP->setFilter(m_filterConfig);
if (!ok) {
ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set filter");
ASTRA_ASSERT(ok);
@@ -215,9 +160,11 @@ bool CCudaFilteredBackProjectionAlgorithm::check()
ASTRA_CONFIG_CHECK(m_pSinogram, "FBP_CUDA", "Invalid Projection Data Object.");
ASTRA_CONFIG_CHECK(m_pReconstruction, "FBP_CUDA", "Invalid Reconstruction Data Object.");
- if((m_eFilter == FILTER_PROJECTION) || (m_eFilter == FILTER_SINOGRAM) || (m_eFilter == FILTER_RPROJECTION) || (m_eFilter == FILTER_RSINOGRAM))
+ ASTRA_CONFIG_CHECK(m_filterConfig.m_eType != FILTER_ERROR, "FBP_CUDA", "Invalid filter name.");
+
+ if((m_filterConfig.m_eType == FILTER_PROJECTION) || (m_filterConfig.m_eType == FILTER_SINOGRAM) || (m_filterConfig.m_eType == FILTER_RPROJECTION) || (m_filterConfig.m_eType == FILTER_RSINOGRAM))
{
- ASTRA_CONFIG_CHECK(m_pfFilter, "FBP_CUDA", "Invalid filter pointer.");
+ ASTRA_CONFIG_CHECK(m_filterConfig.m_pfCustomFilter, "FBP_CUDA", "Invalid filter pointer.");
}
// check initializations
@@ -229,122 +176,12 @@ bool CCudaFilteredBackProjectionAlgorithm::check()
// check pixel supersampling
ASTRA_CONFIG_CHECK(m_iPixelSuperSampling >= 0, "FBP_CUDA", "PixelSuperSampling must be a non-negative integer.");
+ ASTRA_CONFIG_CHECK(checkCustomFilterSize(m_filterConfig, *m_pSinogram->getGeometry()), "FBP_CUDA", "Filter size mismatch");
+
// success
m_bIsInitialized = true;
return true;
}
-static bool stringCompareLowerCase(const char * _stringA, const char * _stringB)
-{
- int iCmpReturn = 0;
-
-#ifdef _MSC_VER
- iCmpReturn = _stricmp(_stringA, _stringB);
-#else
- iCmpReturn = strcasecmp(_stringA, _stringB);
-#endif
-
- return (iCmpReturn == 0);
-}
-
-E_FBPFILTER CCudaFilteredBackProjectionAlgorithm::_convertStringToFilter(const char * _filterType)
-{
- E_FBPFILTER output = FILTER_NONE;
-
- if(stringCompareLowerCase(_filterType, "ram-lak"))
- {
- output = FILTER_RAMLAK;
- }
- else if(stringCompareLowerCase(_filterType, "shepp-logan"))
- {
- output = FILTER_SHEPPLOGAN;
- }
- else if(stringCompareLowerCase(_filterType, "cosine"))
- {
- output = FILTER_COSINE;
- }
- else if(stringCompareLowerCase(_filterType, "hamming"))
- {
- output = FILTER_HAMMING;
- }
- else if(stringCompareLowerCase(_filterType, "hann"))
- {
- output = FILTER_HANN;
- }
- else if(stringCompareLowerCase(_filterType, "none"))
- {
- output = FILTER_NONE;
- }
- else if(stringCompareLowerCase(_filterType, "tukey"))
- {
- output = FILTER_TUKEY;
- }
- else if(stringCompareLowerCase(_filterType, "lanczos"))
- {
- output = FILTER_LANCZOS;
- }
- else if(stringCompareLowerCase(_filterType, "triangular"))
- {
- output = FILTER_TRIANGULAR;
- }
- else if(stringCompareLowerCase(_filterType, "gaussian"))
- {
- output = FILTER_GAUSSIAN;
- }
- else if(stringCompareLowerCase(_filterType, "barlett-hann"))
- {
- output = FILTER_BARTLETTHANN;
- }
- else if(stringCompareLowerCase(_filterType, "blackman"))
- {
- output = FILTER_BLACKMAN;
- }
- else if(stringCompareLowerCase(_filterType, "nuttall"))
- {
- output = FILTER_NUTTALL;
- }
- else if(stringCompareLowerCase(_filterType, "blackman-harris"))
- {
- output = FILTER_BLACKMANHARRIS;
- }
- else if(stringCompareLowerCase(_filterType, "blackman-nuttall"))
- {
- output = FILTER_BLACKMANNUTTALL;
- }
- else if(stringCompareLowerCase(_filterType, "flat-top"))
- {
- output = FILTER_FLATTOP;
- }
- else if(stringCompareLowerCase(_filterType, "kaiser"))
- {
- output = FILTER_KAISER;
- }
- else if(stringCompareLowerCase(_filterType, "parzen"))
- {
- output = FILTER_PARZEN;
- }
- else if(stringCompareLowerCase(_filterType, "projection"))
- {
- output = FILTER_PROJECTION;
- }
- else if(stringCompareLowerCase(_filterType, "sinogram"))
- {
- output = FILTER_SINOGRAM;
- }
- else if(stringCompareLowerCase(_filterType, "rprojection"))
- {
- output = FILTER_RPROJECTION;
- }
- else if(stringCompareLowerCase(_filterType, "rsinogram"))
- {
- output = FILTER_RSINOGRAM;
- }
- else
- {
- ASTRA_ERROR("Failed to convert \"%s\" into a filter.",_filterType);
- }
-
- return output;
-}
diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp
index bb2e722..67a12a2 100644
--- a/src/FilteredBackProjectionAlgorithm.cpp
+++ b/src/FilteredBackProjectionAlgorithm.cpp
@@ -81,6 +81,9 @@ void CFilteredBackProjectionAlgorithm::clear()
m_pSinogram = NULL;
m_pReconstruction = NULL;
m_bIsInitialized = false;
+
+ delete[] m_filterConfig.m_pfCustomFilter;
+ m_filterConfig.m_pfCustomFilter = 0;
}
@@ -151,6 +154,9 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
delete[] projectionAngles;
}
+ m_filterConfig = getFilterConfigForAlgorithm(_cfg, this);
+
+
// TODO: check that the angles are linearly spaced between 0 and pi
// success
@@ -195,11 +201,14 @@ bool CFilteredBackProjectionAlgorithm::initialize(CProjector2D* _pProjector,
CFloat32VolumeData2D* _pVolume,
CFloat32ProjectionData2D* _pSinogram)
{
+ clear();
+
// store classes
m_pProjector = _pProjector;
m_pReconstruction = _pVolume;
m_pSinogram = _pSinogram;
+ m_filterConfig = SFilterConfig();
// TODO: check that the angles are linearly spaced between 0 and pi
@@ -214,6 +223,15 @@ bool CFilteredBackProjectionAlgorithm::_check()
{
ASTRA_CONFIG_CHECK(CReconstructionAlgorithm2D::_check(), "FBP", "Error in ReconstructionAlgorithm2D initialization");
+ ASTRA_CONFIG_CHECK(m_filterConfig.m_eType != FILTER_ERROR, "FBP", "Invalid filter name.");
+
+ if((m_filterConfig.m_eType == FILTER_PROJECTION) || (m_filterConfig.m_eType == FILTER_SINOGRAM) || (m_filterConfig.m_eType == FILTER_RPROJECTION) || (m_filterConfig.m_eType == FILTER_RSINOGRAM))
+ {
+ ASTRA_CONFIG_CHECK(m_filterConfig.m_pfCustomFilter, "FBP", "Invalid filter pointer.");
+ }
+
+ ASTRA_CONFIG_CHECK(checkCustomFilterSize(m_filterConfig, *m_pSinogram->getGeometry()), "FBP", "Filter size mismatch");
+
// success
return true;
}
@@ -249,34 +267,84 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D
ASTRA_ASSERT(_pFilteredSinogram->getAngleCount() == m_pSinogram->getAngleCount());
ASTRA_ASSERT(_pFilteredSinogram->getDetectorCount() == m_pSinogram->getDetectorCount());
+ ASTRA_ASSERT(m_filterConfig.m_eType != FILTER_ERROR);
+ if (m_filterConfig.m_eType == FILTER_NONE)
+ return;
int iAngleCount = m_pProjector->getProjectionGeometry()->getProjectionAngleCount();
int iDetectorCount = m_pProjector->getProjectionGeometry()->getDetectorCount();
- // We'll zero-pad to the smallest power of two at least 64 and
- // at least 2*iDetectorCount
- int zpDetector = 64;
- int nextPow2 = 5;
- while (zpDetector < iDetectorCount*2) {
- zpDetector *= 2;
- nextPow2++;
- }
+ int zpDetector = calcNextPowerOfTwo(2 * m_pSinogram->getDetectorCount());
+ int iHalfFFTSize = astra::calcFFTFourierSize(zpDetector);
- // Create filter
- float32* filter = new float32[zpDetector];
+ // cdft setup
+ int *ip = new int[int(2+sqrt((float)zpDetector)+1)];
+ ip[0] = 0;
+ float32 *w = new float32[zpDetector/2];
- for (int iDetector = 0; iDetector <= zpDetector/2; iDetector++)
- filter[iDetector] = (2.0f * iDetector)/zpDetector;
+ // Create filter
+ bool bFilterMultiAngle = false;
+ bool bFilterComplex = false;
+ float *pfFilter = 0;
+ float *pfFilter_delete = 0;
+ switch (m_filterConfig.m_eType) {
+ case FILTER_ERROR:
+ case FILTER_NONE:
+ // Should have been handled before
+ ASTRA_ASSERT(false);
+ return;
+ case FILTER_PROJECTION:
+ // Fourier space, real, half the coefficients (because symmetric)
+ // 1 x iHalfFFTSize
+ pfFilter = m_filterConfig.m_pfCustomFilter;
+ break;
+ case FILTER_SINOGRAM:
+ bFilterMultiAngle = true;
+ pfFilter = m_filterConfig.m_pfCustomFilter;
+ break;
+ case FILTER_RSINOGRAM:
+ bFilterMultiAngle = true;
+ // fall-through
+ case FILTER_RPROJECTION:
+ {
+ bFilterComplex = true;
+
+ int count = bFilterMultiAngle ? iAngleCount : 1;
+ // Spatial, real, full convolution kernel
+ // Center in center (or right-of-center for even sized.)
+ // I.e., 0 1 0 and 0 0 1 0 both correspond to the identity
+
+ pfFilter = new float[2 * zpDetector * count];
+ pfFilter_delete = pfFilter;
+
+ int iUsedFilterWidth = min(m_filterConfig.m_iCustomFilterWidth, zpDetector);
+ int iStartFilterIndex = (m_filterConfig.m_iCustomFilterWidth - iUsedFilterWidth) / 2;
+ int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth;
+
+ int iFilterShiftSize = m_filterConfig.m_iCustomFilterWidth / 2;
+
+ for (int i = 0; i < count; ++i) {
+ float *rOut = pfFilter + i * 2 * zpDetector;
+ float *rIn = m_filterConfig.m_pfCustomFilter + i * m_filterConfig.m_iCustomFilterWidth;
+ memset(rOut, 0, sizeof(float) * 2 * zpDetector);
+
+ for(int j = iStartFilterIndex; j < iMaxFilterIndex; j++) {
+ int iFFTInFilterIndex = (j + zpDetector - iFilterShiftSize) % zpDetector;
+ rOut[2 * iFFTInFilterIndex] = rIn[j];
+ }
- for (int iDetector = zpDetector/2+1; iDetector < zpDetector; iDetector++)
- filter[iDetector] = (2.0f * (zpDetector - iDetector)) / zpDetector;
+ cdft(2*zpDetector, -1, rOut, ip, w);
+ }
+ break;
+ }
+ default:
+ pfFilter = genFilter(m_filterConfig, zpDetector, iHalfFFTSize);
+ pfFilter_delete = pfFilter;
+ }
float32* pf = new float32[2 * iAngleCount * zpDetector];
- int *ip = new int[int(2+sqrt((float)zpDetector)+1)];
- ip[0]=0;
- float32 *w = new float32[zpDetector/2];
// Copy and zero-pad data
for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) {
@@ -299,11 +367,34 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D
}
// Filter
- for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) {
- float32* pfRow = pf + iAngle * 2 * zpDetector;
- for (int iDetector = 0; iDetector < zpDetector; ++iDetector) {
- pfRow[2*iDetector] *= filter[iDetector];
- pfRow[2*iDetector+1] *= filter[iDetector];
+ if (bFilterComplex) {
+ for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) {
+ float32* pfRow = pf + iAngle * 2 * zpDetector;
+ float *pfFilterRow = pfFilter;
+ if (bFilterMultiAngle)
+ pfFilterRow += iAngle * 2 * zpDetector;
+
+ for (int i = 0; i < zpDetector; ++i) {
+ float re = pfRow[2*i] * pfFilterRow[2*i] - pfRow[2*i+1] * pfFilterRow[2*i+1];
+ float im = pfRow[2*i] * pfFilterRow[2*i+1] + pfRow[2*i+1] * pfFilterRow[2*i];
+ pfRow[2*i] = re;
+ pfRow[2*i+1] = im;
+ }
+ }
+ } else {
+ for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) {
+ float32* pfRow = pf + iAngle * 2 * zpDetector;
+ float *pfFilterRow = pfFilter;
+ if (bFilterMultiAngle)
+ pfFilterRow += iAngle * iHalfFFTSize;
+ for (int iDetector = 0; iDetector < iHalfFFTSize; ++iDetector) {
+ pfRow[2*iDetector] *= pfFilterRow[iDetector];
+ pfRow[2*iDetector+1] *= pfFilterRow[iDetector];
+ }
+ for (int iDetector = iHalfFFTSize; iDetector < zpDetector; ++iDetector) {
+ pfRow[2*iDetector] *= pfFilterRow[zpDetector - iDetector];
+ pfRow[2*iDetector+1] *= pfFilterRow[zpDetector - iDetector];
+ }
}
}
@@ -324,7 +415,7 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D
delete[] pf;
delete[] w;
delete[] ip;
- delete[] filter;
+ delete[] pfFilter_delete;
}
}
diff --git a/src/Filters.cpp b/src/Filters.cpp
new file mode 100644
index 0000000..c13aa6b
--- /dev/null
+++ b/src/Filters.cpp
@@ -0,0 +1,608 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2018, imec Vision Lab, University of Antwerp
+ 2014-2018, CWI, Amsterdam
+
+Contact: astra@astra-toolbox.com
+Website: http://www.astra-toolbox.com/
+
+This file is part of the ASTRA Toolbox.
+
+
+The ASTRA Toolbox is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+The ASTRA Toolbox is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+*/
+
+#include "astra/Globals.h"
+#include "astra/Logging.h"
+#include "astra/Fourier.h"
+#include "astra/Filters.h"
+#include "astra/Config.h"
+#include "astra/Float32ProjectionData2D.h"
+#include "astra/AstraObjectManager.h"
+
+#include <utility>
+#include <cstring>
+
+namespace astra {
+
+float *genFilter(const SFilterConfig &_cfg,
+ int _iFFTRealDetectorCount,
+ 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(_cfg.m_eType)
+ {
+ 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 / _cfg.m_fD) / (pfW[iDetectorIndex] / 2.0f / _cfg.m_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 / _cfg.m_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] / _cfg.m_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] / _cfg.m_fD)) / 2.0f;
+ }
+ break;
+ }
+ case FILTER_TUKEY:
+ {
+ float fAlpha = _cfg.m_fParameter;
+ if(_cfg.m_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 = _cfg.m_fParameter;
+ if(_cfg.m_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 = _cfg.m_fParameter;
+ if(_cfg.m_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 = _cfg.m_fParameter;
+ if(_cfg.m_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 * _cfg.m_fD;
+ for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fWValue = pfW[iDetectorIndex];
+
+ if(fWValue > fPiTimesD)
+ {
+ pfFilt[iDetectorIndex] = 0.0f;
+ }
+ }
+
+ delete[] pfW;
+
+ return pfFilt;
+}
+
+static bool stringCompareLowerCase(const char * _stringA, const char * _stringB)
+{
+ int iCmpReturn = 0;
+
+#ifdef _MSC_VER
+ iCmpReturn = _stricmp(_stringA, _stringB);
+#else
+ iCmpReturn = strcasecmp(_stringA, _stringB);
+#endif
+
+ return (iCmpReturn == 0);
+}
+
+struct FilterNameMapEntry {
+ const char *m_name;
+ E_FBPFILTER m_type;
+};
+
+E_FBPFILTER convertStringToFilter(const char * _filterType)
+{
+
+ static const FilterNameMapEntry map[] = {
+ { "ram-lak", FILTER_RAMLAK },
+ { "shepp-logan", FILTER_SHEPPLOGAN },
+ { "cosine", FILTER_COSINE },
+ { "hamming", FILTER_HAMMING},
+ { "hann", FILTER_HANN},
+ { "tukey", FILTER_TUKEY },
+ { "lanczos", FILTER_LANCZOS},
+ { "triangular", FILTER_TRIANGULAR},
+ { "gaussian", FILTER_GAUSSIAN},
+ { "barlett-hann", FILTER_BARTLETTHANN },
+ { "blackman", FILTER_BLACKMAN},
+ { "nuttall", FILTER_NUTTALL },
+ { "blackman-harris", FILTER_BLACKMANHARRIS },
+ { "blackman-nuttall", FILTER_BLACKMANNUTTALL },
+ { "flat-top", FILTER_FLATTOP },
+ { "kaiser", FILTER_KAISER },
+ { "parzen", FILTER_PARZEN },
+ { "projection", FILTER_PROJECTION },
+ { "sinogram", FILTER_SINOGRAM },
+ { "rprojection", FILTER_RPROJECTION },
+ { "rsinogram", FILTER_RSINOGRAM },
+ { "none", FILTER_NONE },
+ { 0, FILTER_ERROR } };
+
+ const FilterNameMapEntry *i;
+
+ for (i = &map[0]; i->m_name; ++i)
+ if (stringCompareLowerCase(_filterType, i->m_name))
+ return i->m_type;
+
+ ASTRA_ERROR("Failed to convert \"%s\" into a filter.",_filterType);
+
+ return FILTER_ERROR;
+}
+
+
+SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg)
+{
+ ConfigStackCheck<CAlgorithm> CC("getFilterConfig", _alg, _cfg);
+
+ SFilterConfig c;
+
+ // filter type
+ XMLNode node = _cfg.self.getSingleNode("FilterType");
+ if (node)
+ c.m_eType = convertStringToFilter(node.getContent().c_str());
+ else
+ c.m_eType = FILTER_RAMLAK;
+ CC.markNodeParsed("FilterType");
+
+ // filter
+ node = _cfg.self.getSingleNode("FilterSinogramId");
+ if (node)
+ {
+ int id = node.getContentInt();
+ const CFloat32ProjectionData2D * pFilterData = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id));
+ c.m_iCustomFilterWidth = pFilterData->getGeometry()->getDetectorCount();
+ c.m_iCustomFilterHeight = pFilterData->getGeometry()->getProjectionAngleCount();
+
+ c.m_pfCustomFilter = new float[c.m_iCustomFilterWidth * c.m_iCustomFilterHeight];
+ memcpy(c.m_pfCustomFilter, pFilterData->getDataConst(), sizeof(float) * c.m_iCustomFilterWidth * c.m_iCustomFilterHeight);
+ }
+ else
+ {
+ c.m_iCustomFilterWidth = 0;
+ c.m_iCustomFilterHeight = 0;
+ c.m_pfCustomFilter = NULL;
+ }
+ CC.markNodeParsed("FilterSinogramId"); // TODO: Only for some types!
+
+ // filter parameter
+ node = _cfg.self.getSingleNode("FilterParameter");
+ if (node)
+ {
+ float fParameter = node.getContentNumerical();
+ c.m_fParameter = fParameter;
+ }
+ else
+ {
+ c.m_fParameter = -1.0f;
+ }
+ CC.markNodeParsed("FilterParameter"); // TODO: Only for some types!
+
+ // D value
+ node = _cfg.self.getSingleNode("FilterD");
+ if (node)
+ {
+ float fD = node.getContentNumerical();
+ c.m_fD = fD;
+ }
+ else
+ {
+ c.m_fD = 1.0f;
+ }
+ CC.markNodeParsed("FilterD"); // TODO: Only for some types!
+
+ return c;
+}
+
+int calcNextPowerOfTwo(int n)
+{
+ int x = 1;
+ while (x < n && x > 0)
+ x *= 2;
+
+ return x;
+}
+
+// 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;
+}
+
+bool checkCustomFilterSize(const SFilterConfig &_cfg, const CProjectionGeometry2D &_geom) {
+ int iExpectedWidth = -1, iExpectedHeight = 1;
+
+ switch (_cfg.m_eType) {
+ case FILTER_ERROR:
+ ASTRA_ERROR("checkCustomFilterSize: internal error; FILTER_ERROR passed");
+ return false;
+ case FILTER_NONE:
+ return true;
+ case FILTER_SINOGRAM:
+ iExpectedHeight = _geom.getProjectionAngleCount();
+ // fallthrough
+ case FILTER_PROJECTION:
+ {
+ int iPaddedDetCount = calcNextPowerOfTwo(2 * _geom.getDetectorCount());
+ iExpectedWidth = calcFFTFourierSize(iPaddedDetCount);
+ }
+ if (_cfg.m_iCustomFilterWidth != iExpectedWidth ||
+ _cfg.m_iCustomFilterHeight != iExpectedHeight)
+ {
+ ASTRA_ERROR("filter size mismatch: %dx%d (received) is not %dx%d (expected)", _cfg.m_iCustomFilterHeight, _cfg.m_iCustomFilterWidth, iExpectedHeight, iExpectedWidth);
+ return false;
+ }
+ return true;
+ case FILTER_RSINOGRAM:
+ iExpectedHeight = _geom.getProjectionAngleCount();
+ // fallthrough
+ case FILTER_RPROJECTION:
+ if (_cfg.m_iCustomFilterHeight != iExpectedHeight)
+ {
+ ASTRA_ERROR("filter size mismatch: %dx%d (received) is not %dxX (expected)", _cfg.m_iCustomFilterHeight, _cfg.m_iCustomFilterWidth, iExpectedHeight);
+ return false;
+ }
+ return true;
+ default:
+ // Non-custom filters; nothing to check.
+ return true;
+ }
+}
+
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-07 19:27:28 +0000