summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--include/astra/FilteredBackProjectionAlgorithm.h5
-rw-r--r--src/FilteredBackProjectionAlgorithm.cpp69
2 files changed, 54 insertions, 20 deletions
diff --git a/include/astra/FilteredBackProjectionAlgorithm.h b/include/astra/FilteredBackProjectionAlgorithm.h
index 1cd4296..a234845 100644
--- a/include/astra/FilteredBackProjectionAlgorithm.h
+++ b/include/astra/FilteredBackProjectionAlgorithm.h
@@ -35,6 +35,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#include "Projector2D.h"
#include "Float32ProjectionData2D.h"
#include "Float32VolumeData2D.h"
+#include "Filters.h"
namespace astra {
@@ -144,6 +145,10 @@ public:
*/
virtual std::string description() const;
+protected:
+
+ SFilterConfig m_filterConfig;
+
};
// inline functions
diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp
index bb2e722..ed72aa6 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,29 +267,36 @@ 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];
-
- for (int iDetector = 0; iDetector <= zpDetector/2; iDetector++)
- filter[iDetector] = (2.0f * iDetector)/zpDetector;
-
- for (int iDetector = zpDetector/2+1; iDetector < zpDetector; iDetector++)
- filter[iDetector] = (2.0f * (zpDetector - iDetector)) / zpDetector;
-
+ float *pfFilter = 0;
+ switch (m_filterConfig.m_eType) {
+ case FILTER_ERROR:
+ case FILTER_NONE:
+ // Should have been handled before
+ ASTRA_ASSERT(false);
+ return;
+ case FILTER_SINOGRAM:
+ case FILTER_PROJECTION:
+ case FILTER_RSINOGRAM:
+ case FILTER_RPROJECTION:
+ // TODO
+ ASTRA_ERROR("Unimplemented filter type");
+ return;
+
+ default:
+ pfFilter = genFilter(m_filterConfig, zpDetector, iHalfFFTSize);
+ }
float32* pf = new float32[2 * iAngleCount * zpDetector];
int *ip = new int[int(2+sqrt((float)zpDetector)+1)];
@@ -301,9 +326,13 @@ 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];
+ for (int iDetector = 0; iDetector < iHalfFFTSize; ++iDetector) {
+ pfRow[2*iDetector] *= pfFilter[iDetector];
+ pfRow[2*iDetector+1] *= pfFilter[iDetector];
+ }
+ for (int iDetector = iHalfFFTSize; iDetector < zpDetector; ++iDetector) {
+ pfRow[2*iDetector] *= pfFilter[zpDetector - iDetector];
+ pfRow[2*iDetector+1] *= pfFilter[zpDetector - iDetector];
}
}
@@ -324,7 +353,7 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D
delete[] pf;
delete[] w;
delete[] ip;
- delete[] filter;
+ delete[] pfFilter;
}
}