Merge pull request #160 from wjp/FBP_filters

Custom filter support for CPU FBP

4 files changed, 742 insertions, 206 deletions

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