summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--cuda/2d/fbp.cu7
-rw-r--r--cuda/2d/fft.cu11
-rw-r--r--cuda/3d/fdk.cu2
-rw-r--r--include/astra/Filters.h10
-rw-r--r--include/astra/cuda/2d/fft.h2
-rw-r--r--src/CudaFDKAlgorithm3D.cpp4
-rw-r--r--src/CudaFilteredBackProjectionAlgorithm.cpp2
-rw-r--r--src/Filters.cpp66
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;
+ }
+}
+
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-01 06:29:54 +0000