From da434892133fe0979c5e8ef8be277673452e7728 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Fri, 13 Jul 2018 14:48:19 +0200
Subject: Add filter size error reporting

---
 src/CudaFDKAlgorithm3D.cpp                  |  4 +-
 src/CudaFilteredBackProjectionAlgorithm.cpp |  2 +
 src/Filters.cpp                             | 66 +++++++++++++++++++++++++++--
 3 files changed, 67 insertions(+), 5 deletions(-)

(limited to 'src')

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;
+	}
+}
+
 }
-- 
cgit v1.2.3