diff options
Diffstat (limited to 'src')
25 files changed, 746 insertions, 267 deletions
diff --git a/src/AsyncAlgorithm.cpp b/src/AsyncAlgorithm.cpp index fcc4dcb..b265f59 100644 --- a/src/AsyncAlgorithm.cpp +++ b/src/AsyncAlgorithm.cpp @@ -160,32 +160,6 @@ void CAsyncAlgorithm::runWrapped(int _iNrIterations) m_bDone = true; } -void CAsyncAlgorithm::timedJoin(int _milliseconds) -{ -#ifndef USE_PTHREADS - if (m_pThread) { - boost::posix_time::milliseconds rel(_milliseconds); - bool res = m_pThread->timed_join(rel); - if (res) { - delete m_pThread; - m_pThread = 0; - m_bThreadStarted = false; - } - } -#else - if (m_bThreadStarted) { - struct timespec abstime; - clock_gettime(CLOCK_REALTIME, &abstime); - abstime.tv_sec += _milliseconds / 1000; - abstime.tv_nsec += (_milliseconds % 1000) * 1000000L; - int err = pthread_timedjoin_np(m_thread, 0, &abstime); - if (err == 0) { - m_bThreadStarted = false; - } - } -#endif -} - void CAsyncAlgorithm::signalAbort() { if (m_pAlg) diff --git a/src/ConeProjectionGeometry3D.cpp b/src/ConeProjectionGeometry3D.cpp index ec64701..1976901 100644 --- a/src/ConeProjectionGeometry3D.cpp +++ b/src/ConeProjectionGeometry3D.cpp @@ -28,6 +28,8 @@ $Id$ #include "astra/ConeProjectionGeometry3D.h" +#include "astra/Logging.h" + #include <boost/lexical_cast.hpp> #include <cstring> @@ -186,17 +188,22 @@ bool CConeProjectionGeometry3D::isOfType(const std::string& _sType) const } //---------------------------------------------------------------------------------------- -void CConeProjectionGeometry3D::toXML(XMLNode* _sNode) const +// Get the configuration object +Config* CConeProjectionGeometry3D::getConfiguration() const { - _sNode->addAttribute("type", "cone"); - _sNode->addChildNode("DetectorSpacingX", m_fDetectorSpacingX); - _sNode->addChildNode("DetectorSpacingY", m_fDetectorSpacingY); - _sNode->addChildNode("DetectorRowCount", m_iDetectorRowCount); - _sNode->addChildNode("DetectorColCount", m_iDetectorColCount); - _sNode->addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount); - _sNode->addChildNode("DistanceOriginDetector", m_fOriginDetectorDistance); - _sNode->addChildNode("DistanceOriginSource", m_fOriginSourceDistance); + Config* cfg = new Config(); + cfg->initialize("ProjectionGeometry3D"); + cfg->self->addAttribute("type", "cone"); + cfg->self->addChildNode("DetectorSpacingX", m_fDetectorSpacingX); + cfg->self->addChildNode("DetectorSpacingY", m_fDetectorSpacingY); + cfg->self->addChildNode("DetectorRowCount", m_iDetectorRowCount); + cfg->self->addChildNode("DetectorColCount", m_iDetectorColCount); + cfg->self->addChildNode("DistanceOriginDetector", m_fOriginDetectorDistance); + cfg->self->addChildNode("DistanceOriginSource", m_fOriginSourceDistance); + cfg->self->addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount); + return cfg; } + //---------------------------------------------------------------------------------------- CVector3D CConeProjectionGeometry3D::getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex) const @@ -225,4 +232,36 @@ CVector3D CConeProjectionGeometry3D::getProjectionDirection(int _iProjectionInde return ret; } +void CConeProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, + int iAngleIndex, + float32 &fU, float32 &fV) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + float alpha = m_pfProjectionAngles[iAngleIndex]; + + // Project point onto optical axis + + // Projector direction is (cos(alpha), sin(alpha)) + // Vector source->origin is (-sin(alpha), cos(alpha)) + + // Distance from source, projected on optical axis + float fD = -sin(alpha) * fX + cos(alpha) * fY + m_fOriginSourceDistance; + + // Scale fZ to detector plane + fV = detectorOffsetYToRowIndexFloat( (fZ * (m_fOriginSourceDistance + m_fOriginDetectorDistance)) / fD ); + + + // Orthogonal distance in XY-plane to optical axis + float fS = cos(alpha) * fX + sin(alpha) * fY; + + // Scale fS to detector plane + fU = detectorOffsetXToColIndexFloat( (fS * (m_fOriginSourceDistance + m_fOriginDetectorDistance)) / fD ); + + ASTRA_DEBUG("alpha: %f, D: %f, V: %f, S: %f, U: %f", alpha, fD, fV, fS, fU); + +} + + } // end namespace astra diff --git a/src/ConeVecProjectionGeometry3D.cpp b/src/ConeVecProjectionGeometry3D.cpp index 512b26c..9dc273d 100644 --- a/src/ConeVecProjectionGeometry3D.cpp +++ b/src/ConeVecProjectionGeometry3D.cpp @@ -198,18 +198,42 @@ bool CConeVecProjectionGeometry3D::isEqual(const CProjectionGeometry3D * _pGeom2 // is of type bool CConeVecProjectionGeometry3D::isOfType(const std::string& _sType) const { - return (_sType == "cone3d_vec"); + return (_sType == "cone_vec"); } //---------------------------------------------------------------------------------------- -void CConeVecProjectionGeometry3D::toXML(XMLNode* _sNode) const +// Get the configuration object +Config* CConeVecProjectionGeometry3D::getConfiguration() const { - _sNode->addAttribute("type","cone3d_vec"); - _sNode->addChildNode("DetectorRowCount", m_iDetectorRowCount); - _sNode->addChildNode("DetectorColCount", m_iDetectorColCount); - // TODO: - //_sNode->addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount); + Config* cfg = new Config(); + cfg->initialize("ProjectionGeometry3D"); + + cfg->self->addAttribute("type", "cone_vec"); + cfg->self->addChildNode("DetectorRowCount", m_iDetectorRowCount); + cfg->self->addChildNode("DetectorColCount", m_iDetectorColCount); + + std::string vectors = ""; + for (int i = 0; i < m_iProjectionAngleCount; ++i) { + SConeProjection& p = m_pProjectionAngles[i]; + vectors += boost::lexical_cast<string>(p.fSrcX) + ","; + vectors += boost::lexical_cast<string>(p.fSrcY) + ","; + vectors += boost::lexical_cast<string>(p.fSrcZ) + ","; + vectors += boost::lexical_cast<string>(p.fDetSX + 0.5f*m_iDetectorRowCount*p.fDetVX + 0.5f*m_iDetectorColCount*p.fDetUX) + ","; + vectors += boost::lexical_cast<string>(p.fDetSY + 0.5f*m_iDetectorRowCount*p.fDetVY + 0.5f*m_iDetectorColCount*p.fDetUY) + ","; + vectors += boost::lexical_cast<string>(p.fDetSZ + 0.5f*m_iDetectorRowCount*p.fDetVZ + 0.5f*m_iDetectorColCount*p.fDetUZ) + ","; + vectors += boost::lexical_cast<string>(p.fDetUX) + ","; + vectors += boost::lexical_cast<string>(p.fDetUY) + ","; + vectors += boost::lexical_cast<string>(p.fDetUZ) + ","; + vectors += boost::lexical_cast<string>(p.fDetVX) + ","; + vectors += boost::lexical_cast<string>(p.fDetVY) + ","; + vectors += boost::lexical_cast<string>(p.fDetVZ); + if (i < m_iProjectionAngleCount-1) vectors += ';'; + } + cfg->self->addChildNode("Vectors", vectors); + + return cfg; } +//---------------------------------------------------------------------------------------- CVector3D CConeVecProjectionGeometry3D::getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex) const { @@ -220,6 +244,26 @@ CVector3D CConeVecProjectionGeometry3D::getProjectionDirection(int _iProjectionI return CVector3D(p.fDetSX + (u+0.5)*p.fDetUX + (v+0.5)*p.fDetVX - p.fSrcX, p.fDetSY + (u+0.5)*p.fDetUY + (v+0.5)*p.fDetVY - p.fSrcY, p.fDetSZ + (u+0.5)*p.fDetUZ + (v+0.5)*p.fDetVZ - p.fSrcZ); } +void CConeVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, + int iAngleIndex, + float32 &fU, float32 &fV) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + double fUX, fUY, fUZ, fUC; + double fVX, fVY, fVZ, fVC; + double fDX, fDY, fDZ, fDC; + + computeBP_UV_Coeffs(m_pProjectionAngles[iAngleIndex], + fUX, fUY, fUZ, fUC, fVX, fVY, fVZ, fVC, fDX, fDY, fDZ, fDC); + + // The -0.5f shifts from corner to center of detector pixels + double fD = fDX*fX + fDY*fY + fDZ*fZ + fDC; + fU = (fUX*fX + fUY*fY + fUZ*fZ + fUC) / fD - 0.5f; + fV = (fVX*fX + fVY*fY + fVZ*fZ + fVC) / fD - 0.5f; +} + //---------------------------------------------------------------------------------------- diff --git a/src/Config.cpp b/src/Config.cpp index 3a3cb53..d860638 100644 --- a/src/Config.cpp +++ b/src/Config.cpp @@ -37,6 +37,9 @@ $Id$ #include "astra/Projector2D.h" #include "astra/Projector3D.h" +#include "astra/Logging.h" +#include <sstream> + using namespace astra; using namespace std; @@ -54,12 +57,24 @@ Config::Config(XMLNode* _self) self = _self; } +//----------------------------------------------------------------------------- Config::~Config() { delete self; self = 0; } +//----------------------------------------------------------------------------- +void Config::initialize(std::string rootname) +{ + if (self == 0) { + XMLDocument* doc = XMLDocument::createDocument(rootname); + self = doc->getRootNode(); + } +} + + +//----------------------------------------------------------------------------- template <class T> ConfigStackCheck<T>::ConfigStackCheck(const char *_name, T* _obj, const Config& _cfg) : object(_obj), cfg(&_cfg), name(_name) @@ -132,7 +147,9 @@ bool ConfigStackCheck<T>::stopParsing() nodes.clear(); if (!errors.empty()) { - cout << "Warning: " << name << ": unused configuration options: " << errors << std::endl; + ostringstream os; + os << name << ": unused configuration options: " << errors; + ASTRA_WARN(os.str().c_str()); return false; } diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp index 6593313..0badc20 100644 --- a/src/CudaFilteredBackProjectionAlgorithm.cpp +++ b/src/CudaFilteredBackProjectionAlgorithm.cpp @@ -32,8 +32,9 @@ $Id$ #include <cstring> #include "astra/AstraObjectManager.h" +#include "../cuda/2d/astra.h" -#include <astra/Logger.h> +#include "astra/Logging.h" using namespace std; using namespace astra; @@ -104,7 +105,7 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) } CC.markNodeParsed("FilterType"); ASTRA_DELETE(node); - + // filter node = _cfg.self->getSingleNode("FilterSinogramId"); if(node != NULL) @@ -167,7 +168,7 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) CC.markOptionParsed("ShortScan"); } - + m_pFBP = new AstraFBP; @@ -185,7 +186,7 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D * { clear(); } - + // required classes m_pSinogram = _pSinogram; m_pReconstruction = _pReconstruction; @@ -241,28 +242,65 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */) bool ok = true; - // TODO: off-center geometry, non-square pixels + // TODO: non-square pixels? ok &= m_pFBP->setReconstructionGeometry(volgeom.getGridColCount(), volgeom.getGridRowCount(), volgeom.getPixelLengthX()); - // TODO: off-center geometry int iDetectorCount; if (parprojgeom) { + + float *offsets, *angles, detSize, outputScale; + + ok = convertAstraGeometry(&volgeom, parprojgeom, offsets, angles, detSize, outputScale); + + ok &= m_pFBP->setProjectionGeometry(parprojgeom->getProjectionAngleCount(), parprojgeom->getDetectorCount(), - parprojgeom->getProjectionAngles(), + angles, parprojgeom->getDetectorWidth()); iDetectorCount = parprojgeom->getDetectorCount(); + + // TODO: Are detSize and outputScale handled correctly? + + if (offsets) + ok &= m_pFBP->setTOffsets(offsets); + ASTRA_ASSERT(ok); + + delete[] offsets; + delete[] angles; + } else if (fanprojgeom) { + + astraCUDA::SFanProjection* projs; + float outputScale; + + // FIXME: Implement this, and clean up the interface to AstraFBP. + if (abs(volgeom.getWindowMinX() + volgeom.getWindowMaxX()) > 0.00001 * volgeom.getPixelLengthX()) { + // Off-center volume geometry isn't supported yet + ASTRA_ASSERT(false); + } + if (abs(volgeom.getWindowMinY() + volgeom.getWindowMaxY()) > 0.00001 * volgeom.getPixelLengthY()) { + // Off-center volume geometry isn't supported yet + ASTRA_ASSERT(false); + } + + ok = convertAstraGeometry(&volgeom, fanprojgeom, projs, outputScale); + + // CHECKME: outputScale? + ok &= m_pFBP->setFanGeometry(fanprojgeom->getProjectionAngleCount(), - fanprojgeom->getDetectorCount(), - fanprojgeom->getProjectionAngles(), - fanprojgeom->getOriginSourceDistance(), - fanprojgeom->getOriginDetectorDistance(), - fanprojgeom->getDetectorWidth(), - m_bShortScan); + fanprojgeom->getDetectorCount(), + projs, + fanprojgeom->getProjectionAngles(), + fanprojgeom->getOriginSourceDistance(), + fanprojgeom->getOriginDetectorDistance(), + + fanprojgeom->getDetectorWidth(), + m_bShortScan); iDetectorCount = fanprojgeom->getDetectorCount(); + + delete[] projs; } else { assert(false); } @@ -271,11 +309,6 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */) ASTRA_ASSERT(ok); - const float *pfTOffsets = m_pSinogram->getGeometry()->getExtraDetectorOffset(); - if (pfTOffsets) - ok &= m_pFBP->setTOffsets(pfTOffsets); - ASTRA_ASSERT(ok); - ok &= m_pFBP->init(m_iGPUIndex); ASTRA_ASSERT(ok); @@ -293,7 +326,7 @@ void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */) const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); ok &= m_pFBP->getReconstruction(m_pReconstruction->getData(), volgeom.getGridColCount()); - + ASTRA_ASSERT(ok); } @@ -302,7 +335,7 @@ bool CCudaFilteredBackProjectionAlgorithm::check() // check pointers 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_pfFilter, "FBP_CUDA", "Invalid filter pointer."); @@ -354,7 +387,7 @@ static bool stringCompareLowerCase(const char * _stringA, const char * _stringB) #else iCmpReturn = strcasecmp(_stringA, _stringB); #endif - + return (iCmpReturn == 0); } @@ -452,7 +485,7 @@ E_FBPFILTER CCudaFilteredBackProjectionAlgorithm::_convertStringToFilter(const c } else { - cerr << "Failed to convert \"" << _filterType << "\" into a filter." << endl; + ASTRA_ERROR("Failed to convert \"%s\" into a filter.",_filterType); } return output; diff --git a/src/CudaForwardProjectionAlgorithm.cpp b/src/CudaForwardProjectionAlgorithm.cpp index 4865887..95abb62 100644 --- a/src/CudaForwardProjectionAlgorithm.cpp +++ b/src/CudaForwardProjectionAlgorithm.cpp @@ -42,6 +42,8 @@ $Id$ #include "astra/FanFlatVecProjectionGeometry2D.h" #include "astra/CudaProjector2D.h" +#include "astra/Logging.h" + using namespace std; namespace astra { @@ -104,7 +106,7 @@ bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg) 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; + ASTRA_WARN("non-CUDA Projector2D passed to FP_CUDA"); } delete node; } @@ -214,52 +216,45 @@ void CCudaForwardProjectionAlgorithm::run(int) bool ok = false; if (parProjGeom) { + + float *offsets, *angles, detSize, outputScale; + ok = convertAstraGeometry(pVolGeom, parProjGeom, offsets, angles, detSize, outputScale); + + ASTRA_ASSERT(ok); // FIXME + + // FIXME: Output scaling + ok = astraCudaFP(m_pVolume->getDataConst(), m_pSinogram->getData(), pVolGeom->getGridColCount(), pVolGeom->getGridRowCount(), parProjGeom->getProjectionAngleCount(), parProjGeom->getDetectorCount(), - parProjGeom->getProjectionAngles(), - parProjGeom->getExtraDetectorOffset(), parProjGeom->getDetectorWidth() / pVolGeom->getPixelLengthX(), - m_iDetectorSuperSampling, m_iGPUIndex); + angles, offsets, detSize, + m_iDetectorSuperSampling, 1.0f * outputScale, m_iGPUIndex); - } else if (fanProjGeom) { + delete[] offsets; + delete[] angles; - ok = astraCudaFanFP(m_pVolume->getDataConst(), m_pSinogram->getData(), - pVolGeom->getGridColCount(), pVolGeom->getGridRowCount(), - fanProjGeom->getProjectionAngleCount(), - fanProjGeom->getDetectorCount(), - fanProjGeom->getProjectionAngles(), - fanProjGeom->getOriginSourceDistance(), - fanProjGeom->getOriginDetectorDistance(), - pVolGeom->getPixelLengthX(), - fanProjGeom->getDetectorWidth(), - m_iDetectorSuperSampling, m_iGPUIndex); - - } else if (fanVecProjGeom) { - - // Rescale projs to fPixelSize == 1 - float fPixelSize = pVolGeom->getPixelLengthX(); - const astraCUDA::SFanProjection* projs; - projs = fanVecProjGeom->getProjectionVectors(); - - astraCUDA::SFanProjection* scaledProjs = new astraCUDA::SFanProjection[fanVecProjGeom->getProjectionAngleCount()]; -#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 < fanVecProjGeom->getProjectionAngleCount(); ++i) { - SCALE(Src,i,1.0f/fPixelSize); - SCALE(DetS,i,1.0f/fPixelSize); - SCALE(DetU,i,1.0f/fPixelSize); + } else if (fanProjGeom || fanVecProjGeom) { + + astraCUDA::SFanProjection* projs; + float outputScale; + + if (fanProjGeom) { + ok = convertAstraGeometry(pVolGeom, fanProjGeom, projs, outputScale); + } else { + ok = convertAstraGeometry(pVolGeom, fanVecProjGeom, projs, outputScale); } + ASTRA_ASSERT(ok); ok = astraCudaFanFP(m_pVolume->getDataConst(), m_pSinogram->getData(), pVolGeom->getGridColCount(), pVolGeom->getGridRowCount(), - fanVecProjGeom->getProjectionAngleCount(), - fanVecProjGeom->getDetectorCount(), - scaledProjs, - /* 1.0f / pVolGeom->getPixelLengthX(), */ - m_iDetectorSuperSampling, m_iGPUIndex); + m_pSinogram->getGeometry()->getProjectionAngleCount(), + m_pSinogram->getGeometry()->getDetectorCount(), + projs, + m_iDetectorSuperSampling, outputScale, m_iGPUIndex); - delete[] scaledProjs; + delete[] projs; } else { diff --git a/src/CudaForwardProjectionAlgorithm3D.cpp b/src/CudaForwardProjectionAlgorithm3D.cpp index 4ebecef..8e6bab5 100644 --- a/src/CudaForwardProjectionAlgorithm3D.cpp +++ b/src/CudaForwardProjectionAlgorithm3D.cpp @@ -40,6 +40,8 @@ $Id$ #include "astra/ParallelVecProjectionGeometry3D.h" #include "astra/ConeVecProjectionGeometry3D.h" +#include "astra/Logging.h" + #include "../cuda/3d/astra3d.h" using namespace std; @@ -251,6 +253,25 @@ void CCudaForwardProjectionAlgorithm3D::run(int) projKernel = projector->getProjectionKernel(); } +#if 0 + // Debugging code that gives the coordinates of the corners of the volume + // projected on the detector. + { + float fX[] = { volgeom.getWindowMinX(), volgeom.getWindowMaxX() }; + float fY[] = { volgeom.getWindowMinY(), volgeom.getWindowMaxY() }; + float fZ[] = { volgeom.getWindowMinZ(), volgeom.getWindowMaxZ() }; + + for (int a = 0; a < projgeom->getProjectionCount(); ++a) + for (int i = 0; i < 2; ++i) + for (int j = 0; j < 2; ++j) + for (int k = 0; k < 2; ++k) { + float fU, fV; + projgeom->projectPoint(fX[i], fY[j], fZ[k], a, fU, fV); + ASTRA_DEBUG("%3d %c1,%c1,%c1 -> %12f %12f", a, i ? ' ' : '-', j ? ' ' : '-', k ? ' ' : '-', fU, fV); + } + } +#endif + if (conegeom) { astraCudaConeFP(m_pVolume->getDataConst(), m_pProjections->getData(), volgeom.getGridColCount(), diff --git a/src/CudaProjector3D.cpp b/src/CudaProjector3D.cpp index f51433d..2f4c799 100644 --- a/src/CudaProjector3D.cpp +++ b/src/CudaProjector3D.cpp @@ -28,6 +28,9 @@ $Id$ #include "astra/CudaProjector3D.h" +#include "astra/VolumeGeometry3D.h" +#include "astra/ProjectionGeometry3D.h" + namespace astra { @@ -102,21 +105,7 @@ bool CCudaProjector3D::initialize(const Config& _cfg) return false; } - // TODO: These should go to the parent. - - // ProjectionGeometry - XMLNode* node = _cfg.self->getSingleNode("ProjectionGeometry"); - // TODO: Implement - ASTRA_DELETE(node); - CC.markNodeParsed("ProjectionGeometry"); - - // ReconstructionGeometry - node = _cfg.self->getSingleNode("VolumeGeometry"); - // TODO: Implement - ASTRA_DELETE(node); - CC.markNodeParsed("VolumeGeometry"); - - node = _cfg.self->getSingleNode("ProjectionKernel"); + XMLNode* node = _cfg.self->getSingleNode("ProjectionKernel"); m_projectionKernel = ker3d_default; if (node) { std::string sProjKernel = node->getContent(); diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp index 1d66137..1c6b763 100644 --- a/src/CudaReconstructionAlgorithm2D.cpp +++ b/src/CudaReconstructionAlgorithm2D.cpp @@ -37,6 +37,8 @@ $Id$ #include "astra/FanFlatVecProjectionGeometry2D.h" #include "astra/CudaProjector2D.h" +#include "astra/Logging.h" + #include "../cuda/2d/algo.h" #include <ctime> @@ -176,7 +178,7 @@ bool CCudaReconstructionAlgorithm2D::initialize(const Config& _cfg) id = boost::lexical_cast<int>(node->getContent()); CProjector2D *projector = CProjector2DManager::getSingleton().get(id); if (!dynamic_cast<CCudaProjector2D*>(projector)) { - cout << "Warning: non-CUDA Projector2D passed" << std::endl; + ASTRA_WARN("non-CUDA Projector2D passed"); } delete node; } @@ -350,7 +352,7 @@ bool CCudaReconstructionAlgorithm2D::setupGeometry() const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); - // TODO: off-center geometry, non-square pixels + // TODO: non-square pixels? dims.iVolWidth = volgeom.getGridColCount(); dims.iVolHeight = volgeom.getGridRowCount(); float fPixelSize = volgeom.getPixelLengthX(); @@ -365,66 +367,42 @@ bool CCudaReconstructionAlgorithm2D::setupGeometry() if (parProjGeom) { + float *offsets, *angles, detSize, outputScale; + + ok = convertAstraGeometry(&volgeom, parProjGeom, offsets, angles, detSize, outputScale); + dims.iProjAngles = parProjGeom->getProjectionAngleCount(); dims.iProjDets = parProjGeom->getDetectorCount(); dims.fDetScale = parProjGeom->getDetectorWidth() / fPixelSize; ok = m_pAlgo->setGeometry(dims, parProjGeom->getProjectionAngles()); + ok &= m_pAlgo->setTOffsets(offsets); - } else if (fanProjGeom) { + // CHECKME: outputScale? detSize? - dims.iProjAngles = fanProjGeom->getProjectionAngleCount(); - dims.iProjDets = fanProjGeom->getDetectorCount(); - dims.fDetScale = fanProjGeom->getDetectorWidth() / fPixelSize; - float fOriginSourceDistance = fanProjGeom->getOriginSourceDistance(); - float fOriginDetectorDistance = fanProjGeom->getOriginDetectorDistance(); + delete[] offsets; + delete[] angles; - const float* angles = fanProjGeom->getProjectionAngles(); + } else if (fanProjGeom || fanVecProjGeom) { 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]); + float outputScale; + + if (fanProjGeom) { + ok = convertAstraGeometry(&volgeom, fanProjGeom, projs, outputScale); + } else { + ok = convertAstraGeometry(&volgeom, fanVecProjGeom, projs, outputScale); } -#undef ROTATE0 + dims.iProjAngles = m_pSinogram->getGeometry()->getProjectionAngleCount(); + dims.iProjDets = m_pSinogram->getGeometry()->getDetectorCount(); + dims.fDetScale = m_pSinogram->getGeometry()->getDetectorWidth() / fPixelSize; 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(); + // CHECKME: outputScale? - // 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; + delete[] projs; } else { @@ -441,11 +419,6 @@ bool CCudaReconstructionAlgorithm2D::setupGeometry() 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; diff --git a/src/FanFlatProjectionGeometry2D.cpp b/src/FanFlatProjectionGeometry2D.cpp index bfc7fa9..d757f18 100644 --- a/src/FanFlatProjectionGeometry2D.cpp +++ b/src/FanFlatProjectionGeometry2D.cpp @@ -204,6 +204,20 @@ CVector3D CFanFlatProjectionGeometry2D::getProjectionDirection(int _iProjectionI } //---------------------------------------------------------------------------------------- +// Get the configuration object +Config* CFanFlatProjectionGeometry2D::getConfiguration() const +{ + Config* cfg = new Config(); + cfg->initialize("ProjectionGeometry2D"); + cfg->self->addAttribute("type", "fanflat"); + cfg->self->addChildNode("DetectorCount", getDetectorCount()); + cfg->self->addChildNode("DetectorWidth", getDetectorWidth()); + cfg->self->addChildNode("DistanceOriginSource", getOriginSourceDistance()); + cfg->self->addChildNode("DistanceOriginDetector", getOriginDetectorDistance()); + cfg->self->addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount); + return cfg; +} +//---------------------------------------------------------------------------------------- } // namespace astra diff --git a/src/FanFlatVecProjectionGeometry2D.cpp b/src/FanFlatVecProjectionGeometry2D.cpp index 77f9db7..9c7b596 100644 --- a/src/FanFlatVecProjectionGeometry2D.cpp +++ b/src/FanFlatVecProjectionGeometry2D.cpp @@ -194,7 +194,7 @@ bool CFanFlatVecProjectionGeometry2D::isEqual(CProjectionGeometry2D* _pGeom2) co // Is of type bool CFanFlatVecProjectionGeometry2D::isOfType(const std::string& _sType) { - return (_sType == "fanflat_vec"); + return (_sType == "fanflat_vec"); } //---------------------------------------------------------------------------------------- @@ -227,6 +227,28 @@ bool CFanFlatVecProjectionGeometry2D::_check() //---------------------------------------------------------------------------------------- +// Get the configuration object +Config* CFanFlatVecProjectionGeometry2D::getConfiguration() const +{ + Config* cfg = new Config(); + cfg->initialize("ProjectionGeometry2D"); + cfg->self->addAttribute("type", "fanflat_vec"); + cfg->self->addChildNode("DetectorCount", getDetectorCount()); + std::string vectors = ""; + for (int i = 0; i < m_iProjectionAngleCount; ++i) { + SFanProjection& p = m_pProjectionAngles[i]; + vectors += boost::lexical_cast<string>(p.fSrcX) + ","; + vectors += boost::lexical_cast<string>(p.fSrcY) + ","; + vectors += boost::lexical_cast<string>(p.fDetSX + 0.5f * m_iDetectorCount * p.fDetUX) + ","; + vectors += boost::lexical_cast<string>(p.fDetSY + 0.5f * m_iDetectorCount * p.fDetUY) + ","; + vectors += boost::lexical_cast<string>(p.fDetUX) + ","; + vectors += boost::lexical_cast<string>(p.fDetUY); + if (i < m_iProjectionAngleCount-1) vectors += ';'; + } + cfg->self->addChildNode("Vectors", vectors); + return cfg; +} +//---------------------------------------------------------------------------------------- } // namespace astra diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp index 50cf939..4a8e5ac 100644 --- a/src/FilteredBackProjectionAlgorithm.cpp +++ b/src/FilteredBackProjectionAlgorithm.cpp @@ -39,6 +39,8 @@ $Id$ #include "astra/Fourier.h" #include "astra/DataProjector.h" +#include "astra/Logging.h" + using namespace std; namespace astra { @@ -133,7 +135,7 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) for (int i = 0; i < angleCount; i ++) { if (projectionIndex[i] > m_pProjector->getProjectionGeometry()->getProjectionAngleCount() -1 ) { - cout << "Invalid Projection Index" << endl; + ASTRA_ERROR("Invalid Projection Index"); return false; } else { int orgIndex = (int)projectionIndex[i]; diff --git a/src/GeometryUtil3D.cpp b/src/GeometryUtil3D.cpp new file mode 100644 index 0000000..52dd5a9 --- /dev/null +++ b/src/GeometryUtil3D.cpp @@ -0,0 +1,75 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp + 2014-2015, CWI, Amsterdam + +Contact: astra@uantwerpen.be +Website: http://sf.net/projects/astra-toolbox + +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/>. + +----------------------------------------------------------------------- +$Id$ +*/ + +#include "astra/GeometryUtil3D.h" + +namespace astra { + +// (See declaration in header for (mathematical) description of these functions) + + +void computeBP_UV_Coeffs(const SPar3DProjection& proj, double &fUX, double &fUY, double &fUZ, double &fUC, + double &fVX, double &fVY, double &fVZ, double &fVC) +{ + double denom = (proj.fRayX*proj.fDetUY*proj.fDetVZ - proj.fRayX*proj.fDetUZ*proj.fDetVY - proj.fRayY*proj.fDetUX*proj.fDetVZ + proj.fRayY*proj.fDetUZ*proj.fDetVX + proj.fRayZ*proj.fDetUX*proj.fDetVY - proj.fRayZ*proj.fDetUY*proj.fDetVX); + + fUX = ( - (proj.fRayY*proj.fDetVZ - proj.fRayZ*proj.fDetVY)) / denom; + fUY = ( (proj.fRayX*proj.fDetVZ - proj.fRayZ*proj.fDetVX)) / denom; + fUZ = (- (proj.fRayX*proj.fDetVY - proj.fRayY*proj.fDetVX) ) / denom; + fUC = (-(proj.fDetSY*proj.fDetVZ - proj.fDetSZ*proj.fDetVY)*proj.fRayX + (proj.fRayY*proj.fDetVZ - proj.fRayZ*proj.fDetVY)*proj.fDetSX - (proj.fRayY*proj.fDetSZ - proj.fRayZ*proj.fDetSY)*proj.fDetVX) / denom; + + fVX = ((proj.fRayY*proj.fDetUZ - proj.fRayZ*proj.fDetUY) ) / denom; + fVY = (- (proj.fRayX*proj.fDetUZ - proj.fRayZ*proj.fDetUX) ) / denom; + fVZ = ((proj.fRayX*proj.fDetUY - proj.fRayY*proj.fDetUX) ) / denom; + fVC = ((proj.fDetSY*proj.fDetUZ - proj.fDetSZ*proj.fDetUY)*proj.fRayX - (proj.fRayY*proj.fDetUZ - proj.fRayZ*proj.fDetUY)*proj.fDetSX + (proj.fRayY*proj.fDetSZ - proj.fRayZ*proj.fDetSY)*proj.fDetUX ) / denom; +} + + + +void computeBP_UV_Coeffs(const SConeProjection& proj, double &fUX, double &fUY, double &fUZ, double &fUC, + double &fVX, double &fVY, double &fVZ, double &fVC, + double &fDX, double &fDY, double &fDZ, double &fDC) +{ + fUX = (proj.fDetSZ - proj.fSrcZ)*proj.fDetVY - (proj.fDetSY - proj.fSrcY)*proj.fDetVZ; + fUY = (proj.fDetSX - proj.fSrcX)*proj.fDetVZ -(proj.fDetSZ - proj.fSrcZ)*proj.fDetVX; + fUZ = (proj.fDetSY - proj.fSrcY)*proj.fDetVX - (proj.fDetSX - proj.fSrcX)*proj.fDetVY; + fUC = (proj.fDetSY*proj.fDetVZ - proj.fDetSZ*proj.fDetVY)*proj.fSrcX - (proj.fDetSX*proj.fDetVZ - proj.fDetSZ*proj.fDetVX)*proj.fSrcY + (proj.fDetSX*proj.fDetVY - proj.fDetSY*proj.fDetVX)*proj.fSrcZ; + + fVX = (proj.fDetSY - proj.fSrcY)*proj.fDetUZ-(proj.fDetSZ - proj.fSrcZ)*proj.fDetUY; + fVY = (proj.fDetSZ - proj.fSrcZ)*proj.fDetUX - (proj.fDetSX - proj.fSrcX)*proj.fDetUZ; + fVZ = (proj.fDetSX - proj.fSrcX)*proj.fDetUY-(proj.fDetSY - proj.fSrcY)*proj.fDetUX; + fVC = -(proj.fDetSY*proj.fDetUZ - proj.fDetSZ*proj.fDetUY)*proj.fSrcX + (proj.fDetSX*proj.fDetUZ - proj.fDetSZ*proj.fDetUX)*proj.fSrcY - (proj.fDetSX*proj.fDetUY - proj.fDetSY*proj.fDetUX)*proj.fSrcZ; + + fDX = proj.fDetUY*proj.fDetVZ - proj.fDetUZ*proj.fDetVY; + fDY = proj.fDetUZ*proj.fDetVX - proj.fDetUX*proj.fDetVZ; + fDZ = proj.fDetUX*proj.fDetVY - proj.fDetUY*proj.fDetVX; + fDC = -proj.fSrcX * (proj.fDetUY*proj.fDetVZ - proj.fDetUZ*proj.fDetVY) - proj.fSrcY * (proj.fDetUZ*proj.fDetVX - proj.fDetUX*proj.fDetVZ) - proj.fSrcZ * (proj.fDetUX*proj.fDetVY - proj.fDetUY*proj.fDetVX); +} + + +} diff --git a/src/Logger.cpp b/src/Logger.cpp deleted file mode 100644 index 148e18c..0000000 --- a/src/Logger.cpp +++ /dev/null @@ -1,77 +0,0 @@ -/* ------------------------------------------------------------------------ -Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp - 2014-2015, CWI, Amsterdam - -Contact: astra@uantwerpen.be -Website: http://sf.net/projects/astra-toolbox - -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/>. - ------------------------------------------------------------------------ -$Id$ -*/ - -#include <astra/Logger.h> - -using namespace astra; - -const char * g_loggerFileName = "astra_logger.txt"; - -void CLogger::_assureIsInitialized() -{ - if(!m_bInitialized) - { - m_pOutFile = fopen(g_loggerFileName, "r"); - if(m_pOutFile != NULL) - { - // file exists, users wants to log - fclose(m_pOutFile); - m_pOutFile = fopen(g_loggerFileName, "w"); - } - - m_bInitialized = true; - } -} - -void CLogger::writeLine(const char * _text) -{ - _assureIsInitialized(); - - if(m_pOutFile != NULL) - { - fprintf(m_pOutFile, "%s\n", _text); - fflush(m_pOutFile); - } -} - -void CLogger::writeTerminalCUDAError(const char * _fileName, int _iLine, const char * _errString) -{ - char buffer[256]; - - sprintf(buffer, "Cuda error in file '%s' in line %i : %s.", _fileName, _iLine, _errString); - - writeLine(buffer); -} - -CLogger::CLogger() -{ - ; -} - -FILE * CLogger::m_pOutFile = NULL; -bool CLogger::m_bInitialized = false; diff --git a/src/Logging.cpp b/src/Logging.cpp new file mode 100644 index 0000000..8290ca0 --- /dev/null +++ b/src/Logging.cpp @@ -0,0 +1,184 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp + 2014-2015, CWI, Amsterdam + +Contact: astra@uantwerpen.be +Website: http://sf.net/projects/astra-toolbox + +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/>. + +----------------------------------------------------------------------- +$Id$ +*/ + +#define CLOG_MAIN +#include <astra/clog.h> + +#include <astra/Logging.h> + +#include <cstdio> + +using namespace astra; + +void CLogger::enableScreen() +{ + m_bEnabledScreen = true; +} + +void CLogger::enableFile() +{ + m_bEnabledFile = true; +} + +void CLogger::enable() +{ + enableScreen(); + enableFile(); +} + +void CLogger::disableScreen() +{ + m_bEnabledScreen = false; +} + +void CLogger::disableFile() +{ + m_bEnabledFile = false; +} + +void CLogger::disable() +{ + disableScreen(); + disableFile(); +} + +void CLogger::debug(const char *sfile, int sline, const char *fmt, ...) +{ + _assureIsInitialized(); + va_list ap; + va_start(ap, fmt); + if(m_bEnabledScreen) clog_debug(sfile,sline,0,fmt,ap); + if(m_bEnabledFile && m_bFileProvided) clog_debug(sfile,sline,1,fmt,ap); +} + +void CLogger::info(const char *sfile, int sline, const char *fmt, ...) +{ + _assureIsInitialized(); + va_list ap; + va_start(ap, fmt); + if(m_bEnabledScreen) clog_info(sfile,sline,0,fmt,ap); + if(m_bEnabledFile && m_bFileProvided) clog_info(sfile,sline,1,fmt,ap); +} + +void CLogger::warn(const char *sfile, int sline, const char *fmt, ...) +{ + _assureIsInitialized(); + va_list ap; + va_start(ap, fmt); + if(m_bEnabledScreen) clog_warn(sfile,sline,0,fmt,ap); + if(m_bEnabledFile && m_bFileProvided) clog_warn(sfile,sline,1,fmt,ap); +} + +void CLogger::error(const char *sfile, int sline, const char *fmt, ...) +{ + _assureIsInitialized(); + va_list ap; + va_start(ap, fmt); + if(m_bEnabledScreen) clog_error(sfile,sline,0,fmt,ap); + if(m_bEnabledFile && m_bFileProvided) clog_error(sfile,sline,1,fmt,ap); +} + +void CLogger::_setLevel(int id, log_level m_eLevel) +{ + switch(m_eLevel){ + case LOG_DEBUG: + clog_set_level(id,CLOG_DEBUG); + break; + case LOG_INFO: + clog_set_level(id,CLOG_INFO); + break; + case LOG_WARN: + clog_set_level(id,CLOG_WARN); + break; + case LOG_ERROR: + clog_set_level(id,CLOG_ERROR); + break; + } +} + +void CLogger::setOutputScreen(int fd, log_level m_eLevel) +{ + _assureIsInitialized(); + if(fd==1||fd==2){ + clog_set_fd(0, fd); + }else{ + error(__FILE__,__LINE__,"Invalid file descriptor"); + } + _setLevel(0,m_eLevel); +} + +void CLogger::setOutputFile(const char *filename, log_level m_eLevel) +{ + if(m_bFileProvided){ + clog_free(1); + m_bFileProvided=false; + } + if(!clog_init_path(1,filename)){ + m_bFileProvided=true; + _setLevel(1,m_eLevel); + } +} + +void CLogger::_assureIsInitialized() +{ + if(!m_bInitialized) + { + clog_init_fd(0, 2); + clog_set_level(0, CLOG_INFO); + clog_set_fmt(0, "%l: %m\n"); + m_bInitialized = true; + } +} + +void CLogger::setFormatFile(const char *fmt) +{ + if(m_bFileProvided){ + clog_set_fmt(1,fmt); + }else{ + error(__FILE__,__LINE__,"No log file specified"); + } +} +void CLogger::setFormatScreen(const char *fmt) +{ + clog_set_fmt(0,fmt); +} + +CLogger::CLogger() +{ + ; +} + +bool CLogger::setCallbackScreen(void (*cb)(const char *msg, size_t len)){ + _assureIsInitialized(); + return clog_set_cb(0,cb)==0; +} + +bool CLogger::m_bEnabledScreen = true; +bool CLogger::m_bEnabledFile = true; +bool CLogger::m_bFileProvided = false; +bool CLogger::m_bInitialized = false; diff --git a/src/ParallelProjectionGeometry2D.cpp b/src/ParallelProjectionGeometry2D.cpp index 8f8faa4..cac8f30 100644 --- a/src/ParallelProjectionGeometry2D.cpp +++ b/src/ParallelProjectionGeometry2D.cpp @@ -27,6 +27,7 @@ $Id$ */ #include "astra/ParallelProjectionGeometry2D.h" +#include <boost/lexical_cast.hpp> #include <cstring> @@ -168,6 +169,19 @@ bool CParallelProjectionGeometry2D::isOfType(const std::string& _sType) { return (_sType == "parallel"); } + +//---------------------------------------------------------------------------------------- +// Get the configuration object +Config* CParallelProjectionGeometry2D::getConfiguration() const +{ + Config* cfg = new Config(); + cfg->initialize("ProjectionGeometry2D"); + cfg->self->addAttribute("type", "parallel"); + cfg->self->addChildNode("DetectorCount", getDetectorCount()); + cfg->self->addChildNode("DetectorWidth", getDetectorWidth()); + cfg->self->addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount); + return cfg; +} //---------------------------------------------------------------------------------------- CVector3D CParallelProjectionGeometry2D::getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex /* = 0 */) diff --git a/src/ParallelProjectionGeometry3D.cpp b/src/ParallelProjectionGeometry3D.cpp index 104fc6e..eb200f9 100644 --- a/src/ParallelProjectionGeometry3D.cpp +++ b/src/ParallelProjectionGeometry3D.cpp @@ -27,6 +27,7 @@ $Id$ */ #include "astra/ParallelProjectionGeometry3D.h" +#include <boost/lexical_cast.hpp> #include <cstring> @@ -154,19 +155,24 @@ bool CParallelProjectionGeometry3D::isEqual(const CProjectionGeometry3D * _pGeom // is of type bool CParallelProjectionGeometry3D::isOfType(const std::string& _sType) const { - return (_sType == "parallel"); + return (_sType == "parallel3d"); } //---------------------------------------------------------------------------------------- -void CParallelProjectionGeometry3D::toXML(XMLNode* _sNode) const +// Get the configuration object +Config* CParallelProjectionGeometry3D::getConfiguration() const { - _sNode->addAttribute("type","parallel3d"); - _sNode->addChildNode("DetectorSpacingX", m_fDetectorSpacingX); - _sNode->addChildNode("DetectorSpacingY", m_fDetectorSpacingY); - _sNode->addChildNode("DetectorRowCount", m_iDetectorRowCount); - _sNode->addChildNode("DetectorColCount", m_iDetectorColCount); - _sNode->addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount); + Config* cfg = new Config(); + cfg->initialize("ProjectionGeometry3D"); + cfg->self->addAttribute("type", "parallel3d"); + cfg->self->addChildNode("DetectorRowCount", m_iDetectorRowCount); + cfg->self->addChildNode("DetectorColCount", m_iDetectorColCount); + cfg->self->addChildNode("DetectorSpacingX", m_fDetectorSpacingX); + cfg->self->addChildNode("DetectorSpacingY", m_fDetectorSpacingY); + cfg->self->addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount); + return cfg; } +//---------------------------------------------------------------------------------------- CVector3D CParallelProjectionGeometry3D::getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex) const { @@ -179,6 +185,22 @@ CVector3D CParallelProjectionGeometry3D::getProjectionDirection(int _iProjection return CVector3D(fDirX, fDirY, fDirZ); } +void CParallelProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, + int iAngleIndex, + float32 &fU, float32 &fV) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + // V (detector row) + fV = detectorOffsetYToRowIndexFloat(fZ); + + // U (detector column) + float alpha = m_pfProjectionAngles[iAngleIndex]; + // projector direction is (cos(alpha), sin(alpha)) + fU = detectorOffsetXToColIndexFloat(cos(alpha) * fX + sin(alpha) * fY); +} + CParallelProjectionGeometry2D * CParallelProjectionGeometry3D::createProjectionGeometry2D() const { const float32 * pfProjectionAngles = getProjectionAngles(); //new float32[getProjectionCount()]; diff --git a/src/ParallelVecProjectionGeometry3D.cpp b/src/ParallelVecProjectionGeometry3D.cpp index c8071a6..dc325e9 100644 --- a/src/ParallelVecProjectionGeometry3D.cpp +++ b/src/ParallelVecProjectionGeometry3D.cpp @@ -202,14 +202,38 @@ bool CParallelVecProjectionGeometry3D::isOfType(const std::string& _sType) const } //---------------------------------------------------------------------------------------- -void CParallelVecProjectionGeometry3D::toXML(XMLNode* _sNode) const +// Get the configuration object +Config* CParallelVecProjectionGeometry3D::getConfiguration() const { - _sNode->addAttribute("type","parallel3d_vec"); - _sNode->addChildNode("DetectorRowCount", m_iDetectorRowCount); - _sNode->addChildNode("DetectorColCount", m_iDetectorColCount); - // TODO: - //_sNode->addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount); + Config* cfg = new Config(); + cfg->initialize("ProjectionGeometry3D"); + + cfg->self->addAttribute("type", "parallel3d_vec"); + cfg->self->addChildNode("DetectorRowCount", m_iDetectorRowCount); + cfg->self->addChildNode("DetectorColCount", m_iDetectorColCount); + + std::string vectors = ""; + for (int i = 0; i < m_iProjectionAngleCount; ++i) { + SPar3DProjection& p = m_pProjectionAngles[i]; + vectors += boost::lexical_cast<string>(p.fRayX) + ","; + vectors += boost::lexical_cast<string>(p.fRayY) + ","; + vectors += boost::lexical_cast<string>(p.fRayZ) + ","; + vectors += boost::lexical_cast<string>(p.fDetSX + 0.5f*m_iDetectorRowCount*p.fDetVX + 0.5f*m_iDetectorColCount*p.fDetUX) + ","; + vectors += boost::lexical_cast<string>(p.fDetSY + 0.5f*m_iDetectorRowCount*p.fDetVY + 0.5f*m_iDetectorColCount*p.fDetUY) + ","; + vectors += boost::lexical_cast<string>(p.fDetSZ + 0.5f*m_iDetectorRowCount*p.fDetVZ + 0.5f*m_iDetectorColCount*p.fDetUZ) + ","; + vectors += boost::lexical_cast<string>(p.fDetUX) + ","; + vectors += boost::lexical_cast<string>(p.fDetUY) + ","; + vectors += boost::lexical_cast<string>(p.fDetUZ) + ","; + vectors += boost::lexical_cast<string>(p.fDetVX) + ","; + vectors += boost::lexical_cast<string>(p.fDetVY) + ","; + vectors += boost::lexical_cast<string>(p.fDetVZ); + if (i < m_iProjectionAngleCount-1) vectors += ';'; + } + cfg->self->addChildNode("Vectors", vectors); + + return cfg; } +//---------------------------------------------------------------------------------------- CVector3D CParallelVecProjectionGeometry3D::getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex) const { @@ -218,6 +242,24 @@ CVector3D CParallelVecProjectionGeometry3D::getProjectionDirection(int _iProject return CVector3D(p.fRayX, p.fRayY, p.fRayZ); } +void CParallelVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, + int iAngleIndex, + float32 &fU, float32 &fV) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + double fUX, fUY, fUZ, fUC; + double fVX, fVY, fVZ, fVC; + + computeBP_UV_Coeffs(m_pProjectionAngles[iAngleIndex], + fUX, fUY, fUZ, fUC, fVX, fVY, fVZ, fVC); + + // The -0.5f shifts from corner to center of detector pixels + fU = (fUX*fX + fUY*fY + fUZ*fZ + fUC) - 0.5f; + fV = (fVX*fX + fVY*fY + fVZ*fZ + fVC) - 0.5f; + +} //---------------------------------------------------------------------------------------- diff --git a/src/Projector3D.cpp b/src/Projector3D.cpp index 0f5edb7..b546ee9 100644 --- a/src/Projector3D.cpp +++ b/src/Projector3D.cpp @@ -28,6 +28,12 @@ $Id$ #include "astra/Projector3D.h" +#include "astra/VolumeGeometry3D.h" +#include "astra/ParallelProjectionGeometry3D.h" +#include "astra/ParallelVecProjectionGeometry3D.h" +#include "astra/ConeProjectionGeometry3D.h" +#include "astra/ConeVecProjectionGeometry3D.h" + namespace astra { @@ -84,6 +90,37 @@ bool CProjector3D::_check() bool CProjector3D::initialize(const Config& _cfg) { assert(_cfg.self); + ConfigStackCheck<CProjector3D> CC("Projector3D", this, _cfg); + + XMLNode* node; + + node = _cfg.self->getSingleNode("ProjectionGeometry"); + ASTRA_CONFIG_CHECK(node, "Projector3D", "No ProjectionGeometry tag specified."); + std::string type = node->getAttribute("type"); + CProjectionGeometry3D* pProjGeometry = 0; + if (type == "parallel3d") { + pProjGeometry = new CParallelProjectionGeometry3D(); + } else if (type == "parallel3d_vec") { + pProjGeometry = new CParallelVecProjectionGeometry3D(); + } else if (type == "cone") { + pProjGeometry = new CConeProjectionGeometry3D(); + } else if (type == "cone_vec") { + pProjGeometry = new CConeVecProjectionGeometry3D(); + } else { + // Invalid geometry type + } + pProjGeometry->initialize(Config(node)); // this deletes node + m_pProjectionGeometry = pProjGeometry; + ASTRA_CONFIG_CHECK(m_pProjectionGeometry->isInitialized(), "Projector3D", "ProjectionGeometry not initialized."); + CC.markNodeParsed("ProjectionGeometry"); + + node = _cfg.self->getSingleNode("VolumeGeometry"); + ASTRA_CONFIG_CHECK(node, "Projector3D", "No VolumeGeometry tag specified."); + CVolumeGeometry3D* pVolGeometry = new CVolumeGeometry3D(); + pVolGeometry->initialize(Config(node)); // this deletes node + m_pVolumeGeometry = pVolGeometry; + ASTRA_CONFIG_CHECK(m_pVolumeGeometry->isInitialized(), "Projector2D", "VolumeGeometry not initialized."); + CC.markNodeParsed("VolumeGeometry"); return true; } diff --git a/src/ReconstructionAlgorithm3D.cpp b/src/ReconstructionAlgorithm3D.cpp index 0ce9777..13d069d 100644 --- a/src/ReconstructionAlgorithm3D.cpp +++ b/src/ReconstructionAlgorithm3D.cpp @@ -145,7 +145,6 @@ bool CReconstructionAlgorithm3D::initialize(const Config& _cfg) id = boost::lexical_cast<int>(_cfg.self->getOption("SinogramMaskId")); m_pSinogramMask = dynamic_cast<CFloat32ProjectionData3D*>(CData3DManager::getSingleton().get(id)); } - CC.markOptionParsed("SinogramMaskId"); // Constraints - NEW if (_cfg.self->hasOption("MinConstraint")) { diff --git a/src/SparseMatrixProjectionGeometry2D.cpp b/src/SparseMatrixProjectionGeometry2D.cpp index 7f14d9d..86357d2 100644 --- a/src/SparseMatrixProjectionGeometry2D.cpp +++ b/src/SparseMatrixProjectionGeometry2D.cpp @@ -189,7 +189,20 @@ bool CSparseMatrixProjectionGeometry2D::isOfType(const std::string& _sType) return (_sType == "sparse_matrix"); } //---------------------------------------------------------------------------------------- +// Get the configuration object +Config* CSparseMatrixProjectionGeometry2D::getConfiguration() const +{ + Config* cfg = new Config(); + cfg->initialize("ProjectionGeometry2D"); + cfg->self->addAttribute("type", "sparse matrix"); + cfg->self->addChildNode("DetectorCount", getDetectorCount()); + cfg->self->addChildNode("DetectorWidth", getDetectorWidth()); + cfg->self->addChildNode("ProjectionAngles", m_pfProjectionAngles, m_iProjectionAngleCount); + cfg->self->addChildNode("MatrixID", CMatrixManager::getSingleton().getIndex(m_pMatrix)); + return cfg; +} +//---------------------------------------------------------------------------------------- CVector3D CSparseMatrixProjectionGeometry2D::getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex) { CVector3D vOutput(0.0f, 0.0f, 0.0f); diff --git a/src/VolumeGeometry2D.cpp b/src/VolumeGeometry2D.cpp index 19aa905..d412914 100644 --- a/src/VolumeGeometry2D.cpp +++ b/src/VolumeGeometry2D.cpp @@ -152,7 +152,7 @@ CVolumeGeometry2D* CVolumeGeometry2D::clone() } //---------------------------------------------------------------------------------------- -// Initialization witha COnfig object +// Initialization witha Config object bool CVolumeGeometry2D::initialize(const Config& _cfg) { ASTRA_ASSERT(_cfg.self); @@ -277,6 +277,24 @@ bool CVolumeGeometry2D::isEqual(CVolumeGeometry2D* _pGeom2) const return true; } + +//---------------------------------------------------------------------------------------- +// Get the configuration object +Config* CVolumeGeometry2D::getConfiguration() const +{ + Config* cfg = new Config(); + cfg->initialize("VolumeGeometry2D"); + + cfg->self->addChildNode("GridColCount", m_iGridColCount); + cfg->self->addChildNode("GridRowCount", m_iGridRowCount); + + cfg->self->addOption("WindowMinX", m_fWindowMinX); + cfg->self->addOption("WindowMaxX", m_fWindowMaxX); + cfg->self->addOption("WindowMinY", m_fWindowMinY); + cfg->self->addOption("WindowMaxY", m_fWindowMaxY); + + return cfg; +} //---------------------------------------------------------------------------------------- } // namespace astra diff --git a/src/VolumeGeometry3D.cpp b/src/VolumeGeometry3D.cpp index d7a93a9..66e6f0c 100644 --- a/src/VolumeGeometry3D.cpp +++ b/src/VolumeGeometry3D.cpp @@ -380,5 +380,28 @@ CVolumeGeometry2D * CVolumeGeometry3D::createVolumeGeometry2D() const } //---------------------------------------------------------------------------------------- +// Get the configuration object +Config* CVolumeGeometry3D::getConfiguration() const +{ + Config* cfg = new Config(); + cfg->initialize("VolumeGeometry3D"); + + cfg->self->addChildNode("GridColCount", m_iGridColCount); + cfg->self->addChildNode("GridRowCount", m_iGridRowCount); + cfg->self->addChildNode("GridSliceCount", m_iGridSliceCount); + + cfg->self->addOption("WindowMinX", m_fWindowMinX); + cfg->self->addOption("WindowMaxX", m_fWindowMaxX); + cfg->self->addOption("WindowMinY", m_fWindowMinY); + cfg->self->addOption("WindowMaxY", m_fWindowMaxY); + cfg->self->addOption("WindowMinZ", m_fWindowMinZ); + cfg->self->addOption("WindowMaxZ", m_fWindowMaxZ); + + return cfg; +} +//---------------------------------------------------------------------------------------- + + +//---------------------------------------------------------------------------------------- } // namespace astra diff --git a/src/XMLDocument.cpp b/src/XMLDocument.cpp index b39875b..406564f 100644 --- a/src/XMLDocument.cpp +++ b/src/XMLDocument.cpp @@ -109,4 +109,12 @@ void XMLDocument::saveToFile(string sFilename) } //----------------------------------------------------------------------------- +std::string XMLDocument::toString() +{ + std::stringstream ss; + ss << *fDOMDocument->first_node(); + return ss.str(); +} + +//----------------------------------------------------------------------------- diff --git a/src/XMLNode.cpp b/src/XMLNode.cpp index a5b6796..4b2bdf4 100644 --- a/src/XMLNode.cpp +++ b/src/XMLNode.cpp @@ -448,13 +448,11 @@ void XMLNode::setContent(float32 _fValue) // Set content - LIST void XMLNode::setContent(float32* pfList, int _iSize) { - addAttribute("listsize", _iSize); - for (int i = 0; i < _iSize; i++) { - XMLNode* item = addChildNode("ListItem"); - item->addAttribute("index", i); - item->addAttribute("value",pfList[i]); - delete item; + std::string str = (_iSize > 0) ? boost::lexical_cast<std::string>(pfList[0]) : ""; + for (int i = 1; i < _iSize; i++) { + str += "," + boost::lexical_cast<std::string>(pfList[i]); } + setContent(str); } //----------------------------------------------------------------------------- |