Add filter size error reporting

8 files changed, 81 insertions, 23 deletions

diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu
index 1574ccc..e8a224e 100644
--- a/cuda/2d/fbp.cu
+++ b/cuda/2d/fbp.cu
@@ -35,6 +35,7 @@ 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>
 
@@ -55,7 +56,7 @@ static int calcNextPowerOfTwo(int n)
 int FBP::calcFourierFilterSize(int _iDetectorCount)
 {
 	int iFFTRealDetCount = calcNextPowerOfTwo(2 * _iDetectorCount);
-	int iFreqBinCount = calcFFTFourierSize(iFFTRealDetCount);
+	int iFreqBinCount = astra::calcFFTFourierSize(iFFTRealDetCount);
 
 	// CHECKME: Matlab makes this at least 64. Do we also need to?
 	return iFreqBinCount;
@@ -101,7 +102,7 @@ bool FBP::setFilter(const astra::SFilterConfig &_cfg)
 
 
 	int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets);
-	int iFreqBinCount = calcFFTFourierSize(iFFTRealDetCount);
+	int iFreqBinCount = astra::calcFFTFourierSize(iFFTRealDetCount);
 
 	cufftComplex * pHostFilter = new cufftComplex[dims.iProjAngles * iFreqBinCount];
 	memset(pHostFilter, 0, sizeof(cufftComplex) * dims.iProjAngles * iFreqBinCount);
@@ -311,7 +312,7 @@ bool FBP::iterate(unsigned int iterations)
 	if (D_filter) {
 
 		int iFFTRealDetCount = calcNextPowerOfTwo(2 * dims.iProjDets);
-		int iFFTFourDetCount = calcFFTFourierSize(iFFTRealDetCount);
+		int iFFTFourDetCount = astra::calcFFTFourierSize(iFFTRealDetCount);
 
 		cufftComplex * pDevComplexSinogram = NULL;
 
diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu
index 864e325..dc0edc3 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)
 {
diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu
index 014529b..1294721 100644
--- a/cuda/3d/fdk.cu
+++ b/cuda/3d/fdk.cu
@@ -246,7 +246,7 @@ 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];
diff --git a/include/astra/Filters.h b/include/astra/Filters.h
index eec2ba2..bf0fb48 100644
--- a/include/astra/Filters.h
+++ b/include/astra/Filters.h
@@ -32,6 +32,7 @@ namespace astra {
 
 struct Config;
 class CAlgorithm;
+class CProjectionGeometry2D;
 
 enum E_FBPFILTER
 {
@@ -68,9 +69,11 @@ struct SFilterConfig {
 
 	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_pfCustomFilter(0), m_iCustomFilterWidth(0),
+	                  m_iCustomFilterHeight(0) { };
 };
 
 // Generate filter of given size and parameters. Returns newly allocated array.
@@ -83,6 +86,11 @@ E_FBPFILTER convertStringToFilter(const char * _filterType);
 
 SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg);
 
+bool checkCustomFilterSize(const SFilterConfig &_cfg, const CProjectionGeometry2D &_geom);
+
+int calcFFTFourierSize(int _iFFTRealSize);
+
+
 }
 
 #endif
diff --git a/include/astra/cuda/2d/fft.h b/include/astra/cuda/2d/fft.h
index 33612a0..f77cbde 100644
--- a/include/astra/cuda/2d/fft.h
+++ b/include/astra/cuda/2d/fft.h
@@ -58,8 +58,6 @@ bool runCudaIFFT(int _iProjectionCount, const cufftComplex* _pDevSourceComplex,
 void applyFilter(int _iProjectionCount, int _iFreqBinCount,
                  cufftComplex * _pSinogram, cufftComplex * _pFilter);
 
-int calcFFTFourierSize(int _iFFTRealSize);
-
 void genCuFFTFilter(const astra::SFilterConfig &_cfg, int _iProjectionCount,
                    cufftComplex * _pFilter, int _iFFTRealDetectorCount,
                    int _iFFTFourierDetectorCount);
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 90b831e..88e235b 100644
--- a/src/CudaFilteredBackProjectionAlgorithm.cpp
+++ b/src/CudaFilteredBackProjectionAlgorithm.cpp
@@ -176,6 +176,8 @@ 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;
diff --git a/src/Filters.cpp b/src/Filters.cpp
index bbd7cfd..3ee3749 100644
--- a/src/Filters.cpp
+++ b/src/Filters.cpp
@@ -504,14 +504,15 @@ SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg)
 		int id = node.getContentInt();
 		const CFloat32ProjectionData2D * pFilterData = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id));
 		c.m_iCustomFilterWidth = pFilterData->getGeometry()->getDetectorCount();
-		int iFilterProjectionCount = pFilterData->getGeometry()->getProjectionAngleCount();
+		c.m_iCustomFilterHeight = pFilterData->getGeometry()->getProjectionAngleCount();
 
-		c.m_pfCustomFilter = new float[c.m_iCustomFilterWidth * iFilterProjectionCount];
-		memcpy(c.m_pfCustomFilter, pFilterData->getDataConst(), sizeof(float) * c.m_iCustomFilterWidth * iFilterProjectionCount);
+		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!
@@ -545,4 +546,63 @@ SFilterConfig getFilterConfigForAlgorithm(const Config& _cfg, CAlgorithm *_alg)
 	return c;
 }
 
+static int calcNextPowerOfTwo(int n)
+{
+	int x = 1;
+	while (x < n)
+		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;
+	}
+}
+
 }