Merge pull request #160 from wjp/FBP_filters

Custom filter support for CPU FBP

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;
+	}
+}
+
+}