summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/CudaFDKAlgorithm3D.cpp4
-rw-r--r--src/CudaFilteredBackProjectionAlgorithm.cpp199
-rw-r--r--src/FilteredBackProjectionAlgorithm.cpp137
-rw-r--r--src/Filters.cpp608
4 files changed, 742 insertions, 206 deletions
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 944f165..88e235b 100644
--- a/src/CudaFilteredBackProjectionAlgorithm.cpp
+++ b/src/CudaFilteredBackProjectionAlgorithm.cpp
@@ -27,10 +27,10 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#include <astra/CudaFilteredBackProjectionAlgorithm.h>
#include <astra/FanFlatProjectionGeometry2D.h>
-#include <cstring>
#include "astra/AstraObjectManager.h"
#include "astra/CudaProjector2D.h"
+#include "astra/Filters.h"
#include "astra/cuda/2d/astra.h"
#include "astra/cuda/2d/fbp.h"
@@ -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,9 +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.");
- if((m_eFilter == FILTER_PROJECTION) || (m_eFilter == FILTER_SINOGRAM) || (m_eFilter == FILTER_RPROJECTION) || (m_eFilter == FILTER_RSINOGRAM))
+ ASTRA_CONFIG_CHECK(m_filterConfig.m_eType != FILTER_ERROR, "FBP_CUDA", "Invalid filter name.");
+
+ 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
@@ -229,122 +176,12 @@ 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;
return true;
}
-static bool stringCompareLowerCase(const char * _stringA, const char * _stringB)
-{
- int iCmpReturn = 0;
-
-#ifdef _MSC_VER
- iCmpReturn = _stricmp(_stringA, _stringB);
-#else
- iCmpReturn = strcasecmp(_stringA, _stringB);
-#endif
-
- return (iCmpReturn == 0);
-}
-
-E_FBPFILTER CCudaFilteredBackProjectionAlgorithm::_convertStringToFilter(const char * _filterType)
-{
- E_FBPFILTER output = FILTER_NONE;
-
- if(stringCompareLowerCase(_filterType, "ram-lak"))
- {
- output = FILTER_RAMLAK;
- }
- else if(stringCompareLowerCase(_filterType, "shepp-logan"))
- {
- output = FILTER_SHEPPLOGAN;
- }
- else if(stringCompareLowerCase(_filterType, "cosine"))
- {
- output = FILTER_COSINE;
- }
- else if(stringCompareLowerCase(_filterType, "hamming"))
- {
- output = FILTER_HAMMING;
- }
- else if(stringCompareLowerCase(_filterType, "hann"))
- {
- output = FILTER_HANN;
- }
- else if(stringCompareLowerCase(_filterType, "none"))
- {
- output = FILTER_NONE;
- }
- else if(stringCompareLowerCase(_filterType, "tukey"))
- {
- output = FILTER_TUKEY;
- }
- else if(stringCompareLowerCase(_filterType, "lanczos"))
- {
- output = FILTER_LANCZOS;
- }
- else if(stringCompareLowerCase(_filterType, "triangular"))
- {
- output = FILTER_TRIANGULAR;
- }
- else if(stringCompareLowerCase(_filterType, "gaussian"))
- {
- output = FILTER_GAUSSIAN;
- }
- else if(stringCompareLowerCase(_filterType, "barlett-hann"))
- {
- output = FILTER_BARTLETTHANN;
- }
- else if(stringCompareLowerCase(_filterType, "blackman"))
- {
- output = FILTER_BLACKMAN;
- }
- else if(stringCompareLowerCase(_filterType, "nuttall"))
- {
- output = FILTER_NUTTALL;
- }
- else if(stringCompareLowerCase(_filterType, "blackman-harris"))
- {
- output = FILTER_BLACKMANHARRIS;
- }
- else if(stringCompareLowerCase(_filterType, "blackman-nuttall"))
- {
- output = FILTER_BLACKMANNUTTALL;
- }
- else if(stringCompareLowerCase(_filterType, "flat-top"))
- {
- output = FILTER_FLATTOP;
- }
- else if(stringCompareLowerCase(_filterType, "kaiser"))
- {
- output = FILTER_KAISER;
- }
- else if(stringCompareLowerCase(_filterType, "parzen"))
- {
- output = FILTER_PARZEN;
- }
- else if(stringCompareLowerCase(_filterType, "projection"))
- {
- output = FILTER_PROJECTION;
- }
- else if(stringCompareLowerCase(_filterType, "sinogram"))
- {
- output = FILTER_SINOGRAM;
- }
- else if(stringCompareLowerCase(_filterType, "rprojection"))
- {
- output = FILTER_RPROJECTION;
- }
- else if(stringCompareLowerCase(_filterType, "rsinogram"))
- {
- output = FILTER_RSINOGRAM;
- }
- else
- {
- ASTRA_ERROR("Failed to convert \"%s\" into a filter.",_filterType);
- }
-
- return output;
-}
diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp
index bb2e722..67a12a2 100644
--- a/src/FilteredBackProjectionAlgorithm.cpp
+++ b/src/FilteredBackProjectionAlgorithm.cpp
@@ -81,6 +81,9 @@ void CFilteredBackProjectionAlgorithm::clear()
m_pSinogram = NULL;
m_pReconstruction = NULL;
m_bIsInitialized = false;
+
+ delete[] m_filterConfig.m_pfCustomFilter;
+ m_filterConfig.m_pfCustomFilter = 0;
}
@@ -151,6 +154,9 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
delete[] projectionAngles;
}
+ m_filterConfig = getFilterConfigForAlgorithm(_cfg, this);
+
+
// TODO: check that the angles are linearly spaced between 0 and pi
// success
@@ -195,11 +201,14 @@ bool CFilteredBackProjectionAlgorithm::initialize(CProjector2D* _pProjector,
CFloat32VolumeData2D* _pVolume,
CFloat32ProjectionData2D* _pSinogram)
{
+ clear();
+
// store classes
m_pProjector = _pProjector;
m_pReconstruction = _pVolume;
m_pSinogram = _pSinogram;
+ m_filterConfig = SFilterConfig();
// TODO: check that the angles are linearly spaced between 0 and pi
@@ -214,6 +223,15 @@ bool CFilteredBackProjectionAlgorithm::_check()
{
ASTRA_CONFIG_CHECK(CReconstructionAlgorithm2D::_check(), "FBP", "Error in ReconstructionAlgorithm2D initialization");
+ ASTRA_CONFIG_CHECK(m_filterConfig.m_eType != FILTER_ERROR, "FBP", "Invalid filter name.");
+
+ 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_filterConfig.m_pfCustomFilter, "FBP", "Invalid filter pointer.");
+ }
+
+ ASTRA_CONFIG_CHECK(checkCustomFilterSize(m_filterConfig, *m_pSinogram->getGeometry()), "FBP", "Filter size mismatch");
+
// success
return true;
}
@@ -249,34 +267,84 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D
ASTRA_ASSERT(_pFilteredSinogram->getAngleCount() == m_pSinogram->getAngleCount());
ASTRA_ASSERT(_pFilteredSinogram->getDetectorCount() == m_pSinogram->getDetectorCount());
+ ASTRA_ASSERT(m_filterConfig.m_eType != FILTER_ERROR);
+ if (m_filterConfig.m_eType == FILTER_NONE)
+ return;
int iAngleCount = m_pProjector->getProjectionGeometry()->getProjectionAngleCount();
int iDetectorCount = m_pProjector->getProjectionGeometry()->getDetectorCount();
- // We'll zero-pad to the smallest power of two at least 64 and
- // at least 2*iDetectorCount
- int zpDetector = 64;
- int nextPow2 = 5;
- while (zpDetector < iDetectorCount*2) {
- zpDetector *= 2;
- nextPow2++;
- }
+ int zpDetector = calcNextPowerOfTwo(2 * m_pSinogram->getDetectorCount());
+ int iHalfFFTSize = astra::calcFFTFourierSize(zpDetector);
- // Create filter
- float32* filter = new float32[zpDetector];
+ // cdft setup
+ int *ip = new int[int(2+sqrt((float)zpDetector)+1)];
+ ip[0] = 0;
+ float32 *w = new float32[zpDetector/2];
- for (int iDetector = 0; iDetector <= zpDetector/2; iDetector++)
- filter[iDetector] = (2.0f * iDetector)/zpDetector;
+ // Create filter
+ bool bFilterMultiAngle = false;
+ bool bFilterComplex = false;
+ float *pfFilter = 0;
+ float *pfFilter_delete = 0;
+ switch (m_filterConfig.m_eType) {
+ case FILTER_ERROR:
+ case FILTER_NONE:
+ // Should have been handled before
+ ASTRA_ASSERT(false);
+ return;
+ case FILTER_PROJECTION:
+ // Fourier space, real, half the coefficients (because symmetric)
+ // 1 x iHalfFFTSize
+ pfFilter = m_filterConfig.m_pfCustomFilter;
+ break;
+ case FILTER_SINOGRAM:
+ bFilterMultiAngle = true;
+ pfFilter = m_filterConfig.m_pfCustomFilter;
+ break;
+ case FILTER_RSINOGRAM:
+ bFilterMultiAngle = true;
+ // fall-through
+ case FILTER_RPROJECTION:
+ {
+ bFilterComplex = true;
+
+ int count = bFilterMultiAngle ? iAngleCount : 1;
+ // Spatial, real, full convolution kernel
+ // Center in center (or right-of-center for even sized.)
+ // I.e., 0 1 0 and 0 0 1 0 both correspond to the identity
+
+ pfFilter = new float[2 * zpDetector * count];
+ pfFilter_delete = pfFilter;
+
+ int iUsedFilterWidth = min(m_filterConfig.m_iCustomFilterWidth, zpDetector);
+ int iStartFilterIndex = (m_filterConfig.m_iCustomFilterWidth - iUsedFilterWidth) / 2;
+ int iMaxFilterIndex = iStartFilterIndex + iUsedFilterWidth;
+
+ int iFilterShiftSize = m_filterConfig.m_iCustomFilterWidth / 2;
+
+ for (int i = 0; i < count; ++i) {
+ float *rOut = pfFilter + i * 2 * zpDetector;
+ float *rIn = m_filterConfig.m_pfCustomFilter + i * m_filterConfig.m_iCustomFilterWidth;
+ memset(rOut, 0, sizeof(float) * 2 * zpDetector);
+
+ for(int j = iStartFilterIndex; j < iMaxFilterIndex; j++) {
+ int iFFTInFilterIndex = (j + zpDetector - iFilterShiftSize) % zpDetector;
+ rOut[2 * iFFTInFilterIndex] = rIn[j];
+ }
- for (int iDetector = zpDetector/2+1; iDetector < zpDetector; iDetector++)
- filter[iDetector] = (2.0f * (zpDetector - iDetector)) / zpDetector;
+ cdft(2*zpDetector, -1, rOut, ip, w);
+ }
+ break;
+ }
+ default:
+ pfFilter = genFilter(m_filterConfig, zpDetector, iHalfFFTSize);
+ pfFilter_delete = pfFilter;
+ }
float32* pf = new float32[2 * iAngleCount * zpDetector];
- int *ip = new int[int(2+sqrt((float)zpDetector)+1)];
- ip[0]=0;
- float32 *w = new float32[zpDetector/2];
// Copy and zero-pad data
for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) {
@@ -299,11 +367,34 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D
}
// Filter
- for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) {
- float32* pfRow = pf + iAngle * 2 * zpDetector;
- for (int iDetector = 0; iDetector < zpDetector; ++iDetector) {
- pfRow[2*iDetector] *= filter[iDetector];
- pfRow[2*iDetector+1] *= filter[iDetector];
+ if (bFilterComplex) {
+ for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) {
+ float32* pfRow = pf + iAngle * 2 * zpDetector;
+ float *pfFilterRow = pfFilter;
+ if (bFilterMultiAngle)
+ pfFilterRow += iAngle * 2 * zpDetector;
+
+ for (int i = 0; i < zpDetector; ++i) {
+ float re = pfRow[2*i] * pfFilterRow[2*i] - pfRow[2*i+1] * pfFilterRow[2*i+1];
+ float im = pfRow[2*i] * pfFilterRow[2*i+1] + pfRow[2*i+1] * pfFilterRow[2*i];
+ pfRow[2*i] = re;
+ pfRow[2*i+1] = im;
+ }
+ }
+ } else {
+ for (int iAngle = 0; iAngle < iAngleCount; ++iAngle) {
+ float32* pfRow = pf + iAngle * 2 * zpDetector;
+ float *pfFilterRow = pfFilter;
+ if (bFilterMultiAngle)
+ pfFilterRow += iAngle * iHalfFFTSize;
+ for (int iDetector = 0; iDetector < iHalfFFTSize; ++iDetector) {
+ pfRow[2*iDetector] *= pfFilterRow[iDetector];
+ pfRow[2*iDetector+1] *= pfFilterRow[iDetector];
+ }
+ for (int iDetector = iHalfFFTSize; iDetector < zpDetector; ++iDetector) {
+ pfRow[2*iDetector] *= pfFilterRow[zpDetector - iDetector];
+ pfRow[2*iDetector+1] *= pfFilterRow[zpDetector - iDetector];
+ }
}
}
@@ -324,7 +415,7 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D
delete[] pf;
delete[] w;
delete[] ip;
- delete[] filter;
+ delete[] pfFilter_delete;
}
}
diff --git a/src/Filters.cpp b/src/Filters.cpp
new file mode 100644
index 0000000..c13aa6b
--- /dev/null
+++ b/src/Filters.cpp
@@ -0,0 +1,608 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2018, imec Vision Lab, University of Antwerp
+ 2014-2018, CWI, Amsterdam
+
+Contact: astra@astra-toolbox.com
+Website: http://www.astra-toolbox.com/
+
+This file is part of the ASTRA Toolbox.
+
+
+The ASTRA Toolbox is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+The ASTRA Toolbox is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+*/
+
+#include "astra/Globals.h"
+#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(const SFilterConfig &_cfg,
+ int _iFFTRealDetectorCount,
+ int _iFFTFourierDetectorCount)
+{
+ float * pfFilt = new float[_iFFTFourierDetectorCount];
+ float * pfW = new float[_iFFTFourierDetectorCount];
+
+ // We cache one Fourier transform for repeated FBP's of the same size
+ static float *pfData = 0;
+ static int iFilterCacheSize = 0;
+
+ if (!pfData || iFilterCacheSize != _iFFTRealDetectorCount) {
+ // Compute filter in spatial domain
+
+ delete[] pfData;
+ pfData = new float[2*_iFFTRealDetectorCount];
+ int *ip = new int[int(2+sqrt(_iFFTRealDetectorCount)+1)];
+ ip[0] = 0;
+ float32 *w = new float32[_iFFTRealDetectorCount/2];
+
+ for (int i = 0; i < _iFFTRealDetectorCount; ++i) {
+ pfData[2*i+1] = 0.0f;
+
+ if (i & 1) {
+ int j = i;
+ if (2*j > _iFFTRealDetectorCount)
+ j = _iFFTRealDetectorCount - j;
+ float f = M_PI * j;
+ pfData[2*i] = -1 / (f*f);
+ } else {
+ pfData[2*i] = 0.0f;
+ }
+ }
+
+ pfData[0] = 0.25f;
+
+ cdft(2*_iFFTRealDetectorCount, -1, pfData, ip, w);
+ delete[] ip;
+ delete[] w;
+
+ iFilterCacheSize = _iFFTRealDetectorCount;
+ }
+
+ for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount;
+
+ pfFilt[iDetectorIndex] = 2.0f * pfData[2*iDetectorIndex];
+ pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex;
+ }
+
+ switch(_cfg.m_eType)
+ {
+ case FILTER_RAMLAK:
+ {
+ // do nothing
+ break;
+ }
+ case FILTER_SHEPPLOGAN:
+ {
+ // 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 / _cfg.m_fD) / (pfW[iDetectorIndex] / 2.0f / _cfg.m_fD));
+ }
+ break;
+ }
+ case FILTER_COSINE:
+ {
+ // 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 / _cfg.m_fD);
+ }
+ break;
+ }
+ case FILTER_HAMMING:
+ {
+ // 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] / _cfg.m_fD));
+ }
+ break;
+ }
+ case FILTER_HANN:
+ {
+ // 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] / _cfg.m_fD)) / 2.0f;
+ }
+ break;
+ }
+ case FILTER_TUKEY:
+ {
+ 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;
+ float fDenom = (1.0f - fAlpha) * fHalfN;
+ float fBlockStart = fHalfN - fEnumTerm;
+ float fBlockEnd = fHalfN + fEnumTerm;
+
+ for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fAbsSmallN = fabs((float)iDetectorIndex);
+ float fStoredValue = 0.0f;
+
+ if((fBlockStart <= fAbsSmallN) && (fAbsSmallN <= fBlockEnd))
+ {
+ fStoredValue = 1.0f;
+ }
+ else
+ {
+ float fEnum = fAbsSmallN - fEnumTerm;
+ float fCosInput = M_PI * fEnum / fDenom;
+ fStoredValue = 0.5f * (1.0f + cosf(fCosInput));
+ }
+
+ pfFilt[iDetectorIndex] *= fStoredValue;
+ }
+
+ break;
+ }
+ case FILTER_LANCZOS:
+ {
+ float fDenum = (float)(_iFFTFourierDetectorCount - 1);
+
+ for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fSmallN = (float)iDetectorIndex;
+ float fX = 2.0f * fSmallN / fDenum - 1.0f;
+ float fSinInput = M_PI * fX;
+ float fStoredValue = 0.0f;
+
+ if(fabsf(fSinInput) > 0.001f)
+ {
+ fStoredValue = sin(fSinInput)/fSinInput;
+ }
+ else
+ {
+ fStoredValue = 1.0f;
+ }
+
+ pfFilt[iDetectorIndex] *= fStoredValue;
+ }
+
+ break;
+ }
+ case FILTER_TRIANGULAR:
+ {
+ float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1);
+
+ for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fSmallN = (float)iDetectorIndex;
+ float fAbsInput = fSmallN - fNMinusOne / 2.0f;
+ float fParenInput = fNMinusOne / 2.0f - fabsf(fAbsInput);
+ float fStoredValue = 2.0f / fNMinusOne * fParenInput;
+
+ pfFilt[iDetectorIndex] *= fStoredValue;
+ }
+
+ break;
+ }
+ case FILTER_GAUSSIAN:
+ {
+ 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;
+
+ for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fSmallN = (float)iDetectorIndex;
+ float fEnum = fSmallN - fQuotient;
+ float fDenom = fSigma * fQuotient;
+ float fPower = -0.5f * (fEnum / fDenom) * (fEnum / fDenom);
+ float fStoredValue = expf(fPower);
+
+ pfFilt[iDetectorIndex] *= fStoredValue;
+ }
+
+ break;
+ }
+ case FILTER_BARTLETTHANN:
+ {
+ const float fA0 = 0.62f;
+ const float fA1 = 0.48f;
+ const float fA2 = 0.38f;
+ float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f;
+
+ for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fSmallN = (float)iDetectorIndex;
+ float fAbsInput = fSmallN / fNMinusOne - 0.5f;
+ float fFirstTerm = fA1 * fabsf(fAbsInput);
+ float fCosInput = 2.0f * M_PI * fSmallN / fNMinusOne;
+ float fSecondTerm = fA2 * cosf(fCosInput);
+ float fStoredValue = fA0 - fFirstTerm - fSecondTerm;
+
+ pfFilt[iDetectorIndex] *= fStoredValue;
+ }
+
+ break;
+ }
+ case FILTER_BLACKMAN:
+ {
+ 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;
+ float fNMinusOne = (float)(_iFFTFourierDetectorCount - 1);
+
+ for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fSmallN = (float)iDetectorIndex;
+ float fCosInput1 = 2.0f * M_PI * 0.5f * fSmallN / fNMinusOne;
+ float fCosInput2 = 4.0f * M_PI * 0.5f * fSmallN / fNMinusOne;
+ float fStoredValue = fA0 - fA1 * cosf(fCosInput1) + fA2 * cosf(fCosInput2);
+
+ pfFilt[iDetectorIndex] *= fStoredValue;
+ }
+
+ break;
+ }
+ case FILTER_NUTTALL:
+ {
+ const float fA0 = 0.355768f;
+ const float fA1 = 0.487396f;
+ const float fA2 = 0.144232f;
+ const float fA3 = 0.012604f;
+ float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f;
+
+ for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fSmallN = (float)iDetectorIndex;
+ float fBaseCosInput = M_PI * fSmallN / fNMinusOne;
+ float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput);
+ float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput);
+ float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput);
+ float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm;
+
+ pfFilt[iDetectorIndex] *= fStoredValue;
+ }
+
+ break;
+ }
+ case FILTER_BLACKMANHARRIS:
+ {
+ const float fA0 = 0.35875f;
+ const float fA1 = 0.48829f;
+ const float fA2 = 0.14128f;
+ const float fA3 = 0.01168f;
+ float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f;
+
+ for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fSmallN = (float)iDetectorIndex;
+ float fBaseCosInput = M_PI * fSmallN / fNMinusOne;
+ float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput);
+ float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput);
+ float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput);
+ float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm;
+
+ pfFilt[iDetectorIndex] *= fStoredValue;
+ }
+
+ break;
+ }
+ case FILTER_BLACKMANNUTTALL:
+ {
+ const float fA0 = 0.3635819f;
+ const float fA1 = 0.4891775f;
+ const float fA2 = 0.1365995f;
+ const float fA3 = 0.0106411f;
+ float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f;
+
+ for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fSmallN = (float)iDetectorIndex;
+ float fBaseCosInput = M_PI * fSmallN / fNMinusOne;
+ float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput);
+ float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput);
+ float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput);
+ float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm;
+
+ pfFilt[iDetectorIndex] *= fStoredValue;
+ }
+
+ break;
+ }
+ case FILTER_FLATTOP:
+ {
+ const float fA0 = 1.0f;
+ const float fA1 = 1.93f;
+ const float fA2 = 1.29f;
+ const float fA3 = 0.388f;
+ const float fA4 = 0.032f;
+ float fNMinusOne = (float)(_iFFTFourierDetectorCount) - 1.0f;
+
+ for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fSmallN = (float)iDetectorIndex;
+ float fBaseCosInput = M_PI * fSmallN / fNMinusOne;
+ float fFirstTerm = fA1 * cosf(2.0f * fBaseCosInput);
+ float fSecondTerm = fA2 * cosf(4.0f * fBaseCosInput);
+ float fThirdTerm = fA3 * cosf(6.0f * fBaseCosInput);
+ float fFourthTerm = fA4 * cosf(8.0f * fBaseCosInput);
+ float fStoredValue = fA0 - fFirstTerm + fSecondTerm - fThirdTerm + fFourthTerm;
+
+ pfFilt[iDetectorIndex] *= fStoredValue;
+ }
+
+ break;
+ }
+ case FILTER_KAISER:
+ {
+ 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);
+
+ for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fSmallN = (float)iDetectorIndex;
+ float fSquareInput = 2.0f * fSmallN / fNMinusOne - 1;
+ float fSqrtInput = 1.0f - fSquareInput * fSquareInput;
+ float fBesselInput = fPiTimesAlpha * sqrt(fSqrtInput);
+ float fEnum = (float)j0((double)fBesselInput);
+ float fStoredValue = fEnum / fDenom;
+
+ pfFilt[iDetectorIndex] *= fStoredValue;
+ }
+
+ break;
+ }
+ case FILTER_PARZEN:
+ {
+ for(int iDetectorIndex = 1; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fSmallN = (float)iDetectorIndex;
+ float fQ = fSmallN / (float)(_iFFTFourierDetectorCount - 1);
+ float fStoredValue = 0.0f;
+
+ if(fQ <= 0.5f)
+ {
+ fStoredValue = 1.0f - 6.0f * fQ * fQ * (1.0f - fQ);
+ }
+ else
+ {
+ float fCubedValue = 1.0f - fQ;
+ fStoredValue = 2.0f * fCubedValue * fCubedValue * fCubedValue;
+ }
+
+ pfFilt[iDetectorIndex] *= fStoredValue;
+ }
+
+ break;
+ }
+ default:
+ {
+ ASTRA_ERROR("Cannot serve requested filter");
+ }
+ }
+
+ // filt(w>pi*d) = 0;
+ float fPiTimesD = M_PI * _cfg.m_fD;
+ for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
+ {
+ float fWValue = pfW[iDetectorIndex];
+
+ if(fWValue > fPiTimesD)
+ {
+ pfFilt[iDetectorIndex] = 0.0f;
+ }
+ }
+
+ delete[] pfW;
+
+ return pfFilt;
+}
+
+static bool stringCompareLowerCase(const char * _stringA, const char * _stringB)
+{
+ int iCmpReturn = 0;
+
+#ifdef _MSC_VER
+ iCmpReturn = _stricmp(_stringA, _stringB);
+#else
+ iCmpReturn = strcasecmp(_stringA, _stringB);
+#endif
+
+ return (iCmpReturn == 0);
+}
+
+struct FilterNameMapEntry {
+ const char *m_name;
+ E_FBPFILTER m_type;
+};
+
+E_FBPFILTER convertStringToFilter(const char * _filterType)
+{
+
+ static const FilterNameMapEntry map[] = {
+ { "ram-lak", FILTER_RAMLAK },
+ { "shepp-logan", FILTER_SHEPPLOGAN },
+ { "cosine", FILTER_COSINE },
+ { "hamming", FILTER_HAMMING},
+ { "hann", FILTER_HANN},
+ { "tukey", FILTER_TUKEY },
+ { "lanczos", FILTER_LANCZOS},
+ { "triangular", FILTER_TRIANGULAR},
+ { "gaussian", FILTER_GAUSSIAN},
+ { "barlett-hann", FILTER_BARTLETTHANN },
+ { "blackman", FILTER_BLACKMAN},
+ { "nuttall", FILTER_NUTTALL },
+ { "blackman-harris", FILTER_BLACKMANHARRIS },
+ { "blackman-nuttall", FILTER_BLACKMANNUTTALL },
+ { "flat-top", FILTER_FLATTOP },
+ { "kaiser", FILTER_KAISER },
+ { "parzen", FILTER_PARZEN },
+ { "projection", FILTER_PROJECTION },
+ { "sinogram", FILTER_SINOGRAM },
+ { "rprojection", FILTER_RPROJECTION },
+ { "rsinogram", FILTER_RSINOGRAM },
+ { "none", FILTER_NONE },
+ { 0, FILTER_ERROR } };
+
+ const FilterNameMapEntry *i;
+
+ for (i = &map[0]; i->m_name; ++i)
+ if (stringCompareLowerCase(_filterType, i->m_name))
+ return i->m_type;
+
+ ASTRA_ERROR("Failed to convert \"%s\" into a filter.",_filterType);
+
+ return FILTER_ERROR;
+}
+
+
+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();
+ c.m_iCustomFilterHeight = pFilterData->getGeometry()->getProjectionAngleCount();
+
+ 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!
+
+ // 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;
+}
+
+int calcNextPowerOfTwo(int n)
+{
+ int x = 1;
+ while (x < n && x > 0)
+ 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-09-21 02:59:09 +0000