diff options
Diffstat (limited to 'src/CudaReconstructionAlgorithm2D.cpp')
| -rw-r--r-- | src/CudaReconstructionAlgorithm2D.cpp | 518 |
1 files changed, 518 insertions, 0 deletions
diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp new file mode 100644 index 0000000..d567158 --- /dev/null +++ b/src/CudaReconstructionAlgorithm2D.cpp @@ -0,0 +1,518 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("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/>. + +----------------------------------------------------------------------- +$Id$ +*/ + +#ifdef ASTRA_CUDA + +#include "astra/CudaReconstructionAlgorithm2D.h" + +#include <boost/lexical_cast.hpp> + +#include "astra/AstraObjectManager.h" +#include "astra/FanFlatProjectionGeometry2D.h" +#include "astra/FanFlatVecProjectionGeometry2D.h" +#include "astra/CudaProjector2D.h" + +#include "../cuda/2d/algo.h" + +#include <ctime> + +using namespace std; + +namespace astra { + +//---------------------------------------------------------------------------------------- +// Constructor +CCudaReconstructionAlgorithm2D::CCudaReconstructionAlgorithm2D() +{ + _clear(); +} + + + +//---------------------------------------------------------------------------------------- +// Destructor +CCudaReconstructionAlgorithm2D::~CCudaReconstructionAlgorithm2D() +{ + delete m_pAlgo; + m_pAlgo = 0; + m_bAlgoInit = false; +} + +void CCudaReconstructionAlgorithm2D::clear() +{ + delete m_pAlgo; + _clear(); +} + +void CCudaReconstructionAlgorithm2D::_clear() +{ + m_bIsInitialized = false; + m_pAlgo = 0; + m_bAlgoInit = false; + CReconstructionAlgorithm2D::_clear(); + + m_iGPUIndex = 0; + m_iDetectorSuperSampling = 1; + m_iPixelSuperSampling = 1; +} + +//--------------------------------------------------------------------------------------- +// Initialize - Config +bool CCudaReconstructionAlgorithm2D::initialize(const Config& _cfg) +{ + ASTRA_ASSERT(_cfg.self); + ConfigStackCheck<CAlgorithm> CC("CudaReconstructionAlgorithm2D", this, _cfg); + + // if already initialized, clear first + if (m_bIsInitialized) { + clear(); + } + + // sinogram data + XMLNode* node = _cfg.self->getSingleNode("ProjectionDataId"); + ASTRA_CONFIG_CHECK(node, "CudaSirt2", "No ProjectionDataId tag specified."); + int id = boost::lexical_cast<int>(node->getContent()); + m_pSinogram = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id)); + ASTRA_DELETE(node); + CC.markNodeParsed("ProjectionDataId"); + + // reconstruction data + node = _cfg.self->getSingleNode("ReconstructionDataId"); + ASTRA_CONFIG_CHECK(node, "CudaSirt2", "No ReconstructionDataId tag specified."); + id = boost::lexical_cast<int>(node->getContent()); + m_pReconstruction = dynamic_cast<CFloat32VolumeData2D*>(CData2DManager::getSingleton().get(id)); + ASTRA_DELETE(node); + CC.markNodeParsed("ReconstructionDataId"); + + // fixed mask + if (_cfg.self->hasOption("ReconstructionMaskId")) { + m_bUseReconstructionMask = true; + id = boost::lexical_cast<int>(_cfg.self->getOption("ReconstructionMaskId")); + m_pReconstructionMask = dynamic_cast<CFloat32VolumeData2D*>(CData2DManager::getSingleton().get(id)); + } + CC.markOptionParsed("ReconstructionMaskId"); + // fixed mask + if (_cfg.self->hasOption("SinogramMaskId")) { + m_bUseSinogramMask = true; + id = boost::lexical_cast<int>(_cfg.self->getOption("SinogramMaskId")); + m_pSinogramMask = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id)); + } + CC.markOptionParsed("SinogramMaskId"); + + // Constraints - NEW + if (_cfg.self->hasOption("MinConstraint")) { + m_bUseMinConstraint = true; + m_fMinValue = _cfg.self->getOptionNumerical("MinConstraint", 0.0f); + CC.markOptionParsed("MinConstraint"); + } else { + // Constraint - OLD + m_bUseMinConstraint = _cfg.self->getOptionBool("UseMinConstraint", false); + CC.markOptionParsed("UseMinConstraint"); + if (m_bUseMinConstraint) { + m_fMinValue = _cfg.self->getOptionNumerical("MinConstraintValue", 0.0f); + CC.markOptionParsed("MinConstraintValue"); + } + } + if (_cfg.self->hasOption("MaxConstraint")) { + m_bUseMaxConstraint = true; + m_fMaxValue = _cfg.self->getOptionNumerical("MaxConstraint", 255.0f); + CC.markOptionParsed("MaxConstraint"); + } else { + // Constraint - OLD + m_bUseMaxConstraint = _cfg.self->getOptionBool("UseMaxConstraint", false); + CC.markOptionParsed("UseMaxConstraint"); + if (m_bUseMaxConstraint) { + m_fMaxValue = _cfg.self->getOptionNumerical("MaxConstraintValue", 0.0f); + CC.markOptionParsed("MaxConstraintValue"); + } + } + + // GPU number + m_iGPUIndex = (int)_cfg.self->getOptionNumerical("GPUindex", 0); + m_iGPUIndex = (int)_cfg.self->getOptionNumerical("GPUIndex", m_iGPUIndex); + CC.markOptionParsed("GPUindex"); + if (!_cfg.self->hasOption("GPUindex")) + CC.markOptionParsed("GPUIndex"); + + // Detector supersampling factor + m_iDetectorSuperSampling = (int)_cfg.self->getOptionNumerical("DetectorSuperSampling", 1); + CC.markOptionParsed("DetectorSuperSampling"); + + // Pixel supersampling factor + m_iPixelSuperSampling = (int)_cfg.self->getOptionNumerical("PixelSuperSampling", 1); + CC.markOptionParsed("PixelSuperSampling"); + + + // This isn't used yet, but passing it is not something to warn about + node = _cfg.self->getSingleNode("ProjectorId"); + if (node) { + id = boost::lexical_cast<int>(node->getContent()); + CProjector2D *projector = CProjector2DManager::getSingleton().get(id); + if (!dynamic_cast<CCudaProjector2D*>(projector)) { + cout << "Warning: non-CUDA Projector2D passed to FP_CUDA" << std::endl; + } + delete node; + } + CC.markNodeParsed("ProjectorId"); + + + return _check(); +} + +//--------------------------------------------------------------------------------------- +// Initialize - C++ +bool CCudaReconstructionAlgorithm2D::initialize(CProjector2D* _pProjector, + CFloat32ProjectionData2D* _pSinogram, + CFloat32VolumeData2D* _pReconstruction) +{ + return initialize(_pProjector, _pSinogram, _pReconstruction, 0, 1); +} + +//--------------------------------------------------------------------------------------- +// Initialize - C++ +bool CCudaReconstructionAlgorithm2D::initialize(CProjector2D* _pProjector, + CFloat32ProjectionData2D* _pSinogram, + CFloat32VolumeData2D* _pReconstruction, + int _iGPUindex, + int _iDetectorSuperSampling, + int _iPixelSuperSampling) +{ + // if already initialized, clear first + if (m_bIsInitialized) { + clear(); + } + + m_pProjector = 0; + + // required classes + m_pSinogram = _pSinogram; + m_pReconstruction = _pReconstruction; + + m_iDetectorSuperSampling = _iDetectorSuperSampling; + m_iPixelSuperSampling = _iPixelSuperSampling; + m_iGPUIndex = _iGPUindex; + + return _check(); +} + + +//---------------------------------------------------------------------------------------- +// Check +bool CCudaReconstructionAlgorithm2D::_check() +{ + // TODO: CLEAN UP + + + // check pointers + //ASTRA_CONFIG_CHECK(m_pProjector, "Reconstruction2D", "Invalid Projector Object."); + ASTRA_CONFIG_CHECK(m_pSinogram, "SIRT_CUDA", "Invalid Projection Data Object."); + ASTRA_CONFIG_CHECK(m_pReconstruction, "SIRT_CUDA", "Invalid Reconstruction Data Object."); + + // check initializations + //ASTRA_CONFIG_CHECK(m_pProjector->isInitialized(), "Reconstruction2D", "Projector Object Not Initialized."); + ASTRA_CONFIG_CHECK(m_pSinogram->isInitialized(), "SIRT_CUDA", "Projection Data Object Not Initialized."); + ASTRA_CONFIG_CHECK(m_pReconstruction->isInitialized(), "SIRT_CUDA", "Reconstruction Data Object Not Initialized."); + + ASTRA_CONFIG_CHECK(m_iDetectorSuperSampling >= 1, "SIRT_CUDA", "DetectorSuperSampling must be a positive integer."); + ASTRA_CONFIG_CHECK(m_iPixelSuperSampling >= 1, "SIRT_CUDA", "PixelSuperSampling must be a positive integer."); + ASTRA_CONFIG_CHECK(m_iGPUIndex >= 0, "SIRT_CUDA", "GPUIndex must be a non-negative integer."); + + // check compatibility between projector and data classes +// ASTRA_CONFIG_CHECK(m_pSinogram->getGeometry()->isEqual(m_pProjector->getProjectionGeometry()), "SIRT_CUDA", "Projection Data not compatible with the specified Projector."); +// ASTRA_CONFIG_CHECK(m_pReconstruction->getGeometry()->isEqual(m_pProjector->getVolumeGeometry()), "SIRT_CUDA", "Reconstruction Data not compatible with the specified Projector."); + + // todo: turn some of these back on + +// ASTRA_CONFIG_CHECK(m_pProjectionGeometry, "SIRT_CUDA", "ProjectionGeometry not specified."); +// ASTRA_CONFIG_CHECK(m_pProjectionGeometry->isInitialized(), "SIRT_CUDA", "ProjectionGeometry not initialized."); +// ASTRA_CONFIG_CHECK(m_pReconstructionGeometry, "SIRT_CUDA", "ReconstructionGeometry not specified."); +// ASTRA_CONFIG_CHECK(m_pReconstructionGeometry->isInitialized(), "SIRT_CUDA", "ReconstructionGeometry not initialized."); + + // check dimensions + //ASTRA_CONFIG_CHECK(m_pSinogram->getAngleCount() == m_pProjectionGeometry->getProjectionAngleCount(), "SIRT_CUDA", "Sinogram data object size mismatch."); + //ASTRA_CONFIG_CHECK(m_pSinogram->getDetectorCount() == m_pProjectionGeometry->getDetectorCount(), "SIRT_CUDA", "Sinogram data object size mismatch."); + //ASTRA_CONFIG_CHECK(m_pReconstruction->getWidth() == m_pReconstructionGeometry->getGridColCount(), "SIRT_CUDA", "Reconstruction data object size mismatch."); + //ASTRA_CONFIG_CHECK(m_pReconstruction->getHeight() == m_pReconstructionGeometry->getGridRowCount(), "SIRT_CUDA", "Reconstruction data object size mismatch."); + + // check restrictions + // TODO: check restrictions built into cuda code + + + // success + m_bIsInitialized = true; + return true; +} + +void CCudaReconstructionAlgorithm2D::setGPUIndex(int _iGPUIndex) +{ + m_iGPUIndex = _iGPUIndex; +} + + +//--------------------------------------------------------------------------------------- +// Information - All +map<string,boost::any> CCudaReconstructionAlgorithm2D::getInformation() +{ + // TODO: Verify and clean up + + map<string,boost::any> res; + res["ProjectionGeometry"] = getInformation("ProjectionGeometry"); + res["ReconstructionGeometry"] = getInformation("ReconstructionGeometry"); + res["ProjectionDataId"] = getInformation("ProjectionDataId"); + res["ReconstructionDataId"] = getInformation("ReconstructionDataId"); + res["ReconstructionMaskId"] = getInformation("ReconstructionMaskId"); + res["GPUindex"] = getInformation("GPUindex"); + res["DetectorSuperSampling"] = getInformation("DetectorSuperSampling"); + res["PixelSuperSampling"] = getInformation("PixelSuperSampling"); + res["UseMinConstraint"] = getInformation("UseMinConstraint"); + res["MinConstraintValue"] = getInformation("MinConstraintValue"); + res["UseMaxConstraint"] = getInformation("UseMaxConstraint"); + res["MaxConstraintValue"] = getInformation("MaxConstraintValue"); + return mergeMap<string,boost::any>(CReconstructionAlgorithm2D::getInformation(), res); +} + +//--------------------------------------------------------------------------------------- +// Information - Specific +boost::any CCudaReconstructionAlgorithm2D::getInformation(std::string _sIdentifier) +{ + // TODO: Verify and clean up + + if (_sIdentifier == "UseMinConstraint") { return m_bUseMinConstraint ? string("yes") : string("no"); } + if (_sIdentifier == "MinConstraintValue") { return m_fMinValue; } + if (_sIdentifier == "UseMaxConstraint") { return m_bUseMaxConstraint ? string("yes") : string("no"); } + if (_sIdentifier == "MaxConstraintValue") { return m_fMaxValue; } + + // TODO: store these so we can return them? + if (_sIdentifier == "ProjectionGeometry") { return string("not implemented"); } + if (_sIdentifier == "ReconstructionGeometry") { return string("not implemented"); } + if (_sIdentifier == "GPUindex") { return m_iGPUIndex; } + if (_sIdentifier == "DetectorSuperSampling") { return m_iDetectorSuperSampling; } + if (_sIdentifier == "PixelSuperSampling") { return m_iPixelSuperSampling; } + + if (_sIdentifier == "ProjectionDataId") { + int iIndex = CData2DManager::getSingleton().getIndex(m_pSinogram); + if (iIndex != 0) return iIndex; + return std::string("not in manager"); + } + if (_sIdentifier == "ReconstructionDataId") { + int iIndex = CData2DManager::getSingleton().getIndex(m_pReconstruction); + if (iIndex != 0) return iIndex; + return std::string("not in manager"); + } + if (_sIdentifier == "ReconstructionMaskId") { + if (!m_bUseReconstructionMask) return string("not used"); + int iIndex = CData2DManager::getSingleton().getIndex(m_pReconstructionMask); + if (iIndex != 0) return iIndex; + return std::string("not in manager"); + } + return CReconstructionAlgorithm2D::getInformation(_sIdentifier); +} + +bool CCudaReconstructionAlgorithm2D::setupGeometry() +{ + ASTRA_ASSERT(m_bIsInitialized); + ASTRA_ASSERT(!m_bAlgoInit); + + bool ok; + + // TODO: Probably not the best place for this... + ok = m_pAlgo->setGPUIndex(m_iGPUIndex); + if (!ok) return false; + + astraCUDA::SDimensions dims; + + const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); + + // TODO: off-center geometry, non-square pixels + dims.iVolWidth = volgeom.getGridColCount(); + dims.iVolHeight = volgeom.getGridRowCount(); + float fPixelSize = volgeom.getPixelLengthX(); + + dims.iRaysPerDet = m_iDetectorSuperSampling; + dims.iRaysPerPixelDim = m_iPixelSuperSampling; + + + const CParallelProjectionGeometry2D* parProjGeom = dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry()); + const CFanFlatProjectionGeometry2D* fanProjGeom = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pSinogram->getGeometry()); + const CFanFlatVecProjectionGeometry2D* fanVecProjGeom = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pSinogram->getGeometry()); + + if (parProjGeom) { + + dims.iProjAngles = parProjGeom->getProjectionAngleCount(); + dims.iProjDets = parProjGeom->getDetectorCount(); + dims.fDetScale = parProjGeom->getDetectorWidth() / fPixelSize; + + ok = m_pAlgo->setGeometry(dims, parProjGeom->getProjectionAngles()); + + } else if (fanProjGeom) { + + dims.iProjAngles = fanProjGeom->getProjectionAngleCount(); + dims.iProjDets = fanProjGeom->getDetectorCount(); + dims.fDetScale = fanProjGeom->getDetectorWidth() / fPixelSize; + float fOriginSourceDistance = fanProjGeom->getOriginSourceDistance(); + float fOriginDetectorDistance = fanProjGeom->getOriginDetectorDistance(); + + const float* angles = fanProjGeom->getProjectionAngles(); + + astraCUDA::SFanProjection* projs; + projs = new astraCUDA::SFanProjection[dims.iProjAngles]; + + float fSrcX0 = 0.0f; + float fSrcY0 = -fOriginSourceDistance / fPixelSize; + float fDetUX0 = dims.fDetScale; + float fDetUY0 = 0.0f; + float fDetSX0 = dims.iProjDets * fDetUX0 / -2.0f; + float fDetSY0 = fOriginDetectorDistance / fPixelSize; + +#define ROTATE0(name,i,alpha) do { projs[i].f##name##X = f##name##X0 * cos(alpha) - f##name##Y0 * sin(alpha); projs[i].f##name##Y = f##name##X0 * sin(alpha) + f##name##Y0 * cos(alpha); } while(0) + for (unsigned int i = 0; i < dims.iProjAngles; ++i) { + ROTATE0(Src, i, angles[i]); + ROTATE0(DetS, i, angles[i]); + ROTATE0(DetU, i, angles[i]); + } + +#undef ROTATE0 + + ok = m_pAlgo->setFanGeometry(dims, projs); + delete[] projs; + + } else if (fanVecProjGeom) { + + dims.iProjAngles = fanVecProjGeom->getProjectionAngleCount(); + dims.iProjDets = fanVecProjGeom->getDetectorCount(); + dims.fDetScale = fanVecProjGeom->getDetectorWidth() / fPixelSize; + + const astraCUDA::SFanProjection* projs; + projs = fanVecProjGeom->getProjectionVectors(); + + // Rescale projs to fPixelSize == 1 + + astraCUDA::SFanProjection* scaledProjs = new astraCUDA::SFanProjection[dims.iProjAngles]; +#define SCALE(name,i,alpha) do { scaledProjs[i].f##name##X = projs[i].f##name##X * alpha; scaledProjs[i].f##name##Y = projs[i].f##name##Y * alpha; } while (0) + for (unsigned int i = 0; i < dims.iProjAngles; ++i) { + SCALE(Src,i,1.0f/fPixelSize); + SCALE(DetS,i,1.0f/fPixelSize); + SCALE(DetU,i,1.0f/fPixelSize); + } + + ok = m_pAlgo->setFanGeometry(dims, scaledProjs); + + delete[] scaledProjs; + + } else { + + ASTRA_ASSERT(false); + + } + if (!ok) return false; + + + if (m_bUseReconstructionMask) + ok &= m_pAlgo->enableVolumeMask(); + if (!ok) return false; + if (m_bUseSinogramMask) + ok &= m_pAlgo->enableSinogramMask(); + if (!ok) return false; + + const float *pfTOffsets = m_pSinogram->getGeometry()->getExtraDetectorOffset(); + if (pfTOffsets) + ok &= m_pAlgo->setTOffsets(pfTOffsets); + if (!ok) return false; + + ok &= m_pAlgo->init(); + if (!ok) return false; + + + return true; +} + +//---------------------------------------------------------------------------------------- +// Iterate +void CCudaReconstructionAlgorithm2D::run(int _iNrIterations) +{ + // check initialized + ASTRA_ASSERT(m_bIsInitialized); + + bool ok = true; + const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); + + if (!m_bAlgoInit) { + + ok = setupGeometry(); + ASTRA_ASSERT(ok); + + ok = m_pAlgo->allocateBuffers(); + ASTRA_ASSERT(ok); + + m_bAlgoInit = true; + } + + float fPixelSize = volgeom.getPixelLengthX(); + float fSinogramScale = 1.0f/(fPixelSize*fPixelSize); + + ok = m_pAlgo->copyDataToGPU(m_pSinogram->getDataConst(), m_pSinogram->getGeometry()->getDetectorCount(), fSinogramScale, + m_pReconstruction->getDataConst(), volgeom.getGridColCount(), + m_bUseReconstructionMask ? m_pReconstructionMask->getDataConst() : 0, volgeom.getGridColCount(), + m_bUseSinogramMask ? m_pSinogramMask->getDataConst() : 0, m_pSinogram->getGeometry()->getDetectorCount()); + + ASTRA_ASSERT(ok); + + if (m_bUseMinConstraint) + ok &= m_pAlgo->setMinConstraint(m_fMinValue); + if (m_bUseMaxConstraint) + ok &= m_pAlgo->setMaxConstraint(m_fMaxValue); + + ok &= m_pAlgo->iterate(_iNrIterations); + ASTRA_ASSERT(ok); + + ok &= m_pAlgo->getReconstruction(m_pReconstruction->getData(), + volgeom.getGridColCount()); + + ASTRA_ASSERT(ok); +} + +void CCudaReconstructionAlgorithm2D::signalAbort() +{ + if (m_bIsInitialized && m_pAlgo) { + m_pAlgo->signalAbort(); + } +} + +bool CCudaReconstructionAlgorithm2D::getResidualNorm(float32& _fNorm) +{ + if (!m_bIsInitialized || !m_pAlgo) + return false; + + _fNorm = m_pAlgo->computeDiffNorm(); + + return true; +} + +} // namespace astra + +#endif // ASTRA_CUDA |