summaryrefslogtreecommitdiffstats
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/ArtAlgorithm.cpp10
-rw-r--r--src/CudaDataOperationAlgorithm.cpp2
-rw-r--r--src/CudaReconstructionAlgorithm2D.cpp22
-rw-r--r--src/CudaSartAlgorithm.cpp17
-rw-r--r--src/CudaSirtAlgorithm.cpp41
-rw-r--r--src/CudaSirtAlgorithm3D.cpp8
-rw-r--r--src/FilteredBackProjectionAlgorithm.cpp14
-rw-r--r--src/Float32VolumeData3DMemory.cpp1
-rw-r--r--src/SartAlgorithm.cpp32
-rw-r--r--src/SirtAlgorithm.cpp11
-rw-r--r--src/XMLNode.cpp4
11 files changed, 109 insertions, 53 deletions
diff --git a/src/ArtAlgorithm.cpp b/src/ArtAlgorithm.cpp
index b59bd93..526c263 100644
--- a/src/ArtAlgorithm.cpp
+++ b/src/ArtAlgorithm.cpp
@@ -156,8 +156,12 @@ bool CArtAlgorithm::initialize(const Config& _cfg)
return false;
}
+ // "Lambda" is replaced by the more descriptive "Relaxation"
m_fLambda = _cfg.self.getOptionNumerical("Lambda", 1.0f);
- CC.markOptionParsed("Lambda");
+ m_fLambda = _cfg.self.getOptionNumerical("Relaxation", m_fLambda);
+ if (!_cfg.self.hasOption("Relaxation"))
+ CC.markOptionParsed("Lambda");
+ CC.markOptionParsed("Relaxation");
// success
m_bIsInitialized = _check();
@@ -232,7 +236,7 @@ map<string,boost::any> CArtAlgorithm::getInformation()
{
map<string, boost::any> res;
res["RayOrder"] = getInformation("RayOrder");
- res["Lambda"] = getInformation("Lambda");
+ res["Relaxation"] = getInformation("Relaxation");
return mergeMap<string,boost::any>(CReconstructionAlgorithm2D::getInformation(), res);
};
@@ -240,7 +244,7 @@ map<string,boost::any> CArtAlgorithm::getInformation()
// Information - Specific
boost::any CArtAlgorithm::getInformation(std::string _sIdentifier)
{
- if (_sIdentifier == "Lambda") { return m_fLambda; }
+ if (_sIdentifier == "Relaxation") { return m_fLambda; }
if (_sIdentifier == "RayOrder") {
vector<float32> res;
for (int i = 0; i < m_iRayCount; i++) {
diff --git a/src/CudaDataOperationAlgorithm.cpp b/src/CudaDataOperationAlgorithm.cpp
index 15886a4..82b676b 100644
--- a/src/CudaDataOperationAlgorithm.cpp
+++ b/src/CudaDataOperationAlgorithm.cpp
@@ -76,7 +76,7 @@ bool CCudaDataOperationAlgorithm::initialize(const Config& _cfg)
node = _cfg.self.getSingleNode("DataId");
ASTRA_CONFIG_CHECK(node, "CCudaDataOperationAlgorithm", "No DataId tag specified.");
vector<string> data = node.getContentArray();
- for (vector<string>::iterator it = data.begin(); it != data.end(); it++){
+ for (vector<string>::iterator it = data.begin(); it != data.end(); ++it){
int id = StringUtil::stringToInt(*it);
m_pData.push_back(dynamic_cast<CFloat32Data2D*>(CData2DManager::getSingleton().get(id)));
}
diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp
index 5a1910c..2798434 100644
--- a/src/CudaReconstructionAlgorithm2D.cpp
+++ b/src/CudaReconstructionAlgorithm2D.cpp
@@ -329,6 +329,20 @@ bool CCudaReconstructionAlgorithm2D::setupGeometry()
}
//----------------------------------------------------------------------------------------
+
+void CCudaReconstructionAlgorithm2D::initCUDAAlgorithm()
+{
+ bool ok;
+
+ ok = setupGeometry();
+ ASTRA_ASSERT(ok);
+
+ ok = m_pAlgo->allocateBuffers();
+ ASTRA_ASSERT(ok);
+}
+
+
+//----------------------------------------------------------------------------------------
// Iterate
void CCudaReconstructionAlgorithm2D::run(int _iNrIterations)
{
@@ -339,13 +353,7 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations)
const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry();
if (!m_bAlgoInit) {
-
- ok = setupGeometry();
- ASTRA_ASSERT(ok);
-
- ok = m_pAlgo->allocateBuffers();
- ASTRA_ASSERT(ok);
-
+ initCUDAAlgorithm();
m_bAlgoInit = true;
}
diff --git a/src/CudaSartAlgorithm.cpp b/src/CudaSartAlgorithm.cpp
index d202847..bf97224 100644
--- a/src/CudaSartAlgorithm.cpp
+++ b/src/CudaSartAlgorithm.cpp
@@ -107,7 +107,8 @@ bool CCudaSartAlgorithm::initialize(const Config& _cfg)
CC.markOptionParsed("ProjectionOrderList");
}
-
+ m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f);
+ CC.markOptionParsed("Relaxation");
return true;
}
@@ -123,12 +124,26 @@ bool CCudaSartAlgorithm::initialize(CProjector2D* _pProjector,
if (!m_bIsInitialized)
return false;
+ m_fLambda = 1.0f;
+
m_pAlgo = new astraCUDA::SART();
m_bAlgoInit = false;
return true;
}
+//----------------------------------------------------------------------------------------
+
+void CCudaSartAlgorithm::initCUDAAlgorithm()
+{
+ CCudaReconstructionAlgorithm2D::initCUDAAlgorithm();
+
+ astraCUDA::SART* pSart = dynamic_cast<astraCUDA::SART*>(m_pAlgo);
+
+ pSart->setRelaxation(m_fLambda);
+}
+
+
} // namespace astra
diff --git a/src/CudaSirtAlgorithm.cpp b/src/CudaSirtAlgorithm.cpp
index 33e381a..c8dc677 100644
--- a/src/CudaSirtAlgorithm.cpp
+++ b/src/CudaSirtAlgorithm.cpp
@@ -50,6 +50,8 @@ CCudaSirtAlgorithm::CCudaSirtAlgorithm()
m_pMinMask = 0;
m_pMaxMask = 0;
+
+ m_fLambda = 1.0f;
}
//----------------------------------------------------------------------------------------
@@ -86,6 +88,8 @@ bool CCudaSirtAlgorithm::initialize(const Config& _cfg)
}
CC.markOptionParsed("MaxMaskId");
+ m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f);
+ CC.markOptionParsed("Relaxation");
m_pAlgo = new astraCUDA::SIRT();
m_bAlgoInit = false;
@@ -108,41 +112,30 @@ bool CCudaSirtAlgorithm::initialize(CProjector2D* _pProjector,
m_pAlgo = new astraCUDA::SIRT();
m_bAlgoInit = false;
+ m_fLambda = 1.0f;
return true;
}
//----------------------------------------------------------------------------------------
-// Iterate
-void CCudaSirtAlgorithm::run(int _iNrIterations)
+
+void CCudaSirtAlgorithm::initCUDAAlgorithm()
{
- // check initialized
- ASTRA_ASSERT(m_bIsInitialized);
+ CCudaReconstructionAlgorithm2D::initCUDAAlgorithm();
- if (!m_bAlgoInit) {
- // We only override the initialisation step to copy the min/max masks
+ astraCUDA::SIRT* pSirt = dynamic_cast<astraCUDA::SIRT*>(m_pAlgo);
- bool ok = setupGeometry();
+ if (m_pMinMask || m_pMaxMask) {
+ const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry();
+ const float *pfMinMaskData = 0;
+ const float *pfMaxMaskData = 0;
+ if (m_pMinMask) pfMinMaskData = m_pMinMask->getDataConst();
+ if (m_pMaxMask) pfMaxMaskData = m_pMaxMask->getDataConst();
+ bool ok = pSirt->uploadMinMaxMasks(pfMinMaskData, pfMaxMaskData, volgeom.getGridColCount());
ASTRA_ASSERT(ok);
-
- ok = m_pAlgo->allocateBuffers();
- ASTRA_ASSERT(ok);
-
- if (m_pMinMask || m_pMaxMask) {
- const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry();
- astraCUDA::SIRT* pSirt = dynamic_cast<astraCUDA::SIRT*>(m_pAlgo);
- const float *pfMinMaskData = 0;
- const float *pfMaxMaskData = 0;
- if (m_pMinMask) pfMinMaskData = m_pMinMask->getDataConst();
- if (m_pMaxMask) pfMaxMaskData = m_pMaxMask->getDataConst();
- ok = pSirt->uploadMinMaxMasks(pfMinMaskData, pfMaxMaskData, volgeom.getGridColCount());
- ASTRA_ASSERT(ok);
- }
-
- m_bAlgoInit = true;
}
- CCudaReconstructionAlgorithm2D::run(_iNrIterations);
+ pSirt->setRelaxation(m_fLambda);
}
diff --git a/src/CudaSirtAlgorithm3D.cpp b/src/CudaSirtAlgorithm3D.cpp
index 605c470..c819f8e 100644
--- a/src/CudaSirtAlgorithm3D.cpp
+++ b/src/CudaSirtAlgorithm3D.cpp
@@ -56,6 +56,7 @@ CCudaSirtAlgorithm3D::CCudaSirtAlgorithm3D()
m_iGPUIndex = -1;
m_iVoxelSuperSampling = 1;
m_iDetectorSuperSampling = 1;
+ m_fLambda = 1.0f;
}
//----------------------------------------------------------------------------------------
@@ -128,6 +129,8 @@ bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg)
return false;
}
+ m_fLambda = _cfg.self.getOptionNumerical("Relaxation");
+
initializeFromProjector();
// Deprecated options
@@ -135,6 +138,7 @@ bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg)
m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling);
m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex);
m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
+
CC.markOptionParsed("VoxelSuperSampling");
CC.markOptionParsed("DetectorSuperSampling");
CC.markOptionParsed("GPUIndex");
@@ -164,6 +168,8 @@ bool CCudaSirtAlgorithm3D::initialize(CProjector3D* _pProjector,
clear();
}
+ m_fLambda = 1.0f;
+
// required classes
m_pProjector = _pProjector;
m_pSinogram = _pSinogram;
@@ -224,6 +230,8 @@ void CCudaSirtAlgorithm3D::run(int _iNrIterations)
ASTRA_ASSERT(ok);
+ m_pSirt->setRelaxation(m_fLambda);
+
m_bAstraSIRTInit = true;
}
diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp
index c195578..ccbfec6 100644
--- a/src/FilteredBackProjectionAlgorithm.cpp
+++ b/src/FilteredBackProjectionAlgorithm.cpp
@@ -117,12 +117,10 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
int angleCount = projectionIndex.size();
int detectorCount = m_pProjector->getProjectionGeometry()->getDetectorCount();
+ // TODO: There is no need to allocate this. Better just
+ // create the CFloat32ProjectionData2D object directly, and use its
+ // memory.
float32 * sinogramData2D = new float32[angleCount* detectorCount];
- float32 ** sinogramData = new float32*[angleCount];
- for (int i = 0; i < angleCount; i++)
- {
- sinogramData[i] = &(sinogramData2D[i * detectorCount]);
- }
float32 * projectionAngles = new float32[angleCount];
float32 detectorWidth = m_pProjector->getProjectionGeometry()->getDetectorWidth();
@@ -130,6 +128,8 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
for (int i = 0; i < angleCount; i ++) {
if (projectionIndex[i] > m_pProjector->getProjectionGeometry()->getProjectionAngleCount() -1 )
{
+ delete[] sinogramData2D;
+ delete[] projectionAngles;
ASTRA_ERROR("Invalid Projection Index");
return false;
} else {
@@ -139,7 +139,6 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
{
sinogramData2D[i*detectorCount+ iDetector] = m_pSinogram->getData2D()[orgIndex][iDetector];
}
-// sinogramData[i] = m_pSinogram->getSingleProjectionData(projectionIndex[i]);
projectionAngles[i] = m_pProjector->getProjectionGeometry()->getProjectionAngle((int)projectionIndex[i] );
}
@@ -148,6 +147,9 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
CParallelProjectionGeometry2D * pg = new CParallelProjectionGeometry2D(angleCount, detectorCount,detectorWidth,projectionAngles);
m_pProjector = new CParallelBeamLineKernelProjector2D(pg,m_pReconstruction->getGeometry());
m_pSinogram = new CFloat32ProjectionData2D(pg, sinogramData2D);
+
+ delete[] sinogramData2D;
+ delete[] projectionAngles;
}
// TODO: check that the angles are linearly spaced between 0 and pi
diff --git a/src/Float32VolumeData3DMemory.cpp b/src/Float32VolumeData3DMemory.cpp
index af45cb9..14adb1a 100644
--- a/src/Float32VolumeData3DMemory.cpp
+++ b/src/Float32VolumeData3DMemory.cpp
@@ -136,7 +136,6 @@ CFloat32VolumeData2D * CFloat32VolumeData3DMemory::fetchSliceZ(int _iSliceIndex)
CFloat32VolumeData2D* res = new CFloat32VolumeData2D(&volGeom);
// copy data
- int iSliceCount = m_pGeometry->getGridSliceCount();
float * pfTargetData = res->getData();
for(int iRowIndex = 0; iRowIndex < iRowCount; iRowIndex++)
{
diff --git a/src/SartAlgorithm.cpp b/src/SartAlgorithm.cpp
index 9346160..8df3342 100644
--- a/src/SartAlgorithm.cpp
+++ b/src/SartAlgorithm.cpp
@@ -151,6 +151,9 @@ bool CSartAlgorithm::initialize(const Config& _cfg)
CC.markOptionParsed("ProjectionOrderList");
}
+ m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f);
+ CC.markOptionParsed("Relaxation");
+
// create data objects
m_pTotalRayLength = new CFloat32ProjectionData2D(m_pProjector->getProjectionGeometry());
m_pTotalPixelWeight = new CFloat32VolumeData2D(m_pProjector->getVolumeGeometry());
@@ -246,6 +249,7 @@ map<string,boost::any> CSartAlgorithm::getInformation()
{
map<string, boost::any> res;
res["ProjectionOrder"] = getInformation("ProjectionOrder");
+ res["Relaxation"] = getInformation("Relaxation");
return mergeMap<string,boost::any>(CReconstructionAlgorithm2D::getInformation(), res);
};
@@ -253,6 +257,8 @@ map<string,boost::any> CSartAlgorithm::getInformation()
// Information - Specific
boost::any CSartAlgorithm::getInformation(std::string _sIdentifier)
{
+ if (_sIdentifier == "Relaxation")
+ return m_fLambda;
if (_sIdentifier == "ProjectionOrder") {
vector<float32> res;
for (int i = 0; i < m_iProjectionCount; i++) {
@@ -272,9 +278,8 @@ void CSartAlgorithm::run(int _iNrIterations)
m_bShouldAbort = false;
- int iIteration = 0;
-
// data projectors
+ CDataProjectorInterface* pFirstForwardProjector;
CDataProjectorInterface* pForwardProjector;
CDataProjectorInterface* pBackProjector;
@@ -286,13 +291,13 @@ void CSartAlgorithm::run(int _iNrIterations)
m_pProjector,
SinogramMaskPolicy(m_pSinogramMask), // sinogram mask
ReconstructionMaskPolicy(m_pReconstructionMask), // reconstruction mask
- SIRTBPPolicy(m_pReconstruction, m_pDiffSinogram, m_pTotalPixelWeight, m_pTotalRayLength), // SIRT backprojection
+ SIRTBPPolicy(m_pReconstruction, m_pDiffSinogram, m_pTotalPixelWeight, m_pTotalRayLength, m_fLambda), // SIRT backprojection
m_bUseSinogramMask, m_bUseReconstructionMask, true // options on/off
);
// first time forward projection data projector,
// also computes total pixel weight and total ray length
- pForwardProjector = dispatchDataProjector(
+ pFirstForwardProjector = dispatchDataProjector(
m_pProjector,
SinogramMaskPolicy(m_pSinogramMask), // sinogram mask
ReconstructionMaskPolicy(m_pReconstructionMask), // reconstruction mask
@@ -303,16 +308,30 @@ void CSartAlgorithm::run(int _iNrIterations)
m_bUseSinogramMask, m_bUseReconstructionMask, true // options on/off
);
+ // forward projection data projector
+ pForwardProjector = dispatchDataProjector(
+ m_pProjector,
+ SinogramMaskPolicy(m_pSinogramMask), // sinogram mask
+ ReconstructionMaskPolicy(m_pReconstructionMask), // reconstruction mask
+ CombinePolicy<DiffFPPolicy, TotalPixelWeightPolicy>( // 2 basic operations
+ DiffFPPolicy(m_pReconstruction, m_pDiffSinogram, m_pSinogram), // forward projection with difference calculation
+ TotalPixelWeightPolicy(m_pTotalPixelWeight)), // calculate the total pixel weights
+ m_bUseSinogramMask, m_bUseReconstructionMask, true // options on/off
+ );
+
// iteration loop
- for (; iIteration < _iNrIterations && !m_bShouldAbort; ++iIteration) {
+ for (int iIteration = 0; iIteration < _iNrIterations && !m_bShouldAbort; ++iIteration) {
int iProjection = m_piProjectionOrder[m_iIterationCount % m_iProjectionCount];
// forward projection and difference calculation
m_pTotalPixelWeight->setData(0.0f);
- pForwardProjector->projectSingleProjection(iProjection);
+ if (iIteration < m_iProjectionCount)
+ pFirstForwardProjector->projectSingleProjection(iProjection);
+ else
+ pForwardProjector->projectSingleProjection(iProjection);
// backprojection
pBackProjector->projectSingleProjection(iProjection);
// update iteration count
@@ -325,6 +344,7 @@ void CSartAlgorithm::run(int _iNrIterations)
}
+ ASTRA_DELETE(pFirstForwardProjector);
ASTRA_DELETE(pForwardProjector);
ASTRA_DELETE(pBackProjector);
diff --git a/src/SirtAlgorithm.cpp b/src/SirtAlgorithm.cpp
index d9f3a65..ff25648 100644
--- a/src/SirtAlgorithm.cpp
+++ b/src/SirtAlgorithm.cpp
@@ -76,6 +76,7 @@ void CSirtAlgorithm::_clear()
m_pDiffSinogram = NULL;
m_pTmpVolume = NULL;
+ m_fLambda = 1.0f;
m_iIterationCount = 0;
}
@@ -91,6 +92,7 @@ void CSirtAlgorithm::clear()
ASTRA_DELETE(m_pDiffSinogram);
ASTRA_DELETE(m_pTmpVolume);
+ m_fLambda = 1.0f;
m_iIterationCount = 0;
}
@@ -128,6 +130,9 @@ bool CSirtAlgorithm::initialize(const Config& _cfg)
return false;
}
+ m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f);
+ CC.markOptionParsed("Relaxation");
+
// init data objects and data projectors
_init();
@@ -152,6 +157,8 @@ bool CSirtAlgorithm::initialize(CProjector2D* _pProjector,
m_pSinogram = _pSinogram;
m_pReconstruction = _pReconstruction;
+ m_fLambda = 1.0f;
+
// init data objects and data projectors
_init();
@@ -248,7 +255,7 @@ void CSirtAlgorithm::run(int _iNrIterations)
x = 1.0f / x;
else
x = 0.0f;
- pfT[i] = x;
+ pfT[i] = m_fLambda * x;
}
pfT = m_pTotalRayLength->getData();
for (int i = 0; i < m_pTotalRayLength->getSize(); ++i) {
@@ -296,7 +303,7 @@ void CSirtAlgorithm::run(int _iNrIterations)
m_pTmpVolume->setData(0.0f);
pBackProjector->project();
- // divide by pixel weights
+ // multiply with relaxation factor divided by pixel weights
(*m_pTmpVolume) *= (*m_pTotalPixelWeight);
(*m_pReconstruction) += (*m_pTmpVolume);
diff --git a/src/XMLNode.cpp b/src/XMLNode.cpp
index 40a9b22..cf268c2 100644
--- a/src/XMLNode.cpp
+++ b/src/XMLNode.cpp
@@ -158,7 +158,7 @@ vector<string> XMLNode::getContentArray() const
vector<string> res(iSize);
// loop all list item nodes
list<XMLNode> nodes = getNodes("ListItem");
- for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); it++) {
+ for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); ++it) {
int iIndex = it->getAttributeNumerical("index");
string sValue = it->getAttribute("value");
ASTRA_ASSERT(iIndex < iSize);
@@ -290,7 +290,7 @@ vector<float32> XMLNode::getOptionNumericalArray(string _sKey) const
if (!hasOption(_sKey)) return vector<float32>();
list<XMLNode> nodes = getNodes("Option");
- for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); it++) {
+ for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); ++it) {
if (it->getAttribute("key") == _sKey) {
vector<float32> vals = it->getContentNumericalArray();
return vals;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 08:41:09 +0000