summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
authorWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2018-07-12 12:11:52 +0200
committerWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2018-07-17 11:30:05 +0200
commitd32df3bcd1910f56195e828a0f7fba8fc04b90ab (patch)
tree56963d80f825af10639dee0e1c7c439bbffd1f76 /src
parentd58f6b51493c7931cbc02feb3ffeb0eca6ea3a4e (diff)
downloadastra-d32df3bcd1910f56195e828a0f7fba8fc04b90ab.tar.gz
astra-d32df3bcd1910f56195e828a0f7fba8fc04b90ab.tar.bz2
astra-d32df3bcd1910f56195e828a0f7fba8fc04b90ab.tar.xz
astra-d32df3bcd1910f56195e828a0f7fba8fc04b90ab.zip
Refactor filter config
Diffstat (limited to 'src')
-rw-r--r--src/CudaFilteredBackProjectionAlgorithm.cpp83
-rw-r--r--src/Filters.cpp96
2 files changed, 94 insertions, 85 deletions
diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp
index 9ebbd60..90b831e 100644
--- a/src/CudaFilteredBackProjectionAlgorithm.cpp
+++ b/src/CudaFilteredBackProjectionAlgorithm.cpp
@@ -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,11 +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.");
- ASTRA_CONFIG_CHECK(m_eFilter != FILTER_ERROR, "FBP_CUDA", "Invalid filter name.");
+ ASTRA_CONFIG_CHECK(m_filterConfig.m_eType != FILTER_ERROR, "FBP_CUDA", "Invalid filter name.");
- if((m_eFilter == FILTER_PROJECTION) || (m_eFilter == FILTER_SINOGRAM) || (m_eFilter == FILTER_RPROJECTION) || (m_eFilter == FILTER_RSINOGRAM))
+ 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
diff --git a/src/Filters.cpp b/src/Filters.cpp
index 30b2f24..bbd7cfd 100644
--- a/src/Filters.cpp
+++ b/src/Filters.cpp
@@ -29,15 +29,18 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#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(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
+float *genFilter(const SFilterConfig &_cfg, int _iProjectionCount,
int _iFFTRealDetectorCount,
- int _iFFTFourierDetectorCount, float _fParameter /* = -1.0f */)
+ int _iFFTFourierDetectorCount)
{
float * pfFilt = new float[_iFFTFourierDetectorCount];
float * pfW = new float[_iFFTFourierDetectorCount];
@@ -86,7 +89,7 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex;
}
- switch(_eFilter)
+ switch(_cfg.m_eType)
{
case FILTER_RAMLAK:
{
@@ -98,7 +101,7 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
// 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 / _fD) / (pfW[iDetectorIndex] / 2.0f / _fD));
+ pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (sinf(pfW[iDetectorIndex] / 2.0f / _cfg.m_fD) / (pfW[iDetectorIndex] / 2.0f / _cfg.m_fD));
}
break;
}
@@ -107,7 +110,7 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
// 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 / _fD);
+ pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * cosf(pfW[iDetectorIndex] / 2.0f / _cfg.m_fD);
}
break;
}
@@ -116,7 +119,7 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
// 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] / _fD));
+ pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * ( 0.54f + 0.46f * cosf(pfW[iDetectorIndex] / _cfg.m_fD));
}
break;
}
@@ -125,14 +128,14 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
// 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] / _fD)) / 2.0f;
+ pfFilt[iDetectorIndex] = pfFilt[iDetectorIndex] * (1.0f + cosf(pfW[iDetectorIndex] / _cfg.m_fD)) / 2.0f;
}
break;
}
case FILTER_TUKEY:
{
- float fAlpha = _fParameter;
- if(_fParameter < 0.0f) fAlpha = 0.5f;
+ 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;
@@ -204,8 +207,8 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
}
case FILTER_GAUSSIAN:
{
- float fSigma = _fParameter;
- if(_fParameter < 0.0f) fSigma = 0.4f;
+ 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;
@@ -245,8 +248,8 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
}
case FILTER_BLACKMAN:
{
- float fAlpha = _fParameter;
- if(_fParameter < 0.0f) fAlpha = 0.16f;
+ 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;
@@ -356,8 +359,8 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
}
case FILTER_KAISER:
{
- float fAlpha = _fParameter;
- if(_fParameter < 0.0f) fAlpha = 3.0f;
+ 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);
@@ -406,7 +409,7 @@ float *genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
}
// filt(w>pi*d) = 0;
- float fPiTimesD = M_PI * _fD;
+ float fPiTimesD = M_PI * _cfg.m_fD;
for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
{
float fWValue = pfW[iDetectorIndex];
@@ -480,5 +483,66 @@ E_FBPFILTER convertStringToFilter(const char * _filterType)
}
+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();
+ int iFilterProjectionCount = pFilterData->getGeometry()->getProjectionAngleCount();
+
+ c.m_pfCustomFilter = new float[c.m_iCustomFilterWidth * iFilterProjectionCount];
+ memcpy(c.m_pfCustomFilter, pFilterData->getDataConst(), sizeof(float) * c.m_iCustomFilterWidth * iFilterProjectionCount);
+ }
+ else
+ {
+ c.m_iCustomFilterWidth = 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;
+}
}