From da434892133fe0979c5e8ef8be277673452e7728 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 13 Jul 2018 14:48:19 +0200 Subject: Add filter size error reporting --- cuda/2d/fbp.cu | 7 +-- cuda/2d/fft.cu | 11 ----- cuda/3d/fdk.cu | 2 +- include/astra/Filters.h | 10 ++++- include/astra/cuda/2d/fft.h | 2 - src/CudaFDKAlgorithm3D.cpp | 4 +- src/CudaFilteredBackProjectionAlgorithm.cpp | 2 + src/Filters.cpp | 66 +++++++++++++++++++++++++++-- 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 . #include "astra/cuda/3d/fdk.h" #include "astra/Logging.h" +#include "astra/Filters.h" #include @@ -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 . #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(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