From 5304d08cd1ab7b8d778c367912934376eb92370f Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 9 Mar 2015 15:43:56 +0100 Subject: Allow non-centered volume geometry in SIRT3D and CGLS3D --- src/CudaCglsAlgorithm3D.cpp | 39 +-------------------------------------- src/CudaSirtAlgorithm3D.cpp | 38 +------------------------------------- 2 files changed, 2 insertions(+), 75 deletions(-) (limited to 'src') diff --git a/src/CudaCglsAlgorithm3D.cpp b/src/CudaCglsAlgorithm3D.cpp index a5500d6..3677458 100644 --- a/src/CudaCglsAlgorithm3D.cpp +++ b/src/CudaCglsAlgorithm3D.cpp @@ -171,9 +171,6 @@ void CCudaCglsAlgorithm3D::run(int _iNrIterations) ASTRA_ASSERT(m_bIsInitialized); const CProjectionGeometry3D* projgeom = m_pSinogram->getGeometry(); - const CConeProjectionGeometry3D* conegeom = dynamic_cast(projgeom); - const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast(projgeom); - const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast(projgeom); const CVolumeGeometry3D& volgeom = *m_pReconstruction->getGeometry(); bool ok = true; @@ -182,41 +179,7 @@ void CCudaCglsAlgorithm3D::run(int _iNrIterations) ok &= m_pCgls->setGPUIndex(m_iGPUIndex); - ok &= m_pCgls->setReconstructionGeometry(volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount()); -/* - unsigned int iProjAngles, - unsigned int iProjU, - unsigned int iProjV, - float fOriginSourceDistance, - float fOriginDetectorDistance, - float fDetUSize, - float fDetVSize, - const float *pfAngles) -*/ - if (conegeom) { - ok &= m_pCgls->setConeGeometry(conegeom->getProjectionCount(), - conegeom->getDetectorColCount(), - conegeom->getDetectorRowCount(), - conegeom->getOriginSourceDistance(), - conegeom->getOriginDetectorDistance(), - conegeom->getDetectorSpacingX(), - conegeom->getDetectorSpacingY(), - conegeom->getProjectionAngles()); - } else if (parvec3dgeom) { - ok &= m_pCgls->setPar3DGeometry(parvec3dgeom->getProjectionCount(), - parvec3dgeom->getDetectorColCount(), - parvec3dgeom->getDetectorRowCount(), - parvec3dgeom->getProjectionVectors()); - } else if (conevec3dgeom) { - ok &= m_pCgls->setConeGeometry(conevec3dgeom->getProjectionCount(), - conevec3dgeom->getDetectorColCount(), - conevec3dgeom->getDetectorRowCount(), - conevec3dgeom->getProjectionVectors()); - } else { - ASTRA_ASSERT(false); - } + ok &= m_pCgls->setGeometry(&volgeom, projgeom); ok &= m_pCgls->enableSuperSampling(m_iVoxelSuperSampling, m_iDetectorSuperSampling); diff --git a/src/CudaSirtAlgorithm3D.cpp b/src/CudaSirtAlgorithm3D.cpp index da83c7e..d67778f 100644 --- a/src/CudaSirtAlgorithm3D.cpp +++ b/src/CudaSirtAlgorithm3D.cpp @@ -172,10 +172,6 @@ void CCudaSirtAlgorithm3D::run(int _iNrIterations) ASTRA_ASSERT(m_bIsInitialized); const CProjectionGeometry3D* projgeom = m_pSinogram->getGeometry(); - const CConeProjectionGeometry3D* conegeom = dynamic_cast(projgeom); - const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast(projgeom); - const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast(projgeom); - const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast(projgeom); const CVolumeGeometry3D& volgeom = *m_pReconstruction->getGeometry(); bool ok = true; @@ -184,39 +180,7 @@ void CCudaSirtAlgorithm3D::run(int _iNrIterations) ok &= m_pSirt->setGPUIndex(m_iGPUIndex); - ok &= m_pSirt->setReconstructionGeometry(volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount()); - - if (conegeom) { - ok &= m_pSirt->setConeGeometry(conegeom->getProjectionCount(), - conegeom->getDetectorColCount(), - conegeom->getDetectorRowCount(), - conegeom->getOriginSourceDistance(), - conegeom->getOriginDetectorDistance(), - conegeom->getDetectorSpacingX(), - conegeom->getDetectorSpacingY(), - conegeom->getProjectionAngles()); - } else if (par3dgeom) { - ok &= m_pSirt->setPar3DGeometry(par3dgeom->getProjectionCount(), - par3dgeom->getDetectorColCount(), - par3dgeom->getDetectorRowCount(), - par3dgeom->getDetectorSpacingX(), - par3dgeom->getDetectorSpacingY(), - par3dgeom->getProjectionAngles()); - } else if (parvec3dgeom) { - ok &= m_pSirt->setPar3DGeometry(parvec3dgeom->getProjectionCount(), - parvec3dgeom->getDetectorColCount(), - parvec3dgeom->getDetectorRowCount(), - parvec3dgeom->getProjectionVectors()); - } else if (conevec3dgeom) { - ok &= m_pSirt->setConeGeometry(conevec3dgeom->getProjectionCount(), - conevec3dgeom->getDetectorColCount(), - conevec3dgeom->getDetectorRowCount(), - conevec3dgeom->getProjectionVectors()); - } else { - ASTRA_ASSERT(false); - } + ok &= m_pSirt->setGeometry(&volgeom, projgeom); ok &= m_pSirt->enableSuperSampling(m_iVoxelSuperSampling, m_iDetectorSuperSampling); -- cgit v1.2.3 From 140f64028a6c06895ba7dad8997e14b7a05aadab Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 11 Mar 2015 12:07:48 +0100 Subject: Let astraCudaFDK use utility functions --- src/CudaFDKAlgorithm3D.cpp | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) (limited to 'src') diff --git a/src/CudaFDKAlgorithm3D.cpp b/src/CudaFDKAlgorithm3D.cpp index 7638696..0a46ff6 100644 --- a/src/CudaFDKAlgorithm3D.cpp +++ b/src/CudaFDKAlgorithm3D.cpp @@ -171,17 +171,7 @@ void CCudaFDKAlgorithm3D::run(int _iNrIterations) bool ok = true; ok = astraCudaFDK(pReconMem->getData(), pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - conegeom->getProjectionCount(), - conegeom->getDetectorColCount(), - conegeom->getDetectorRowCount(), - conegeom->getOriginSourceDistance(), - conegeom->getOriginDetectorDistance(), - conegeom->getDetectorSpacingX(), - conegeom->getDetectorSpacingY(), - conegeom->getProjectionAngles(), + &volgeom, conegeom, m_bShortScan, m_iGPUIndex, m_iVoxelSuperSampling); ASTRA_ASSERT(ok); -- cgit v1.2.3 From 18d12242207d1113c3015b451f522531168e626a Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 11 Mar 2015 17:27:44 +0100 Subject: Add flexible volgeom3d support to astraCudaBP_SIRTWeighted --- src/CudaBackProjectionAlgorithm3D.cpp | 87 +++++++++++------------------------ 1 file changed, 28 insertions(+), 59 deletions(-) (limited to 'src') diff --git a/src/CudaBackProjectionAlgorithm3D.cpp b/src/CudaBackProjectionAlgorithm3D.cpp index abcf096..7117cfc 100644 --- a/src/CudaBackProjectionAlgorithm3D.cpp +++ b/src/CudaBackProjectionAlgorithm3D.cpp @@ -107,16 +107,8 @@ bool CCudaBackProjectionAlgorithm3D::initialize(const Config& _cfg) m_iVoxelSuperSampling = (int)_cfg.self->getOptionNumerical("VoxelSuperSampling", 1); CC.markOptionParsed("VoxelSuperSampling"); - CFloat32ProjectionData3DMemory* pSinoMem = dynamic_cast(m_pSinogram); - ASTRA_ASSERT(pSinoMem); - const CProjectionGeometry3D* projgeom = pSinoMem->getGeometry(); -const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast(projgeom); - const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast(projgeom); - if (parvec3dgeom || par3dgeom) { - // This option is only supported for Par3D currently - m_bSIRTWeighting = _cfg.self->getOptionBool("SIRTWeighting", false); - CC.markOptionParsed("SIRTWeighting"); - } + m_bSIRTWeighting = _cfg.self->getOptionBool("SIRTWeighting", false); + CC.markOptionParsed("SIRTWeighting"); // success m_bIsInitialized = _check(); @@ -178,7 +170,12 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations) const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast(projgeom); const CVolumeGeometry3D& volgeom = *pReconMem->getGeometry(); - if (conegeom) { + if (m_bSIRTWeighting) { + astraCudaBP_SIRTWeighted(pReconMem->getData(), + pSinoMem->getDataConst(), + &volgeom, projgeom, + m_iGPUIndex, m_iVoxelSuperSampling); + } else if (conegeom) { astraCudaConeBP(pReconMem->getData(), pSinoMem->getDataConst(), volgeom.getGridColCount(), volgeom.getGridRowCount(), @@ -193,55 +190,27 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations) conegeom->getProjectionAngles(), m_iGPUIndex, m_iVoxelSuperSampling); } else if (par3dgeom) { - if (!m_bSIRTWeighting) { - astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - par3dgeom->getProjectionCount(), - par3dgeom->getDetectorColCount(), - par3dgeom->getDetectorRowCount(), - par3dgeom->getDetectorSpacingX(), - par3dgeom->getDetectorSpacingY(), - par3dgeom->getProjectionAngles(), - m_iGPUIndex, m_iVoxelSuperSampling); - } else { - astraCudaPar3DBP_SIRTWeighted(pReconMem->getData(), - pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - par3dgeom->getProjectionCount(), - par3dgeom->getDetectorColCount(), - par3dgeom->getDetectorRowCount(), - par3dgeom->getDetectorSpacingX(), - par3dgeom->getDetectorSpacingY(), - par3dgeom->getProjectionAngles(), - m_iGPUIndex, m_iVoxelSuperSampling); - } + astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(), + volgeom.getGridColCount(), + volgeom.getGridRowCount(), + volgeom.getGridSliceCount(), + par3dgeom->getProjectionCount(), + par3dgeom->getDetectorColCount(), + par3dgeom->getDetectorRowCount(), + par3dgeom->getDetectorSpacingX(), + par3dgeom->getDetectorSpacingY(), + par3dgeom->getProjectionAngles(), + m_iGPUIndex, m_iVoxelSuperSampling); } else if (parvec3dgeom) { - if (!m_bSIRTWeighting) { - astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - parvec3dgeom->getProjectionCount(), - parvec3dgeom->getDetectorColCount(), - parvec3dgeom->getDetectorRowCount(), - parvec3dgeom->getProjectionVectors(), - m_iGPUIndex, m_iVoxelSuperSampling); - } else { - astraCudaPar3DBP_SIRTWeighted(pReconMem->getData(), - pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - parvec3dgeom->getProjectionCount(), - parvec3dgeom->getDetectorColCount(), - parvec3dgeom->getDetectorRowCount(), - parvec3dgeom->getProjectionVectors(), - m_iGPUIndex, m_iVoxelSuperSampling); - } + astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(), + volgeom.getGridColCount(), + volgeom.getGridRowCount(), + volgeom.getGridSliceCount(), + parvec3dgeom->getProjectionCount(), + parvec3dgeom->getDetectorColCount(), + parvec3dgeom->getDetectorRowCount(), + parvec3dgeom->getProjectionVectors(), + m_iGPUIndex, m_iVoxelSuperSampling); } else if (conevecgeom) { astraCudaConeBP(pReconMem->getData(), pSinoMem->getDataConst(), volgeom.getGridColCount(), -- cgit v1.2.3 From 6909836555afe155ffc3897ef2189ed0562bb045 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 11 Mar 2015 18:44:53 +0100 Subject: Add flexible volgeom3d support to astraCudaBP --- src/CudaBackProjectionAlgorithm3D.cpp | 54 ++--------------------------------- 1 file changed, 3 insertions(+), 51 deletions(-) (limited to 'src') diff --git a/src/CudaBackProjectionAlgorithm3D.cpp b/src/CudaBackProjectionAlgorithm3D.cpp index 7117cfc..a8a1b0a 100644 --- a/src/CudaBackProjectionAlgorithm3D.cpp +++ b/src/CudaBackProjectionAlgorithm3D.cpp @@ -164,10 +164,6 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations) ASTRA_ASSERT(pReconMem); const CProjectionGeometry3D* projgeom = pSinoMem->getGeometry(); - const CConeProjectionGeometry3D* conegeom = dynamic_cast(projgeom); - const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast(projgeom); - const CConeVecProjectionGeometry3D* conevecgeom = dynamic_cast(projgeom); - const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast(projgeom); const CVolumeGeometry3D& volgeom = *pReconMem->getGeometry(); if (m_bSIRTWeighting) { @@ -175,54 +171,10 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations) pSinoMem->getDataConst(), &volgeom, projgeom, m_iGPUIndex, m_iVoxelSuperSampling); - } else if (conegeom) { - astraCudaConeBP(pReconMem->getData(), pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - conegeom->getProjectionCount(), - conegeom->getDetectorColCount(), - conegeom->getDetectorRowCount(), - conegeom->getOriginSourceDistance(), - conegeom->getOriginDetectorDistance(), - conegeom->getDetectorSpacingX(), - conegeom->getDetectorSpacingY(), - conegeom->getProjectionAngles(), - m_iGPUIndex, m_iVoxelSuperSampling); - } else if (par3dgeom) { - astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - par3dgeom->getProjectionCount(), - par3dgeom->getDetectorColCount(), - par3dgeom->getDetectorRowCount(), - par3dgeom->getDetectorSpacingX(), - par3dgeom->getDetectorSpacingY(), - par3dgeom->getProjectionAngles(), - m_iGPUIndex, m_iVoxelSuperSampling); - } else if (parvec3dgeom) { - astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - parvec3dgeom->getProjectionCount(), - parvec3dgeom->getDetectorColCount(), - parvec3dgeom->getDetectorRowCount(), - parvec3dgeom->getProjectionVectors(), - m_iGPUIndex, m_iVoxelSuperSampling); - } else if (conevecgeom) { - astraCudaConeBP(pReconMem->getData(), pSinoMem->getDataConst(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - conevecgeom->getProjectionCount(), - conevecgeom->getDetectorColCount(), - conevecgeom->getDetectorRowCount(), - conevecgeom->getProjectionVectors(), - m_iGPUIndex, m_iVoxelSuperSampling); } else { - ASTRA_ASSERT(false); + astraCudaBP(pReconMem->getData(), pSinoMem->getDataConst(), + &volgeom, projgeom, + m_iGPUIndex, m_iVoxelSuperSampling); } } -- cgit v1.2.3 From 57ee6b85884b8226b26b7415ef151b4a6e63337c Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 12 Mar 2015 11:53:40 +0100 Subject: Add flexible volgeom3d support to astraCudaFP --- src/CudaForwardProjectionAlgorithm3D.cpp | 59 ++------------------------------ 1 file changed, 3 insertions(+), 56 deletions(-) (limited to 'src') diff --git a/src/CudaForwardProjectionAlgorithm3D.cpp b/src/CudaForwardProjectionAlgorithm3D.cpp index bb122e0..914ee2f 100644 --- a/src/CudaForwardProjectionAlgorithm3D.cpp +++ b/src/CudaForwardProjectionAlgorithm3D.cpp @@ -239,10 +239,6 @@ void CCudaForwardProjectionAlgorithm3D::run(int) assert(m_bIsInitialized); const CProjectionGeometry3D* projgeom = m_pProjections->getGeometry(); - const CConeProjectionGeometry3D* conegeom = dynamic_cast(projgeom); - const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast(projgeom); - const CConeVecProjectionGeometry3D* conevecgeom = dynamic_cast(projgeom); - const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast(projgeom); const CVolumeGeometry3D& volgeom = *m_pVolume->getGeometry(); Cuda3DProjectionKernel projKernel = ker3d_default; @@ -270,58 +266,9 @@ void CCudaForwardProjectionAlgorithm3D::run(int) } #endif - if (conegeom) { - astraCudaConeFP(m_pVolume->getDataConst(), m_pProjections->getData(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - conegeom->getProjectionCount(), - conegeom->getDetectorColCount(), - conegeom->getDetectorRowCount(), - conegeom->getOriginSourceDistance(), - conegeom->getOriginDetectorDistance(), - conegeom->getDetectorSpacingX(), - conegeom->getDetectorSpacingY(), - conegeom->getProjectionAngles(), - m_iGPUIndex, m_iDetectorSuperSampling); - } else if (par3dgeom) { - astraCudaPar3DFP(m_pVolume->getDataConst(), m_pProjections->getData(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - par3dgeom->getProjectionCount(), - par3dgeom->getDetectorColCount(), - par3dgeom->getDetectorRowCount(), - par3dgeom->getDetectorSpacingX(), - par3dgeom->getDetectorSpacingY(), - par3dgeom->getProjectionAngles(), - m_iGPUIndex, m_iDetectorSuperSampling, - projKernel); - } else if (parvec3dgeom) { - astraCudaPar3DFP(m_pVolume->getDataConst(), m_pProjections->getData(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - parvec3dgeom->getProjectionCount(), - parvec3dgeom->getDetectorColCount(), - parvec3dgeom->getDetectorRowCount(), - parvec3dgeom->getProjectionVectors(), - m_iGPUIndex, m_iDetectorSuperSampling, - projKernel); - } else if (conevecgeom) { - astraCudaConeFP(m_pVolume->getDataConst(), m_pProjections->getData(), - volgeom.getGridColCount(), - volgeom.getGridRowCount(), - volgeom.getGridSliceCount(), - conevecgeom->getProjectionCount(), - conevecgeom->getDetectorColCount(), - conevecgeom->getDetectorRowCount(), - conevecgeom->getProjectionVectors(), - m_iGPUIndex, m_iDetectorSuperSampling); - } else { - ASTRA_ASSERT(false); - } - + astraCudaFP(m_pVolume->getDataConst(), m_pProjections->getData(), + &volgeom, projgeom, + m_iGPUIndex, m_iDetectorSuperSampling, projKernel); } -- cgit v1.2.3 From 18b6d25f7e4f0943b3592f3bb4f6ca5ed9c285d3 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Fri, 19 Jun 2015 22:28:06 +0200 Subject: Add support for Python algorithm plugins --- src/PluginAlgorithm.cpp | 294 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 294 insertions(+) create mode 100644 src/PluginAlgorithm.cpp (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp new file mode 100644 index 0000000..df13f31 --- /dev/null +++ b/src/PluginAlgorithm.cpp @@ -0,0 +1,294 @@ +/* +----------------------------------------------------------------------- +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 . + +----------------------------------------------------------------------- +$Id$ +*/ + +#ifdef ASTRA_PYTHON + +#include "astra/PluginAlgorithm.h" +#include +#include +#include +#include +#include +#include + +namespace astra { + +CPluginAlgorithm::CPluginAlgorithm(PyObject* pyclass){ + instance = PyObject_CallObject(pyclass, NULL); +} + +CPluginAlgorithm::~CPluginAlgorithm(){ + if(instance!=NULL){ + Py_DECREF(instance); + instance = NULL; + } +} + +bool CPluginAlgorithm::initialize(const Config& _cfg){ + if(instance==NULL) return false; + PyObject *cfgDict = XMLNode2dict(_cfg.self); + PyObject *retVal = PyObject_CallMethod(instance, "astra_init", "O",cfgDict); + Py_DECREF(cfgDict); + if(retVal==NULL) return false; + m_bIsInitialized = true; + Py_DECREF(retVal); + return m_bIsInitialized; +} + +void CPluginAlgorithm::run(int _iNrIterations){ + if(instance==NULL) return; + PyObject *retVal = PyObject_CallMethod(instance, "run", "i",_iNrIterations); + if(retVal==NULL) return; + Py_DECREF(retVal); +} + +const char ps = +#ifdef _WIN32 + '\\'; +#else + '/'; +#endif + +std::vector CPluginAlgorithmFactory::getPluginPathList(){ + std::vector list; + list.push_back("/etc/astra-toolbox"); + PyObject *ret, *retb; + ret = PyObject_CallMethod(inspect,"getfile","O",astra); + if(ret!=NULL){ + retb = PyObject_CallMethod(six,"b","O",ret); + Py_DECREF(ret); + if(retb!=NULL){ + std::string astra_inst (PyBytes_AsString(retb)); + Py_DECREF(retb); + ret = PyObject_CallMethod(ospath,"dirname","s",astra_inst.c_str()); + if(ret!=NULL){ + retb = PyObject_CallMethod(six,"b","O",ret); + Py_DECREF(ret); + if(retb!=NULL){ + list.push_back(std::string(PyBytes_AsString(retb))); + Py_DECREF(retb); + } + } + } + } + ret = PyObject_CallMethod(ospath,"expanduser","s","~"); + if(ret!=NULL){ + retb = PyObject_CallMethod(six,"b","O",ret); + Py_DECREF(ret); + if(retb!=NULL){ + list.push_back(std::string(PyBytes_AsString(retb)) + ps + ".astra-toolbox"); + Py_DECREF(retb); + } + } + const char *envval = getenv("ASTRA_PLUGIN_PATH"); + if(envval!=NULL){ + list.push_back(std::string(envval)); + } + return list; +} + +CPluginAlgorithmFactory::CPluginAlgorithmFactory(){ + Py_Initialize(); + pluginDict = PyDict_New(); + ospath = PyImport_ImportModule("os.path"); + inspect = PyImport_ImportModule("inspect"); + six = PyImport_ImportModule("six"); + astra = PyImport_ImportModule("astra"); + std::vector fls = getPluginPathList(); + std::vector items; + for(unsigned int i=0;i items; + boost::split(items, str, boost::is_any_of(".")); + PyObject *pyclass = PyImport_ImportModule(items[0].c_str()); + if(pyclass==NULL) return NULL; + PyObject *submod = pyclass; + for(unsigned int i=1;i= 3 +PyObject * pyStringFromString(std::string str){ + return PyUnicode_FromString(str.c_str()); +} +#else +PyObject * pyStringFromString(std::string str){ + return PyBytes_FromString(str.c_str()); +} +#endif + +PyObject* stringToPythonValue(std::string str){ + if(str.find(";")!=std::string::npos){ + std::vector rows, row; + boost::split(rows, str, boost::is_any_of(";")); + PyObject *mat = PyList_New(rows.size()); + for(unsigned int i=0; i(row[j]))); + } + PyList_SetItem(mat, i, rowlist); + } + return mat; + } + if(str.find(",")!=std::string::npos){ + std::vector vec; + boost::split(vec, str, boost::is_any_of(",")); + PyObject *veclist = PyList_New(vec.size()); + for(unsigned int i=0;i(vec[i]))); + } + return veclist; + } + try{ + return PyLong_FromLong(boost::lexical_cast(str)); + }catch(const boost::bad_lexical_cast &){ + try{ + return PyFloat_FromDouble(boost::lexical_cast(str)); + }catch(const boost::bad_lexical_cast &){ + return pyStringFromString(str); + } + } +} + +PyObject* XMLNode2dict(XMLNode node){ + PyObject *dct = PyDict_New(); + PyObject *opts = PyDict_New(); + if(node.hasAttribute("type")){ + PyObject *obj = pyStringFromString(node.getAttribute("type").c_str()); + PyDict_SetItemString(dct, "type", obj); + Py_DECREF(obj); + } + std::list nodes = node.getNodes(); + std::list::iterator it = nodes.begin(); + while(it!=nodes.end()){ + XMLNode subnode = *it; + if(subnode.getName()=="Option"){ + PyObject *obj = stringToPythonValue(subnode.getAttribute("value")); + PyDict_SetItemString(opts, subnode.getAttribute("key").c_str(), obj); + Py_DECREF(obj); + }else{ + PyObject *obj = stringToPythonValue(subnode.getContent()); + PyDict_SetItemString(dct, subnode.getName().c_str(), obj); + Py_DECREF(obj); + } + ++it; + } + PyDict_SetItemString(dct, "options", opts); + Py_DECREF(opts); + return dct; +} + +} +#endif \ No newline at end of file -- cgit v1.2.3 From 4c9e432ae4581fdc110e9a9c45267227be1c7c31 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Wed, 24 Jun 2015 20:43:05 +0200 Subject: Fix config to dict translation for array options --- src/PluginAlgorithm.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index df13f31..a27ce2c 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -275,7 +275,12 @@ PyObject* XMLNode2dict(XMLNode node){ while(it!=nodes.end()){ XMLNode subnode = *it; if(subnode.getName()=="Option"){ - PyObject *obj = stringToPythonValue(subnode.getAttribute("value")); + PyObject *obj; + if(subnode.hasAttribute("value")){ + obj = stringToPythonValue(subnode.getAttribute("value")); + }else{ + obj = stringToPythonValue(subnode.getContent()); + } PyDict_SetItemString(opts, subnode.getAttribute("key").c_str(), obj); Py_DECREF(obj); }else{ -- cgit v1.2.3 From edae78481cf0e9cbffe335de1e541821758c5da1 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Wed, 24 Jun 2015 21:36:04 +0200 Subject: Log error when running Python plugin algorithm --- src/PluginAlgorithm.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index a27ce2c..7dcaf68 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -62,7 +62,7 @@ bool CPluginAlgorithm::initialize(const Config& _cfg){ void CPluginAlgorithm::run(int _iNrIterations){ if(instance==NULL) return; - PyObject *retVal = PyObject_CallMethod(instance, "run", "i",_iNrIterations); + PyObject *retVal = PyObject_CallMethod(instance, "astra_run", "i",_iNrIterations); if(retVal==NULL) return; Py_DECREF(retVal); } -- cgit v1.2.3 From 2f871bc7068d6c87a7d950ae044ba66b0b8dcd3f Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Fri, 17 Jul 2015 12:05:46 +0200 Subject: Remove config text file loading for plugins --- src/PluginAlgorithm.cpp | 72 ++++--------------------------------------------- 1 file changed, 5 insertions(+), 67 deletions(-) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index 7dcaf68..8ba6631 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -67,79 +67,19 @@ void CPluginAlgorithm::run(int _iNrIterations){ Py_DECREF(retVal); } -const char ps = -#ifdef _WIN32 - '\\'; -#else - '/'; -#endif - -std::vector CPluginAlgorithmFactory::getPluginPathList(){ - std::vector list; - list.push_back("/etc/astra-toolbox"); - PyObject *ret, *retb; - ret = PyObject_CallMethod(inspect,"getfile","O",astra); - if(ret!=NULL){ - retb = PyObject_CallMethod(six,"b","O",ret); - Py_DECREF(ret); - if(retb!=NULL){ - std::string astra_inst (PyBytes_AsString(retb)); - Py_DECREF(retb); - ret = PyObject_CallMethod(ospath,"dirname","s",astra_inst.c_str()); - if(ret!=NULL){ - retb = PyObject_CallMethod(six,"b","O",ret); - Py_DECREF(ret); - if(retb!=NULL){ - list.push_back(std::string(PyBytes_AsString(retb))); - Py_DECREF(retb); - } - } - } - } - ret = PyObject_CallMethod(ospath,"expanduser","s","~"); - if(ret!=NULL){ - retb = PyObject_CallMethod(six,"b","O",ret); - Py_DECREF(ret); - if(retb!=NULL){ - list.push_back(std::string(PyBytes_AsString(retb)) + ps + ".astra-toolbox"); - Py_DECREF(retb); - } - } - const char *envval = getenv("ASTRA_PLUGIN_PATH"); - if(envval!=NULL){ - list.push_back(std::string(envval)); - } - return list; -} - CPluginAlgorithmFactory::CPluginAlgorithmFactory(){ Py_Initialize(); pluginDict = PyDict_New(); - ospath = PyImport_ImportModule("os.path"); inspect = PyImport_ImportModule("inspect"); six = PyImport_ImportModule("six"); - astra = PyImport_ImportModule("astra"); - std::vector fls = getPluginPathList(); - std::vector items; - for(unsigned int i=0;i Date: Fri, 17 Jul 2015 16:22:05 +0200 Subject: Fix numpy lapack loading when running in Matlab --- src/Globals.cpp | 3 +++ src/PluginAlgorithm.cpp | 29 +++++++++++++++++++++++++++++ 2 files changed, 32 insertions(+) (limited to 'src') diff --git a/src/Globals.cpp b/src/Globals.cpp index 813f9c9..904a459 100644 --- a/src/Globals.cpp +++ b/src/Globals.cpp @@ -28,5 +28,8 @@ $Id$ #include "astra/Globals.h" +namespace astra{ + bool running_in_matlab=false; +} // nothing to see here :) diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index 8ba6631..c26ee3f 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -67,8 +67,37 @@ void CPluginAlgorithm::run(int _iNrIterations){ Py_DECREF(retVal); } +void fixLapackLoading(){ + // When running in Matlab, we need to force numpy + // to use its internal lapack library instead of + // Matlab's MKL library to avoid errors. To do this, + // we set Python's dlopen flags to RTLD_NOW|RTLD_DEEPBIND + // and import 'numpy.linalg.lapack_lite' here. We reset + // Python's dlopen flags afterwards. + PyObject *sys = PyImport_ImportModule("sys"); + if(sys!=NULL){ + PyObject *curFlags = PyObject_CallMethod(sys,"getdlopenflags",NULL); + if(curFlags!=NULL){ + PyObject *retVal = PyObject_CallMethod(sys, "setdlopenflags", "i",10); + if(retVal!=NULL){ + PyObject *lapack = PyImport_ImportModule("numpy.linalg.lapack_lite"); + if(lapack!=NULL){ + Py_DECREF(lapack); + } + PyObject_CallMethod(sys, "setdlopenflags", "O",curFlags); + Py_DECREF(retVal); + } + Py_DECREF(curFlags); + } + Py_DECREF(sys); + } +} + CPluginAlgorithmFactory::CPluginAlgorithmFactory(){ Py_Initialize(); +#ifndef _MSC_VER + if(astra::running_in_matlab) fixLapackLoading(); +#endif pluginDict = PyDict_New(); inspect = PyImport_ImportModule("inspect"); six = PyImport_ImportModule("six"); -- cgit v1.2.3 From ef9eb1dc7eb494e87f728af7caff8e5291cf320c Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Mon, 20 Jul 2015 10:34:55 +0200 Subject: Also log Python errors when importing and creating Python plugins --- src/PluginAlgorithm.cpp | 48 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index c26ee3f..a118f54 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -29,6 +29,7 @@ $Id$ #ifdef ASTRA_PYTHON #include "astra/PluginAlgorithm.h" +#include "astra/Logging.h" #include #include #include @@ -38,8 +39,53 @@ $Id$ namespace astra { + +void logPythonError(){ + if(PyErr_Occurred()){ + PyObject *ptype, *pvalue, *ptraceback; + PyErr_Fetch(&ptype, &pvalue, &ptraceback); + PyObject *traceback = PyImport_ImportModule("traceback"); + if(traceback!=NULL){ + PyObject *exc; + if(ptraceback==NULL){ + exc = PyObject_CallMethod(traceback,"format_exception_only","OO",ptype, pvalue); + }else{ + exc = PyObject_CallMethod(traceback,"format_exception","OOO",ptype, pvalue, ptraceback); + } + if(exc!=NULL){ + PyObject *six = PyImport_ImportModule("six"); + if(six!=NULL){ + PyObject *iter = PyObject_GetIter(exc); + if(iter!=NULL){ + PyObject *line; + std::string errStr = ""; + while(line = PyIter_Next(iter)){ + PyObject *retb = PyObject_CallMethod(six,"b","O",line); + if(retb!=NULL){ + errStr += std::string(PyBytes_AsString(retb)); + Py_DECREF(retb); + } + Py_DECREF(line); + } + ASTRA_ERROR("%s",errStr.c_str()); + Py_DECREF(iter); + } + Py_DECREF(six); + } + Py_DECREF(exc); + } + Py_DECREF(traceback); + } + if(ptype!=NULL) Py_DECREF(ptype); + if(pvalue!=NULL) Py_DECREF(pvalue); + if(ptraceback!=NULL) Py_DECREF(ptraceback); + } +} + + CPluginAlgorithm::CPluginAlgorithm(PyObject* pyclass){ instance = PyObject_CallObject(pyclass, NULL); + if(instance==NULL) logPythonError(); } CPluginAlgorithm::~CPluginAlgorithm(){ @@ -148,6 +194,8 @@ CPluginAlgorithm * CPluginAlgorithmFactory::getPlugin(std::string name){ if(pyclass!=NULL){ alg = new CPluginAlgorithm(pyclass); Py_DECREF(pyclass); + }else{ + logPythonError(); } }else{ alg = new CPluginAlgorithm(className); -- cgit v1.2.3 From 37abc22cf8d26fa3f7e282a1ee50a2a129d5a295 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Mon, 20 Jul 2015 11:26:39 +0200 Subject: Always log Python errors when importing/creating plugins --- src/PluginAlgorithm.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index a118f54..d6cf731 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -173,13 +173,19 @@ PyObject * getClassFromString(std::string str){ std::vector items; boost::split(items, str, boost::is_any_of(".")); PyObject *pyclass = PyImport_ImportModule(items[0].c_str()); - if(pyclass==NULL) return NULL; + if(pyclass==NULL){ + logPythonError(); + return NULL; + } PyObject *submod = pyclass; for(unsigned int i=1;i Date: Mon, 20 Jul 2015 14:07:21 +0200 Subject: Allow registering plugins without explicit name, and fix exception handling when running in Matlab --- src/PluginAlgorithm.cpp | 95 +++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 76 insertions(+), 19 deletions(-) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index d6cf731..7f7ff61 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -100,7 +100,10 @@ bool CPluginAlgorithm::initialize(const Config& _cfg){ PyObject *cfgDict = XMLNode2dict(_cfg.self); PyObject *retVal = PyObject_CallMethod(instance, "astra_init", "O",cfgDict); Py_DECREF(cfgDict); - if(retVal==NULL) return false; + if(retVal==NULL){ + logPythonError(); + return false; + } m_bIsInitialized = true; Py_DECREF(retVal); return m_bIsInitialized; @@ -108,8 +111,11 @@ bool CPluginAlgorithm::initialize(const Config& _cfg){ void CPluginAlgorithm::run(int _iNrIterations){ if(instance==NULL) return; - PyObject *retVal = PyObject_CallMethod(instance, "astra_run", "i",_iNrIterations); - if(retVal==NULL) return; + PyObject *retVal = PyObject_CallMethod(instance, "run", "i",_iNrIterations); + if(retVal==NULL){ + logPythonError(); + return; + } Py_DECREF(retVal); } @@ -157,18 +163,6 @@ CPluginAlgorithmFactory::~CPluginAlgorithmFactory(){ if(six!=NULL) Py_DECREF(six); } -bool CPluginAlgorithmFactory::registerPlugin(std::string name, std::string className){ - PyObject *str = PyBytes_FromString(className.c_str()); - PyDict_SetItemString(pluginDict, name.c_str(), str); - Py_DECREF(str); - return true; -} - -bool CPluginAlgorithmFactory::registerPluginClass(std::string name, PyObject * className){ - PyDict_SetItemString(pluginDict, name.c_str(), className); - return true; -} - PyObject * getClassFromString(std::string str){ std::vector items; boost::split(items, str, boost::is_any_of(".")); @@ -190,6 +184,43 @@ PyObject * getClassFromString(std::string str){ return pyclass; } +bool CPluginAlgorithmFactory::registerPlugin(std::string name, std::string className){ + PyObject *str = PyBytes_FromString(className.c_str()); + PyDict_SetItemString(pluginDict, name.c_str(), str); + Py_DECREF(str); + return true; +} + +bool CPluginAlgorithmFactory::registerPlugin(std::string className){ + PyObject *pyclass = getClassFromString(className); + if(pyclass==NULL) return false; + bool ret = registerPluginClass(pyclass); + Py_DECREF(pyclass); + return ret; +} + +bool CPluginAlgorithmFactory::registerPluginClass(std::string name, PyObject * className){ + PyDict_SetItemString(pluginDict, name.c_str(), className); + return true; +} + +bool CPluginAlgorithmFactory::registerPluginClass(PyObject * className){ + PyObject *astra_name = PyObject_GetAttrString(className,"astra_name"); + if(astra_name==NULL){ + logPythonError(); + return false; + } + PyObject *retb = PyObject_CallMethod(six,"b","O",astra_name); + if(retb!=NULL){ + PyDict_SetItemString(pluginDict,PyBytes_AsString(retb),className); + Py_DECREF(retb); + }else{ + logPythonError(); + } + Py_DECREF(astra_name); + return true; +} + CPluginAlgorithm * CPluginAlgorithmFactory::getPlugin(std::string name){ PyObject *className = PyDict_GetItemString(pluginDict, name.c_str()); if(className==NULL) return NULL; @@ -212,12 +243,34 @@ PyObject * CPluginAlgorithmFactory::getRegistered(){ return pluginDict; } +std::map CPluginAlgorithmFactory::getRegisteredMap(){ + std::map ret; + PyObject *key, *value; + Py_ssize_t pos = 0; + while (PyDict_Next(pluginDict, &pos, &key, &value)) { + PyObject * keyb = PyObject_Bytes(key); + PyObject * valb = PyObject_Bytes(value); + ret[PyBytes_AsString(keyb)] = PyBytes_AsString(valb); + Py_DECREF(keyb); + Py_DECREF(valb); + } + return ret; +} + std::string CPluginAlgorithmFactory::getHelp(std::string name){ PyObject *className = PyDict_GetItemString(pluginDict, name.c_str()); - if(className==NULL) return ""; - std::string str = std::string(PyBytes_AsString(className)); + if(className==NULL){ + ASTRA_ERROR("Plugin %s not found!",name.c_str()); + return ""; + } std::string ret = ""; - PyObject *pyclass = getClassFromString(str); + PyObject *pyclass; + if(PyBytes_Check(className)){ + std::string str = std::string(PyBytes_AsString(className)); + pyclass = getClassFromString(str); + }else{ + pyclass = className; + } if(pyclass==NULL) return ""; if(inspect!=NULL && six!=NULL){ PyObject *retVal = PyObject_CallMethod(inspect,"getdoc","O",pyclass); @@ -228,9 +281,13 @@ std::string CPluginAlgorithmFactory::getHelp(std::string name){ ret = std::string(PyBytes_AsString(retb)); Py_DECREF(retb); } + }else{ + logPythonError(); } } - Py_DECREF(pyclass); + if(PyBytes_Check(className)){ + Py_DECREF(pyclass); + } return ret; } -- cgit v1.2.3 From 3808967cfaa6beb9d93d2035ebc72fa010fdab11 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Mon, 20 Jul 2015 16:41:55 +0200 Subject: Normalize Python exceptions (needed for some) --- src/PluginAlgorithm.cpp | 1 + 1 file changed, 1 insertion(+) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index 7f7ff61..56c4e4d 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -44,6 +44,7 @@ void logPythonError(){ if(PyErr_Occurred()){ PyObject *ptype, *pvalue, *ptraceback; PyErr_Fetch(&ptype, &pvalue, &ptraceback); + PyErr_NormalizeException(&ptype, &pvalue, &ptraceback); PyObject *traceback = PyImport_ImportModule("traceback"); if(traceback!=NULL){ PyObject *exc; -- cgit v1.2.3 From dc3bed557603d4735ddc20961c28e5e868fc315c Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Tue, 21 Jul 2015 11:44:29 +0200 Subject: Clear Python error when plugin is not find in getHelp --- src/PluginAlgorithm.cpp | 1 + 1 file changed, 1 insertion(+) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index 56c4e4d..5c779fd 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -262,6 +262,7 @@ std::string CPluginAlgorithmFactory::getHelp(std::string name){ PyObject *className = PyDict_GetItemString(pluginDict, name.c_str()); if(className==NULL){ ASTRA_ERROR("Plugin %s not found!",name.c_str()); + PyErr_Clear(); return ""; } std::string ret = ""; -- cgit v1.2.3 From 645122f4b365ce44849afda2ed8a711ae649ed76 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Tue, 21 Jul 2015 13:56:18 +0200 Subject: Fix 'get_registered' in Matlab with Python 3 --- src/PluginAlgorithm.cpp | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index 5c779fd..5d6d733 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -249,11 +249,24 @@ std::map CPluginAlgorithmFactory::getRegisteredMap(){ PyObject *key, *value; Py_ssize_t pos = 0; while (PyDict_Next(pluginDict, &pos, &key, &value)) { - PyObject * keyb = PyObject_Bytes(key); - PyObject * valb = PyObject_Bytes(value); - ret[PyBytes_AsString(keyb)] = PyBytes_AsString(valb); - Py_DECREF(keyb); - Py_DECREF(valb); + PyObject *keystr = PyObject_Str(key); + if(keystr!=NULL){ + PyObject *valstr = PyObject_Str(value); + if(valstr!=NULL){ + PyObject * keyb = PyObject_CallMethod(six,"b","O",keystr); + if(keyb!=NULL){ + PyObject * valb = PyObject_CallMethod(six,"b","O",valstr); + if(valb!=NULL){ + ret[PyBytes_AsString(keyb)] = PyBytes_AsString(valb); + Py_DECREF(valb); + } + Py_DECREF(keyb); + } + Py_DECREF(valstr); + } + Py_DECREF(keystr); + } + logPythonError(); } return ret; } -- cgit v1.2.3 From ab980d9f088c0f4e28d61b94c32788c30a9c4cb9 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Wed, 5 Aug 2015 16:26:01 +0200 Subject: Fix get_help for classes without docstring --- src/PluginAlgorithm.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index 5d6d733..4066e30 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -290,12 +290,14 @@ std::string CPluginAlgorithmFactory::getHelp(std::string name){ if(inspect!=NULL && six!=NULL){ PyObject *retVal = PyObject_CallMethod(inspect,"getdoc","O",pyclass); if(retVal!=NULL){ - PyObject *retb = PyObject_CallMethod(six,"b","O",retVal); - Py_DECREF(retVal); - if(retb!=NULL){ - ret = std::string(PyBytes_AsString(retb)); - Py_DECREF(retb); + if(retVal!=Py_None){ + PyObject *retb = PyObject_CallMethod(six,"b","O",retVal); + if(retb!=NULL){ + ret = std::string(PyBytes_AsString(retb)); + Py_DECREF(retb); + } } + Py_DECREF(retVal); }else{ logPythonError(); } -- cgit v1.2.3 From 0d5947a0e8e7d6f86c7591a96d877dfe14b187e4 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Mon, 10 Aug 2015 17:08:34 +0200 Subject: Ensure we have acquired the GIL before calling Python plugin 'run' method --- src/PluginAlgorithm.cpp | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index 4066e30..e79c77b 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -112,12 +112,14 @@ bool CPluginAlgorithm::initialize(const Config& _cfg){ void CPluginAlgorithm::run(int _iNrIterations){ if(instance==NULL) return; + PyGILState_STATE state = PyGILState_Ensure(); PyObject *retVal = PyObject_CallMethod(instance, "run", "i",_iNrIterations); if(retVal==NULL){ logPythonError(); - return; + }else{ + Py_DECREF(retVal); } - Py_DECREF(retVal); + PyGILState_Release(state); } void fixLapackLoading(){ @@ -147,7 +149,10 @@ void fixLapackLoading(){ } CPluginAlgorithmFactory::CPluginAlgorithmFactory(){ - Py_Initialize(); + if(!Py_IsInitialized()){ + Py_Initialize(); + PyEval_InitThreads(); + } #ifndef _MSC_VER if(astra::running_in_matlab) fixLapackLoading(); #endif -- cgit v1.2.3 From 43a38c117405f99e3a1b498f899de4ba6d01a044 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 7 Oct 2015 18:14:39 +0200 Subject: Improve option passing through CudaProjector3D Not all constructors were reading options from the projector. Also allow passing GPUIndex via CudaProjector3D. Thanks to Nicola Vigano for part of the patch. --- src/CudaBackProjectionAlgorithm3D.cpp | 43 ++++++++++++++++++++-------- src/CudaCglsAlgorithm3D.cpp | 47 +++++++++++++++++++++---------- src/CudaFDKAlgorithm3D.cpp | 40 ++++++++++++++++++-------- src/CudaForwardProjectionAlgorithm3D.cpp | 46 ++++++++++++++++++++---------- src/CudaProjector3D.cpp | 7 +++++ src/CudaSirtAlgorithm3D.cpp | 48 +++++++++++++++++++++----------- 6 files changed, 161 insertions(+), 70 deletions(-) (limited to 'src') diff --git a/src/CudaBackProjectionAlgorithm3D.cpp b/src/CudaBackProjectionAlgorithm3D.cpp index e8e0433..c9d9447 100644 --- a/src/CudaBackProjectionAlgorithm3D.cpp +++ b/src/CudaBackProjectionAlgorithm3D.cpp @@ -38,6 +38,8 @@ $Id$ #include "astra/ParallelVecProjectionGeometry3D.h" #include "astra/ConeVecProjectionGeometry3D.h" +#include "astra/Logging.h" + #include "../cuda/3d/astra3d.h" using namespace std; @@ -86,6 +88,24 @@ bool CCudaBackProjectionAlgorithm3D::_check() return true; } +//--------------------------------------------------------------------------------------- +void CCudaBackProjectionAlgorithm3D::initializeFromProjector() +{ + m_iVoxelSuperSampling = 1; + m_iGPUIndex = -1; + + CCudaProjector3D* pCudaProjector = dynamic_cast(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector3D passed to BP3D_CUDA"); + } + } else { + m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } + +} + //--------------------------------------------------------------------------------------- // Initialize - Config bool CCudaBackProjectionAlgorithm3D::initialize(const Config& _cfg) @@ -103,21 +123,18 @@ bool CCudaBackProjectionAlgorithm3D::initialize(const Config& _cfg) return false; } - CCudaProjector3D* pCudaProjector = 0; - pCudaProjector = dynamic_cast(m_pProjector); - if (!pCudaProjector) { - // TODO: Report - } - - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - CC.markOptionParsed("GPUindex"); - + initializeFromProjector(); - m_iVoxelSuperSampling = 1; - if (pCudaProjector) - m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + // Deprecated options m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); CC.markOptionParsed("VoxelSuperSampling"); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); + + CFloat32ProjectionData3DMemory* pSinoMem = dynamic_cast(m_pSinogram); ASTRA_ASSERT(pSinoMem); @@ -151,6 +168,8 @@ bool CCudaBackProjectionAlgorithm3D::initialize(CProjector3D* _pProjector, m_pSinogram = _pSinogram; m_pReconstruction = _pReconstruction; + initializeFromProjector(); + // success m_bIsInitialized = _check(); return m_bIsInitialized; diff --git a/src/CudaCglsAlgorithm3D.cpp b/src/CudaCglsAlgorithm3D.cpp index f527dc5..1cccb6a 100644 --- a/src/CudaCglsAlgorithm3D.cpp +++ b/src/CudaCglsAlgorithm3D.cpp @@ -37,6 +37,8 @@ $Id$ #include "astra/ParallelVecProjectionGeometry3D.h" #include "astra/ConeVecProjectionGeometry3D.h" +#include "astra/Logging.h" + #include "../cuda/3d/astra3d.h" using namespace std; @@ -89,6 +91,26 @@ bool CCudaCglsAlgorithm3D::_check() return true; } +//--------------------------------------------------------------------------------------- +void CCudaCglsAlgorithm3D::initializeFromProjector() +{ + m_iVoxelSuperSampling = 1; + m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; + + CCudaProjector3D* pCudaProjector = dynamic_cast(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector3D passed to CGLS3D_CUDA"); + } + } else { + m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } + +} + //--------------------------------------------------------------------------------------- // Initialize - Config bool CCudaCglsAlgorithm3D::initialize(const Config& _cfg) @@ -107,27 +129,20 @@ bool CCudaCglsAlgorithm3D::initialize(const Config& _cfg) return false; } - CCudaProjector3D* pCudaProjector = 0; - pCudaProjector = dynamic_cast(m_pProjector); - if (!pCudaProjector) { - // TODO: Report - } + initializeFromProjector(); - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - CC.markOptionParsed("GPUindex"); - - m_iVoxelSuperSampling = 1; - m_iDetectorSuperSampling = 1; - if (pCudaProjector) { - // New interface - m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); - m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); - } // Deprecated options m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling); 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"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); + + m_pCgls = new AstraCGLS3d(); @@ -155,6 +170,8 @@ bool CCudaCglsAlgorithm3D::initialize(CProjector3D* _pProjector, m_pSinogram = _pSinogram; m_pReconstruction = _pReconstruction; + initializeFromProjector(); + m_pCgls = new AstraCGLS3d; m_bAstraCGLSInit = false; diff --git a/src/CudaFDKAlgorithm3D.cpp b/src/CudaFDKAlgorithm3D.cpp index 667d926..625d02a 100644 --- a/src/CudaFDKAlgorithm3D.cpp +++ b/src/CudaFDKAlgorithm3D.cpp @@ -35,6 +35,8 @@ $Id$ #include "astra/CudaProjector3D.h" #include "astra/ConeProjectionGeometry3D.h" +#include "astra/Logging.h" + #include "../cuda/3d/astra3d.h" using namespace std; @@ -84,6 +86,24 @@ bool CCudaFDKAlgorithm3D::_check() return true; } +//--------------------------------------------------------------------------------------- +void CCudaFDKAlgorithm3D::initializeFromProjector() +{ + m_iVoxelSuperSampling = 1; + m_iGPUIndex = -1; + + CCudaProjector3D* pCudaProjector = dynamic_cast(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector3D passed to FDK_CUDA"); + } + } else { + m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } + +} + //--------------------------------------------------------------------------------------- // Initialize - Config bool CCudaFDKAlgorithm3D::initialize(const Config& _cfg) @@ -101,20 +121,18 @@ bool CCudaFDKAlgorithm3D::initialize(const Config& _cfg) return false; } - CCudaProjector3D* pCudaProjector = 0; - pCudaProjector = dynamic_cast(m_pProjector); - if (!pCudaProjector) { - // TODO: Report - } - - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - CC.markOptionParsed("GPUindex"); + initializeFromProjector(); - m_iVoxelSuperSampling = 1; - if (pCudaProjector) - m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + // Deprecated options m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); CC.markOptionParsed("VoxelSuperSampling"); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); + + m_bShortScan = _cfg.self.getOptionBool("ShortScan", false); CC.markOptionParsed("ShortScan"); diff --git a/src/CudaForwardProjectionAlgorithm3D.cpp b/src/CudaForwardProjectionAlgorithm3D.cpp index 46dab12..6498885 100644 --- a/src/CudaForwardProjectionAlgorithm3D.cpp +++ b/src/CudaForwardProjectionAlgorithm3D.cpp @@ -71,6 +71,23 @@ CCudaForwardProjectionAlgorithm3D::~CCudaForwardProjectionAlgorithm3D() } +//--------------------------------------------------------------------------------------- +void CCudaForwardProjectionAlgorithm3D::initializeFromProjector() +{ + m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; + + CCudaProjector3D* pCudaProjector = dynamic_cast(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector3D passed to FP3D_CUDA"); + } + } else { + m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } +} + //--------------------------------------------------------------------------------------- // Initialize - Config bool CCudaForwardProjectionAlgorithm3D::initialize(const Config& _cfg) @@ -97,29 +114,21 @@ bool CCudaForwardProjectionAlgorithm3D::initialize(const Config& _cfg) // optional: projector node = _cfg.self.getSingleNode("ProjectorId"); - CCudaProjector3D* pCudaProjector = 0; m_pProjector = 0; if (node) { id = boost::lexical_cast(node.getContent()); m_pProjector = CProjector3DManager::getSingleton().get(id); - pCudaProjector = dynamic_cast(CProjector3DManager::getSingleton().get(id)); - m_pProjector = pCudaProjector; - if (!pCudaProjector) { - // TODO: Report - } } CC.markNodeParsed("ProjectorId"); - // GPU number - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - CC.markOptionParsed("GPUindex"); - + initializeFromProjector(); - m_iDetectorSuperSampling = 1; - if (pCudaProjector) - m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); + // Deprecated options m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex); CC.markOptionParsed("DetectorSuperSampling"); + CC.markOptionParsed("GPUindex"); + // success m_bIsInitialized = check(); @@ -142,8 +151,15 @@ bool CCudaForwardProjectionAlgorithm3D::initialize(CProjector3D* _pProjector, m_pProjections = _pProjections; m_pVolume = _pVolume; - m_iDetectorSuperSampling = _iDetectorSuperSampling; - m_iGPUIndex = _iGPUindex; + CCudaProjector3D* pCudaProjector = dynamic_cast(m_pProjector); + if (!pCudaProjector) { + // TODO: Report + m_iDetectorSuperSampling = _iDetectorSuperSampling; + m_iGPUIndex = _iGPUindex; + } else { + m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } // success m_bIsInitialized = check(); diff --git a/src/CudaProjector3D.cpp b/src/CudaProjector3D.cpp index d2fd74c..bbfbd34 100644 --- a/src/CudaProjector3D.cpp +++ b/src/CudaProjector3D.cpp @@ -64,6 +64,7 @@ void CCudaProjector3D::_clear() m_projectionKernel = ker3d_default; m_iVoxelSuperSampling = 1; m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; } //---------------------------------------------------------------------------------------- @@ -128,6 +129,12 @@ bool CCudaProjector3D::initialize(const Config& _cfg) m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", 1); CC.markOptionParsed("DetectorSuperSampling"); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); + m_bIsInitialized = _check(); return m_bIsInitialized; } diff --git a/src/CudaSirtAlgorithm3D.cpp b/src/CudaSirtAlgorithm3D.cpp index abbb9fd..67594f4 100644 --- a/src/CudaSirtAlgorithm3D.cpp +++ b/src/CudaSirtAlgorithm3D.cpp @@ -38,6 +38,8 @@ $Id$ #include "astra/ConeVecProjectionGeometry3D.h" #include "astra/CudaProjector3D.h" +#include "astra/Logging.h" + #include "../cuda/3d/astra3d.h" using namespace std; @@ -90,7 +92,27 @@ bool CCudaSirtAlgorithm3D::_check() return true; } -//--------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------------- +void CCudaSirtAlgorithm3D::initializeFromProjector() +{ + m_iVoxelSuperSampling = 1; + m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; + + CCudaProjector3D* pCudaProjector = dynamic_cast(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector3D passed to SIRT3D_CUDA"); + } + } else { + m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } + +} + +//-------------------------------------------------------------------------------------- // Initialize - Config bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg) { @@ -108,28 +130,20 @@ bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg) return false; } - CCudaProjector3D* pCudaProjector = 0; - pCudaProjector = dynamic_cast(m_pProjector); - if (!pCudaProjector) { - // TODO: Report - } + initializeFromProjector(); - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - CC.markOptionParsed("GPUindex"); - - - m_iVoxelSuperSampling = 1; - m_iDetectorSuperSampling = 1; - if (pCudaProjector) { - // New interface - m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling(); - m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); - } // Deprecated options m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling); 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"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); + + m_pSirt = new AstraSIRT3d(); -- cgit v1.2.3 From 003663649a191fc5bc011d6e5424496576b5e793 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 8 Oct 2015 11:24:49 +0200 Subject: Improve option passing through CudaProjector2D Not all constructors were reading options from the projector. Also allow passing GPUIndex via CudaProjector2D. Also refactor CudaReconstructionAlgorithm::initialize/check to avoid code duplication with ReconstructionAlgorithm. --- src/CudaBackProjectionAlgorithm.cpp | 5 +- src/CudaCglsAlgorithm.cpp | 6 +- src/CudaEMAlgorithm.cpp | 6 +- src/CudaFilteredBackProjectionAlgorithm.cpp | 43 ++++--- src/CudaForwardProjectionAlgorithm.cpp | 60 +++++----- src/CudaProjector2D.cpp | 17 +-- src/CudaReconstructionAlgorithm2D.cpp | 169 ++++++---------------------- src/CudaSartAlgorithm.cpp | 5 +- src/CudaSirtAlgorithm.cpp | 6 +- src/ReconstructionAlgorithm2D.cpp | 25 ++-- 10 files changed, 134 insertions(+), 208 deletions(-) (limited to 'src') diff --git a/src/CudaBackProjectionAlgorithm.cpp b/src/CudaBackProjectionAlgorithm.cpp index 365e058..a73f895 100644 --- a/src/CudaBackProjectionAlgorithm.cpp +++ b/src/CudaBackProjectionAlgorithm.cpp @@ -76,10 +76,9 @@ bool CCudaBackProjectionAlgorithm::initialize(const Config& _cfg) // Initialize - C++ bool CCudaBackProjectionAlgorithm::initialize(CProjector2D* _pProjector, CFloat32ProjectionData2D* _pSinogram, - CFloat32VolumeData2D* _pReconstruction, - int _iGPUindex, int _iPixelSuperSampling) + CFloat32VolumeData2D* _pReconstruction) { - m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, 1, _iPixelSuperSampling); + m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction); if (!m_bIsInitialized) return false; diff --git a/src/CudaCglsAlgorithm.cpp b/src/CudaCglsAlgorithm.cpp index 0cedff6..9dd4f78 100644 --- a/src/CudaCglsAlgorithm.cpp +++ b/src/CudaCglsAlgorithm.cpp @@ -77,11 +77,9 @@ bool CCudaCglsAlgorithm::initialize(const Config& _cfg) // Initialize - C++ bool CCudaCglsAlgorithm::initialize(CProjector2D* _pProjector, CFloat32ProjectionData2D* _pSinogram, - CFloat32VolumeData2D* _pReconstruction, - int _iGPUindex, int _iDetectorSuperSampling, - int _iPixelSuperSampling) + CFloat32VolumeData2D* _pReconstruction) { - m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, _iDetectorSuperSampling, _iPixelSuperSampling); + m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction); if (!m_bIsInitialized) return false; diff --git a/src/CudaEMAlgorithm.cpp b/src/CudaEMAlgorithm.cpp index 5c71f3d..d0afd80 100644 --- a/src/CudaEMAlgorithm.cpp +++ b/src/CudaEMAlgorithm.cpp @@ -76,11 +76,9 @@ bool CCudaEMAlgorithm::initialize(const Config& _cfg) // Initialize - C++ bool CCudaEMAlgorithm::initialize(CProjector2D* _pProjector, CFloat32ProjectionData2D* _pSinogram, - CFloat32VolumeData2D* _pReconstruction, - int _iGPUindex, int _iDetectorSuperSampling, - int _iPixelSuperSampling) + CFloat32VolumeData2D* _pReconstruction) { - m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, _iDetectorSuperSampling, _iPixelSuperSampling); + m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction); if (!m_bIsInitialized) return false; diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp index aac96d6..8c0659d 100644 --- a/src/CudaFilteredBackProjectionAlgorithm.cpp +++ b/src/CudaFilteredBackProjectionAlgorithm.cpp @@ -67,6 +67,24 @@ CCudaFilteredBackProjectionAlgorithm::~CCudaFilteredBackProjectionAlgorithm() } } +void CCudaFilteredBackProjectionAlgorithm::initializeFromProjector() +{ + m_iPixelSuperSampling = 1; + m_iGPUIndex = -1; + + // Projector + CCudaProjector2D* pCudaProjector = dynamic_cast(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector2D passed to FBP_CUDA"); + } + } else { + m_iPixelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } + +} + bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) { ASTRA_ASSERT(_cfg.self); @@ -163,27 +181,24 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) } CC.markNodeParsed("FilterD"); // TODO: Only for some types! - // GPU number - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); - CC.markOptionParsed("GPUindex"); - - m_iPixelSuperSampling = 1; - if (pCudaProjector) { - // New interface - m_iPixelSuperSampling = pCudaProjector->getVoxelSuperSampling(); - } - // Deprecated options - m_iPixelSuperSampling = (int)_cfg.self.getOptionNumerical("PixelSuperSampling", m_iPixelSuperSampling); - CC.markOptionParsed("PixelSuperSampling"); - - // Fan beam short scan mode if (m_pSinogram && dynamic_cast(m_pSinogram->getGeometry())) { m_bShortScan = (int)_cfg.self.getOptionBool("ShortScan", false); CC.markOptionParsed("ShortScan"); } + initializeFromProjector(); + // Deprecated options + m_iPixelSuperSampling = (int)_cfg.self.getOptionNumerical("PixelSuperSampling", m_iPixelSuperSampling); + CC.markOptionParsed("PixelSuperSampling"); + + // GPU number + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); m_pFBP = new AstraFBP; diff --git a/src/CudaForwardProjectionAlgorithm.cpp b/src/CudaForwardProjectionAlgorithm.cpp index b382f2e..9ca13ae 100644 --- a/src/CudaForwardProjectionAlgorithm.cpp +++ b/src/CudaForwardProjectionAlgorithm.cpp @@ -38,8 +38,11 @@ $Id$ #include #include "astra/AstraObjectManager.h" +#include "astra/ParallelProjectionGeometry2D.h" #include "astra/FanFlatProjectionGeometry2D.h" #include "astra/FanFlatVecProjectionGeometry2D.h" +#include "astra/Float32ProjectionData2D.h" +#include "astra/Float32VolumeData2D.h" #include "astra/CudaProjector2D.h" #include "astra/Logging.h" @@ -65,6 +68,24 @@ CCudaForwardProjectionAlgorithm::~CCudaForwardProjectionAlgorithm() } +//--------------------------------------------------------------------------------------- +void CCudaForwardProjectionAlgorithm::initializeFromProjector() +{ + m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; + + // Projector + CCudaProjector2D* pCudaProjector = dynamic_cast(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { + ASTRA_WARN("non-CUDA Projector2D passed to FP_CUDA"); + } + } else { + m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); + } +} + //--------------------------------------------------------------------------------------- // Initialize - Config bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg) @@ -74,14 +95,9 @@ bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg) // Projector XMLNode node = _cfg.self.getSingleNode("ProjectorId"); - CCudaProjector2D* pCudaProjector = 0; if (node) { int id = boost::lexical_cast(node.getContent()); - CProjector2D *projector = CProjector2DManager::getSingleton().get(id); - pCudaProjector = dynamic_cast(projector); - if (!pCudaProjector) { - ASTRA_WARN("non-CUDA Projector2D passed to FP_CUDA"); - } + m_pProjector = CProjector2DManager::getSingleton().get(id); } CC.markNodeParsed("ProjectorId"); @@ -101,22 +117,18 @@ bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg) m_pVolume = dynamic_cast(CData2DManager::getSingleton().get(id)); CC.markNodeParsed("VolumeDataId"); + initializeFromProjector(); + + // Deprecated options + m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling); + CC.markOptionParsed("DetectorSuperSampling"); // GPU number m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); - CC.markOptionParsed("GPUindex"); - if (!_cfg.self.hasOption("GPUindex")) - CC.markOptionParsed("GPUIndex"); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); - // Detector supersampling factor - m_iDetectorSuperSampling = 1; - if (pCudaProjector) { - // New interface - m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); - } - // Deprecated option - m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling); - CC.markOptionParsed("DetectorSuperSampling"); // return success @@ -125,20 +137,16 @@ bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg) //---------------------------------------------------------------------------------------- // Initialize - C++ -bool CCudaForwardProjectionAlgorithm::initialize(CProjectionGeometry2D* _pProjectionGeometry, - CVolumeGeometry2D* _pReconstructionGeometry, +bool CCudaForwardProjectionAlgorithm::initialize(CProjector2D* _pProjector, CFloat32VolumeData2D* _pVolume, - CFloat32ProjectionData2D* _pSinogram, - int _iGPUindex, int _iDetectorSuperSampling) + CFloat32ProjectionData2D* _pSinogram) { // store classes - //m_pProjectionGeometry = _pProjectionGeometry; - //m_pReconstructionGeometry = _pReconstructionGeometry; + m_pProjector = _pProjector; m_pVolume = _pVolume; m_pSinogram = _pSinogram; - m_iDetectorSuperSampling = _iDetectorSuperSampling; - m_iGPUIndex = _iGPUindex; + initializeFromProjector(); // return success return check(); diff --git a/src/CudaProjector2D.cpp b/src/CudaProjector2D.cpp index a26e32d..acf6000 100644 --- a/src/CudaProjector2D.cpp +++ b/src/CudaProjector2D.cpp @@ -61,6 +61,7 @@ void CCudaProjector2D::_clear() m_projectionKernel = ker2d_default; m_iVoxelSuperSampling = 1; m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; } //---------------------------------------------------------------------------------------- @@ -125,18 +126,18 @@ bool CCudaProjector2D::initialize(const Config& _cfg) m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", 1); CC.markOptionParsed("DetectorSuperSampling"); + // GPU number + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); + + m_bIsInitialized = _check(); return m_bIsInitialized; } -/* -bool CProjector2D::initialize(astra::CProjectionGeometry2D *, astra::CVolumeGeometry2D *) -{ - ASTRA_ASSERT(false); - - return false; -} -*/ std::string CCudaProjector2D::description() const { diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp index 71b6637..bccdb43 100644 --- a/src/CudaReconstructionAlgorithm2D.cpp +++ b/src/CudaReconstructionAlgorithm2D.cpp @@ -84,111 +84,51 @@ void CCudaReconstructionAlgorithm2D::_clear() } //--------------------------------------------------------------------------------------- -// Initialize - Config -bool CCudaReconstructionAlgorithm2D::initialize(const Config& _cfg) +void CCudaReconstructionAlgorithm2D::initializeFromProjector() { - ASTRA_ASSERT(_cfg.self); - ConfigStackCheck CC("CudaReconstructionAlgorithm2D", this, _cfg); - - // if already initialized, clear first - if (m_bIsInitialized) { - clear(); - } + m_iPixelSuperSampling = 1; + m_iDetectorSuperSampling = 1; + m_iGPUIndex = -1; // Projector - XMLNode node = _cfg.self.getSingleNode("ProjectorId"); - CCudaProjector2D* pCudaProjector = 0; - if (node) { - int id = boost::lexical_cast(node.getContent()); - CProjector2D *projector = CProjector2DManager::getSingleton().get(id); - pCudaProjector = dynamic_cast(projector); - if (!pCudaProjector) { + CCudaProjector2D* pCudaProjector = dynamic_cast(m_pProjector); + if (!pCudaProjector) { + if (m_pProjector) { ASTRA_WARN("non-CUDA Projector2D passed"); } - } - CC.markNodeParsed("ProjectorId"); - - - // sinogram data - node = _cfg.self.getSingleNode("ProjectionDataId"); - ASTRA_CONFIG_CHECK(node, "CudaSirt2", "No ProjectionDataId tag specified."); - int id = boost::lexical_cast(node.getContent()); - m_pSinogram = dynamic_cast(CData2DManager::getSingleton().get(id)); - CC.markNodeParsed("ProjectionDataId"); - - // reconstruction data - node = _cfg.self.getSingleNode("ReconstructionDataId"); - ASTRA_CONFIG_CHECK(node, "CudaSirt2", "No ReconstructionDataId tag specified."); - id = boost::lexical_cast(node.getContent()); - m_pReconstruction = dynamic_cast(CData2DManager::getSingleton().get(id)); - CC.markNodeParsed("ReconstructionDataId"); - - // fixed mask - if (_cfg.self.hasOption("ReconstructionMaskId")) { - m_bUseReconstructionMask = true; - id = boost::lexical_cast(_cfg.self.getOption("ReconstructionMaskId")); - m_pReconstructionMask = dynamic_cast(CData2DManager::getSingleton().get(id)); - ASTRA_CONFIG_CHECK(m_pReconstructionMask, "CudaReconstruction2D", "Invalid ReconstructionMaskId."); - } - CC.markOptionParsed("ReconstructionMaskId"); - // fixed mask - if (_cfg.self.hasOption("SinogramMaskId")) { - m_bUseSinogramMask = true; - id = boost::lexical_cast(_cfg.self.getOption("SinogramMaskId")); - m_pSinogramMask = dynamic_cast(CData2DManager::getSingleton().get(id)); - ASTRA_CONFIG_CHECK(m_pSinogramMask, "CudaReconstruction2D", "Invalid SinogramMaskId."); - } - 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", -1); - m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); - CC.markOptionParsed("GPUindex"); - if (!_cfg.self.hasOption("GPUindex")) - CC.markOptionParsed("GPUIndex"); - - // Supersampling factors - m_iDetectorSuperSampling = 1; - m_iPixelSuperSampling = 1; - if (pCudaProjector) { - // New interface m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling(); m_iPixelSuperSampling = pCudaProjector->getVoxelSuperSampling(); + m_iGPUIndex = pCudaProjector->getGPUIndex(); } +} + +//--------------------------------------------------------------------------------------- +// Initialize - Config +bool CCudaReconstructionAlgorithm2D::initialize(const Config& _cfg) +{ + ASTRA_ASSERT(_cfg.self); + ConfigStackCheck CC("CudaReconstructionAlgorithm2D", this, _cfg); + + m_bIsInitialized = CReconstructionAlgorithm2D::initialize(_cfg); + + if (!m_bIsInitialized) + return false; + + initializeFromProjector(); + // Deprecated options m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling); m_iPixelSuperSampling = (int)_cfg.self.getOptionNumerical("PixelSuperSampling", m_iPixelSuperSampling); CC.markOptionParsed("DetectorSuperSampling"); CC.markOptionParsed("PixelSuperSampling"); + // GPU number + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1); + m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); + CC.markOptionParsed("GPUIndex"); + if (!_cfg.self.hasOption("GPUIndex")) + CC.markOptionParsed("GPUindex"); return _check(); } @@ -198,33 +138,19 @@ bool CCudaReconstructionAlgorithm2D::initialize(const Config& _cfg) 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; + m_pProjector = _pProjector; // required classes m_pSinogram = _pSinogram; m_pReconstruction = _pReconstruction; - m_iDetectorSuperSampling = _iDetectorSuperSampling; - m_iPixelSuperSampling = _iPixelSuperSampling; - m_iGPUIndex = _iGPUindex; + initializeFromProjector(); return _check(); } @@ -234,40 +160,13 @@ bool CCudaReconstructionAlgorithm2D::initialize(CProjector2D* _pProjector, // 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."); + if (!CReconstructionAlgorithm2D::_check()) + return false; 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 >= -1, "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 diff --git a/src/CudaSartAlgorithm.cpp b/src/CudaSartAlgorithm.cpp index 8c0c6d7..d202847 100644 --- a/src/CudaSartAlgorithm.cpp +++ b/src/CudaSartAlgorithm.cpp @@ -116,10 +116,9 @@ bool CCudaSartAlgorithm::initialize(const Config& _cfg) // Initialize - C++ bool CCudaSartAlgorithm::initialize(CProjector2D* _pProjector, CFloat32ProjectionData2D* _pSinogram, - CFloat32VolumeData2D* _pReconstruction, - int _iGPUindex, int _iDetectorSuperSampling) + CFloat32VolumeData2D* _pReconstruction) { - m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, _iDetectorSuperSampling, 1); + m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction); if (!m_bIsInitialized) return false; diff --git a/src/CudaSirtAlgorithm.cpp b/src/CudaSirtAlgorithm.cpp index d424915..ab0a418 100644 --- a/src/CudaSirtAlgorithm.cpp +++ b/src/CudaSirtAlgorithm.cpp @@ -98,11 +98,9 @@ bool CCudaSirtAlgorithm::initialize(const Config& _cfg) // Initialize - C++ bool CCudaSirtAlgorithm::initialize(CProjector2D* _pProjector, CFloat32ProjectionData2D* _pSinogram, - CFloat32VolumeData2D* _pReconstruction, - int _iGPUindex, int _iDetectorSuperSampling, - int _iPixelSuperSampling) + CFloat32VolumeData2D* _pReconstruction) { - m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, _iDetectorSuperSampling, _iPixelSuperSampling); + m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction); if (!m_bIsInitialized) return false; diff --git a/src/ReconstructionAlgorithm2D.cpp b/src/ReconstructionAlgorithm2D.cpp index 767efe6..4575ff7 100644 --- a/src/ReconstructionAlgorithm2D.cpp +++ b/src/ReconstructionAlgorithm2D.cpp @@ -85,9 +85,16 @@ bool CReconstructionAlgorithm2D::initialize(const Config& _cfg) // projector XMLNode node = _cfg.self.getSingleNode("ProjectorId"); - ASTRA_CONFIG_CHECK(node, "Reconstruction2D", "No ProjectorId tag specified."); - int id = boost::lexical_cast(node.getContent()); - m_pProjector = CProjector2DManager::getSingleton().get(id); + if (requiresProjector()) { + ASTRA_CONFIG_CHECK(node, "Reconstruction2D", "No ProjectorId tag specified."); + } + int id; + if (node) { + id = boost::lexical_cast(node.getContent()); + m_pProjector = CProjector2DManager::getSingleton().get(id); + } else { + m_pProjector = 0; + } CC.markNodeParsed("ProjectorId"); // sinogram data @@ -205,18 +212,22 @@ void CReconstructionAlgorithm2D::setSinogramMask(CFloat32ProjectionData2D* _pMas bool CReconstructionAlgorithm2D::_check() { // check pointers - ASTRA_CONFIG_CHECK(m_pProjector, "Reconstruction2D", "Invalid Projector Object."); + if (requiresProjector()) + ASTRA_CONFIG_CHECK(m_pProjector, "Reconstruction2D", "Invalid Projector Object."); ASTRA_CONFIG_CHECK(m_pSinogram, "Reconstruction2D", "Invalid Projection Data Object."); ASTRA_CONFIG_CHECK(m_pReconstruction, "Reconstruction2D", "Invalid Reconstruction Data Object."); // check initializations - ASTRA_CONFIG_CHECK(m_pProjector->isInitialized(), "Reconstruction2D", "Projector Object Not Initialized."); + if (requiresProjector()) + ASTRA_CONFIG_CHECK(m_pProjector->isInitialized(), "Reconstruction2D", "Projector Object Not Initialized."); ASTRA_CONFIG_CHECK(m_pSinogram->isInitialized(), "Reconstruction2D", "Projection Data Object Not Initialized."); ASTRA_CONFIG_CHECK(m_pReconstruction->isInitialized(), "Reconstruction2D", "Reconstruction Data Object Not Initialized."); // check compatibility between projector and data classes - ASTRA_CONFIG_CHECK(m_pSinogram->getGeometry()->isEqual(m_pProjector->getProjectionGeometry()), "Reconstruction2D", "Projection Data not compatible with the specified Projector."); - ASTRA_CONFIG_CHECK(m_pReconstruction->getGeometry()->isEqual(m_pProjector->getVolumeGeometry()), "Reconstruction2D", "Reconstruction Data not compatible with the specified Projector."); + if (requiresProjector()) { + ASTRA_CONFIG_CHECK(m_pSinogram->getGeometry()->isEqual(m_pProjector->getProjectionGeometry()), "Reconstruction2D", "Projection Data not compatible with the specified Projector."); + ASTRA_CONFIG_CHECK(m_pReconstruction->getGeometry()->isEqual(m_pProjector->getVolumeGeometry()), "Reconstruction2D", "Reconstruction Data not compatible with the specified Projector."); + } // success return true; -- cgit v1.2.3 From f7e01f5a3ca7780a29d1fbc3790e527c310cc7f8 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 9 Oct 2015 15:55:11 +0200 Subject: Fix loop bounds in (unused) Float32ProjectionData3D arithmetic functions --- src/Float32ProjectionData3D.cpp | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) (limited to 'src') diff --git a/src/Float32ProjectionData3D.cpp b/src/Float32ProjectionData3D.cpp index 2bd0447..680ad55 100644 --- a/src/Float32ProjectionData3D.cpp +++ b/src/Float32ProjectionData3D.cpp @@ -53,13 +53,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const CFloat32Pro CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); #ifdef _DEBUG CProjectionGeometry3D * pDataGeometry = _data.getGeometry(); - int iThisProjectionDetectorCount = pThisGeometry->getDetectorRowCount() * pThisGeometry->getDetectorColCount(); - int iDataProjectionDetectorCount = pDataGeometry->getDetectorRowCount() * pDataGeometry->getDetectorColCount(); + int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount(); ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount()); - ASTRA_ASSERT(iThisProjectionDetectorCount == iDataProjectionDetectorCount); + ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount); #endif for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) @@ -67,7 +67,7 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const CFloat32Pro CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex]; @@ -91,13 +91,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const CFloat32Pro CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); #ifdef _DEBUG CProjectionGeometry3D * pDataGeometry = _data.getGeometry(); - int iThisProjectionDetectorCount = pThisGeometry->getDetectorRowCount() * pThisGeometry->getDetectorColCount(); - int iDataProjectionDetectorCount = pDataGeometry->getDetectorRowCount() * pDataGeometry->getDetectorColCount(); + int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount(); ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount()); - ASTRA_ASSERT(iThisProjectionDetectorCount == iDataProjectionDetectorCount); + ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount); #endif for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) @@ -105,7 +105,7 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const CFloat32Pro CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex]; @@ -129,13 +129,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const CFloat32Pro CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); #ifdef _DEBUG CProjectionGeometry3D * pDataGeometry = _data.getGeometry(); - int iThisProjectionDetectorCount = pThisGeometry->getDetectorRowCount() * pThisGeometry->getDetectorColCount(); - int iDataProjectionDetectorCount = pDataGeometry->getDetectorRowCount() * pDataGeometry->getDetectorColCount(); + int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount(); ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount()); - ASTRA_ASSERT(iThisProjectionDetectorCount == iDataProjectionDetectorCount); + ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount); #endif for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) @@ -143,7 +143,7 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const CFloat32Pro CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex]; @@ -167,12 +167,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const float32& _f CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) { CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; @@ -194,12 +195,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator/=(const float32& _f CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) { CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; @@ -221,12 +223,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const float32& _f CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) { CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; @@ -248,12 +251,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const float32& _f CProjectionGeometry3D * pThisGeometry = getGeometry(); int iProjectionCount = pThisGeometry->getProjectionCount(); + int iDetectorCount = pThisGeometry->getDetectorTotCount(); for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++) { CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex); - for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++) + for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++) { float32 fThisValue = pThisProjection->getData()[iDetectorIndex]; -- cgit v1.2.3 From c7128284fdbbfa0d4a5cbc951b9cdeaf8f9b41e0 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 9 Oct 2015 16:10:45 +0200 Subject: Call check() function after initializing CUDA_FBP This would cause crashes when specifying invalid data. --- src/CudaFilteredBackProjectionAlgorithm.cpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp index aac96d6..6353c46 100644 --- a/src/CudaFilteredBackProjectionAlgorithm.cpp +++ b/src/CudaFilteredBackProjectionAlgorithm.cpp @@ -189,9 +189,7 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) m_pFBP = new AstraFBP; m_bAstraFBPInit = false; - // success - m_bIsInitialized = true; - return m_bIsInitialized; + return check(); } bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D * _pSinogram, CFloat32VolumeData2D * _pReconstruction, E_FBPFILTER _eFilter, const float * _pfFilter /* = NULL */, int _iFilterWidth /* = 0 */, int _iGPUIndex /* = 0 */, float _fFilterParameter /* = -1.0f */) @@ -241,7 +239,7 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D * m_fFilterParameter = _fFilterParameter; - return m_bIsInitialized; + return check(); } void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */) @@ -361,7 +359,7 @@ bool CCudaFilteredBackProjectionAlgorithm::check() ASTRA_CONFIG_CHECK(m_pReconstruction->isInitialized(), "FBP_CUDA", "Reconstruction Data Object Not Initialized."); // check gpu index - ASTRA_CONFIG_CHECK(m_iGPUIndex >= -1, "FBP_CUDA", "GPUIndex must be a non-negative integer."); + ASTRA_CONFIG_CHECK(m_iGPUIndex >= -1, "FBP_CUDA", "GPUIndex must be a non-negative integer or -1."); // check pixel supersampling ASTRA_CONFIG_CHECK(m_iPixelSuperSampling >= 0, "FBP_CUDA", "PixelSuperSampling must be a non-negative integer."); -- cgit v1.2.3 From fb44faa449990400861f1869b52f5afc8fefe01b Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 9 Oct 2015 16:19:54 +0200 Subject: Fix warning text --- src/CudaReconstructionAlgorithm2D.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp index bccdb43..2d023b7 100644 --- a/src/CudaReconstructionAlgorithm2D.cpp +++ b/src/CudaReconstructionAlgorithm2D.cpp @@ -163,9 +163,9 @@ bool CCudaReconstructionAlgorithm2D::_check() if (!CReconstructionAlgorithm2D::_check()) return false; - 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 >= -1, "SIRT_CUDA", "GPUIndex must be a non-negative integer."); + ASTRA_CONFIG_CHECK(m_iDetectorSuperSampling >= 1, "CudaReconstructionAlgorithm2D", "DetectorSuperSampling must be a positive integer."); + ASTRA_CONFIG_CHECK(m_iPixelSuperSampling >= 1, "CudaReconstructionAlgorithm2D", "PixelSuperSampling must be a positive integer."); + ASTRA_CONFIG_CHECK(m_iGPUIndex >= -1, "CudaReconstructionAlgorithm2D", "GPUIndex must be a non-negative integer or -1."); // check restrictions // TODO: check restrictions built into cuda code -- cgit v1.2.3 From 4298c2f212aac1e76f1f123ab199749a9a668415 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 9 Oct 2015 16:40:39 +0200 Subject: Give a warning on ignored Min/MaxContraint in some CUDA algorithms. Previously it would fail an assertion. --- src/CudaReconstructionAlgorithm2D.cpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp index 71b6637..18627fc 100644 --- a/src/CudaReconstructionAlgorithm2D.cpp +++ b/src/CudaReconstructionAlgorithm2D.cpp @@ -462,10 +462,18 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations) ASTRA_ASSERT(ok); - if (m_bUseMinConstraint) - ok &= m_pAlgo->setMinConstraint(m_fMinValue); - if (m_bUseMaxConstraint) - ok &= m_pAlgo->setMaxConstraint(m_fMaxValue); + if (m_bUseMinConstraint) { + bool ret = m_pAlgo->setMinConstraint(m_fMinValue); + if (!ret) { + ASTRA_WARN("This algorithm ignores MinConstraint"); + } + } + if (m_bUseMaxConstraint) { + bool ret= m_pAlgo->setMaxConstraint(m_fMaxValue); + if (!ret) { + ASTRA_WARN("This algorithm ignores MaxConstraint"); + } + } ok &= m_pAlgo->iterate(_iNrIterations); ASTRA_ASSERT(ok); -- cgit v1.2.3 From 21d08656ead6f974f83b0a02b03b105a7cd617a8 Mon Sep 17 00:00:00 2001 From: "Daniel M. Pelt" Date: Tue, 13 Oct 2015 17:02:09 +0200 Subject: Do not reuse va_list when logging both to screen and file --- src/Logging.cpp | 60 ++++++++++++++++++++++++++++++++++++++++++--------------- 1 file changed, 44 insertions(+), 16 deletions(-) (limited to 'src') diff --git a/src/Logging.cpp b/src/Logging.cpp index 8290ca0..cd7e3f0 100644 --- a/src/Logging.cpp +++ b/src/Logging.cpp @@ -70,37 +70,65 @@ void CLogger::disable() 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); + va_list ap, apf; + if(m_bEnabledScreen){ + va_start(ap, fmt); + clog_debug(sfile,sline,0,fmt,ap); + va_end(ap); + } + if(m_bEnabledFile && m_bFileProvided){ + va_start(apf, fmt); + clog_debug(sfile,sline,1,fmt,apf); + va_end(apf); + } } 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); + va_list ap, apf; + if(m_bEnabledScreen){ + va_start(ap, fmt); + clog_info(sfile,sline,0,fmt,ap); + va_end(ap); + } + if(m_bEnabledFile && m_bFileProvided){ + va_start(apf, fmt); + clog_info(sfile,sline,1,fmt,apf); + va_end(apf); + } } 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); + va_list ap, apf; + if(m_bEnabledScreen){ + va_start(ap, fmt); + clog_warn(sfile,sline,0,fmt,ap); + va_end(ap); + } + if(m_bEnabledFile && m_bFileProvided){ + va_start(apf, fmt); + clog_warn(sfile,sline,1,fmt,apf); + va_end(apf); + } } 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); + va_list ap, apf; + if(m_bEnabledScreen){ + va_start(ap, fmt); + clog_error(sfile,sline,0,fmt,ap); + va_end(ap); + } + if(m_bEnabledFile && m_bFileProvided){ + va_start(apf, fmt); + clog_error(sfile,sline,1,fmt,apf); + va_end(apf); + } } void CLogger::_setLevel(int id, log_level m_eLevel) -- cgit v1.2.3 From 07c31b932078544205d61551edd4a66f69be30ae Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 2 Dec 2015 11:25:59 +0100 Subject: Avoid unnecessary include in header --- src/PluginAlgorithm.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index e79c77b..8f7dfc5 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -37,9 +37,13 @@ $Id$ #include #include +#include +#include "bytesobject.h" + namespace astra { + void logPythonError(){ if(PyErr_Occurred()){ PyObject *ptype, *pvalue, *ptraceback; @@ -394,4 +398,4 @@ PyObject* XMLNode2dict(XMLNode node){ } } -#endif \ No newline at end of file +#endif -- cgit v1.2.3 From 3ea35516aceec4f5817871a00008b109777ebb13 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 30 Nov 2015 16:07:52 +0100 Subject: Disable error-prone checks --- src/VolumeGeometry3D.cpp | 2 ++ 1 file changed, 2 insertions(+) (limited to 'src') diff --git a/src/VolumeGeometry3D.cpp b/src/VolumeGeometry3D.cpp index a1cf424..3de146f 100644 --- a/src/VolumeGeometry3D.cpp +++ b/src/VolumeGeometry3D.cpp @@ -45,6 +45,7 @@ bool CVolumeGeometry3D::_check() ASTRA_CONFIG_CHECK(m_fWindowMinZ < m_fWindowMaxZ, "VolumeGeometry3D", "WindowMinZ should be lower than WindowMaxZ."); ASTRA_CONFIG_CHECK(m_iGridTotCount == (m_iGridColCount * m_iGridRowCount * m_iGridSliceCount), "VolumeGeometry3D", "Internal configuration error."); +#if 0 ASTRA_CONFIG_CHECK(m_fWindowLengthX == (m_fWindowMaxX - m_fWindowMinX), "VolumeGeometry3D", "Internal configuration error."); ASTRA_CONFIG_CHECK(m_fWindowLengthY == (m_fWindowMaxY - m_fWindowMinY), "VolumeGeometry3D", "Internal configuration error."); ASTRA_CONFIG_CHECK(m_fWindowLengthZ == (m_fWindowMaxZ - m_fWindowMinZ), "VolumeGeometry3D", "Internal configuration error."); @@ -57,6 +58,7 @@ bool CVolumeGeometry3D::_check() ASTRA_CONFIG_CHECK(m_fDivPixelLengthX == (1.0f / m_fPixelLengthX), "VolumeGeometry3D", "Internal configuration error."); ASTRA_CONFIG_CHECK(m_fDivPixelLengthY == (1.0f / m_fPixelLengthY), "VolumeGeometry3D", "Internal configuration error."); ASTRA_CONFIG_CHECK(m_fDivPixelLengthZ == (1.0f / m_fPixelLengthZ), "VolumeGeometry3D", "Internal configuration error."); +#endif return true; } -- cgit v1.2.3 From b14fb531ad9ae3d565f2cf28f5506408ab10dbed Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 18 Nov 2015 11:26:15 +0100 Subject: Add CompositeGeometryManager This handles FP and BP operations on multiple data objects at once, splitting them to fit in GPU memory where necessary. --- src/CompositeGeometryManager.cpp | 884 +++++++++++++++++++++++++++++++ src/ConeProjectionGeometry3D.cpp | 92 +++- src/ConeVecProjectionGeometry3D.cpp | 58 +- src/CudaBackProjectionAlgorithm3D.cpp | 8 + src/CudaForwardProjectionAlgorithm3D.cpp | 9 + src/GeometryUtil3D.cpp | 172 ++++++ src/ParallelProjectionGeometry3D.cpp | 81 ++- src/ParallelVecProjectionGeometry3D.cpp | 61 ++- 8 files changed, 1351 insertions(+), 14 deletions(-) create mode 100644 src/CompositeGeometryManager.cpp (limited to 'src') diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp new file mode 100644 index 0000000..fc8bc2e --- /dev/null +++ b/src/CompositeGeometryManager.cpp @@ -0,0 +1,884 @@ +/* +----------------------------------------------------------------------- +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 . + +----------------------------------------------------------------------- +*/ + +#include "astra/CompositeGeometryManager.h" + +#ifdef ASTRA_CUDA + +#include "astra/GeometryUtil3D.h" +#include "astra/VolumeGeometry3D.h" +#include "astra/ConeProjectionGeometry3D.h" +#include "astra/ConeVecProjectionGeometry3D.h" +#include "astra/ParallelProjectionGeometry3D.h" +#include "astra/ParallelVecProjectionGeometry3D.h" +#include "astra/Projector3D.h" +#include "astra/CudaProjector3D.h" +#include "astra/Float32ProjectionData3DMemory.h" +#include "astra/Float32VolumeData3DMemory.h" +#include "astra/Logging.h" + +#include "../cuda/3d/mem3d.h" + +#include + +namespace astra { + +// JOB: +// +// VolumePart +// ProjectionPart +// FP-or-BP +// SET-or-ADD + + +// Running a set of jobs: +// +// [ Assume OUTPUT Parts in a single JobSet don't alias?? ] +// Group jobs by output Part +// One thread per group? + +// Automatically split parts if too large +// Performance model for odd-sized tasks? +// Automatically split parts if not enough tasks to fill available GPUs + + +// Splitting: +// Constraints: +// number of sub-parts divisible by N +// max size of sub-parts + +// For splitting on both input and output side: +// How to divide up memory? (Optimization problem; compute/benchmark) +// (First approach: 0.5/0.5) + + + +bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div, TJobSet &split) +{ + split.clear(); + + for (TJobSet::const_iterator i = jobs.begin(); i != jobs.end(); ++i) + { + CPart* pOutput = i->first; + const TJobList &L = i->second; + + // 1. Split output part + // 2. Per sub-part: + // a. reduce input part + // b. split input part + // c. create jobs for new (input,output) subparts + + TPartList splitOutput = pOutput->split(maxSize/3, div); + + for (TJobList::const_iterator j = L.begin(); j != L.end(); ++j) + { + const SJob &job = *j; + + for (TPartList::iterator i_out = splitOutput.begin(); + i_out != splitOutput.end(); ++i_out) + { + boost::shared_ptr outputPart = *i_out; + split[outputPart.get()] = TJobList(); + + SJob newjob; + newjob.pOutput = outputPart; + newjob.eType = j->eType; + newjob.eMode = j->eMode; + newjob.pProjector = j->pProjector; + + CPart* input = job.pInput->reduce(outputPart.get()); + + if (input->getSize() == 0) { + ASTRA_DEBUG("Empty input"); + newjob.eType = SJob::JOB_NOP; + split[outputPart.get()].push_back(newjob); + continue; + } + + size_t remainingSize = ( maxSize - outputPart->getSize() ) / 2; + + TPartList splitInput = input->split(remainingSize, 1); + delete input; + ASTRA_DEBUG("Input split into %d parts", splitInput.size()); + + for (TPartList::iterator i_in = splitInput.begin(); + i_in != splitInput.end(); ++i_in) + { + newjob.pInput = *i_in; + + split[outputPart.get()].push_back(newjob); + + // Second and later (input) parts should always be added to + // output of first (input) part. + newjob.eMode = SJob::MODE_ADD; + } + + + } + + } + } + + return true; +} + +CCompositeGeometryManager::CPart::CPart(const CPart& other) +{ + eType = other.eType; + pData = other.pData; + subX = other.subX; + subY = other.subY; + subZ = other.subZ; +} + +CCompositeGeometryManager::CVolumePart::CVolumePart(const CVolumePart& other) + : CPart(other) +{ + pGeom = other.pGeom->clone(); +} + +CCompositeGeometryManager::CVolumePart::~CVolumePart() +{ + delete pGeom; +} + +void CCompositeGeometryManager::CVolumePart::getDims(size_t &x, size_t &y, size_t &z) +{ + if (!pGeom) { + x = y = z = 0; + return; + } + + x = pGeom->getGridColCount(); + y = pGeom->getGridRowCount(); + z = pGeom->getGridSliceCount(); +} + +size_t CCompositeGeometryManager::CPart::getSize() +{ + size_t x, y, z; + getDims(x, y, z); + return x * y * z; +} + + + +CCompositeGeometryManager::CPart* CCompositeGeometryManager::CVolumePart::reduce(const CPart *_other) +{ + const CProjectionPart *other = dynamic_cast(_other); + assert(other); + + // TODO: Is 0.5 sufficient? + double umin = -0.5; + double umax = other->pGeom->getDetectorColCount() + 0.5; + double vmin = -0.5; + double vmax = other->pGeom->getDetectorRowCount() + 0.5; + + double uu[4]; + double vv[4]; + uu[0] = umin; vv[0] = vmin; + uu[1] = umin; vv[1] = vmax; + uu[2] = umax; vv[2] = vmin; + uu[3] = umax; vv[3] = vmax; + + double pixx = pGeom->getPixelLengthX(); + double pixy = pGeom->getPixelLengthY(); + double pixz = pGeom->getPixelLengthZ(); + + double xmin = pGeom->getWindowMinX() - 0.5 * pixx; + double xmax = pGeom->getWindowMaxX() + 0.5 * pixx; + double ymin = pGeom->getWindowMinY() - 0.5 * pixy; + double ymax = pGeom->getWindowMaxY() + 0.5 * pixy; + + // NB: Flipped + double zmax = pGeom->getWindowMinZ() - 2.5 * pixz; + double zmin = pGeom->getWindowMaxZ() + 2.5 * pixz; + + // TODO: This isn't as tight as it could be. + // In particular it won't detect the detector being + // missed entirely on the u side. + + for (int i = 0; i < other->pGeom->getProjectionCount(); ++i) { + for (int j = 0; j < 4; ++j) { + double px, py, pz; + + other->pGeom->backprojectPointX(i, uu[j], vv[j], xmin, py, pz); + //ASTRA_DEBUG("%f %f (%f - %f)", py, pz, ymin, ymax); + if (pz < zmin) zmin = pz; + if (pz > zmax) zmax = pz; + other->pGeom->backprojectPointX(i, uu[j], vv[j], xmax, py, pz); + //ASTRA_DEBUG("%f %f (%f - %f)", py, pz, ymin, ymax); + if (pz < zmin) zmin = pz; + if (pz > zmax) zmax = pz; + + other->pGeom->backprojectPointY(i, uu[j], vv[j], ymin, px, pz); + //ASTRA_DEBUG("%f %f (%f - %f)", px, pz, xmin, xmax); + if (pz < zmin) zmin = pz; + if (pz > zmax) zmax = pz; + other->pGeom->backprojectPointY(i, uu[j], vv[j], ymax, px, pz); + //ASTRA_DEBUG("%f %f (%f - %f)", px, pz, xmin, xmax); + if (pz < zmin) zmin = pz; + if (pz > zmax) zmax = pz; + } + } + + //ASTRA_DEBUG("coord extent: %f - %f", zmin, zmax); + + zmin = (zmin - pixz - pGeom->getWindowMinZ()) / pixz; + zmax = (zmax + pixz - pGeom->getWindowMinZ()) / pixz; + + int _zmin = (int)floor(zmin); + int _zmax = (int)ceil(zmax); + + //ASTRA_DEBUG("index extent: %d - %d", _zmin, _zmax); + + if (_zmin < 0) + _zmin = 0; + if (_zmax > pGeom->getGridSliceCount()) + _zmax = pGeom->getGridSliceCount(); + + if (_zmax <= _zmin) { + _zmin = _zmax = 0; + } + //ASTRA_DEBUG("adjusted extent: %d - %d", _zmin, _zmax); + + CVolumePart *sub = new CVolumePart(); + sub->subX = this->subX; + sub->subY = this->subY; + sub->subZ = this->subZ + _zmin; + sub->pData = pData; + + if (_zmin == _zmax) { + sub->pGeom = 0; + } else { + sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(), + pGeom->getGridRowCount(), + _zmax - _zmin, + pGeom->getWindowMinX(), + pGeom->getWindowMinY(), + pGeom->getWindowMinZ() + _zmin * pixz, + pGeom->getWindowMaxX(), + pGeom->getWindowMaxY(), + pGeom->getWindowMinZ() + _zmax * pixz); + } + + ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d", this->subZ, this->subZ + pGeom->getGridSliceCount(), this->subZ + _zmin, this->subZ + _zmax); + + return sub; +} + + + +static size_t ceildiv(size_t a, size_t b) { + return (a + b - 1) / b; +} + +static size_t computeVerticalSplit(size_t maxBlock, int div, size_t sliceCount) +{ + size_t blockSize = maxBlock; + size_t blockCount = ceildiv(sliceCount, blockSize); + + // Increase number of blocks to be divisible by div + size_t divCount = div * ceildiv(blockCount, div); + + // If divCount is above sqrt(number of slices), then + // we can't guarantee divisibility by div, but let's try anyway + if (ceildiv(sliceCount, ceildiv(sliceCount, divCount)) % div == 0) { + blockCount = divCount; + } else { + // If divisibility isn't achievable, we may want to optimize + // differently. + // TODO: Figure out how to model and optimize this. + } + + // Final adjustment to make blocks more evenly sized + // (This can't make the blocks larger) + blockSize = ceildiv(sliceCount, blockCount); + + ASTRA_DEBUG("%ld %ld -> %ld * %ld\n", sliceCount, maxBlock, blockCount, blockSize); + + assert(blockSize <= maxBlock); + assert((divCount * divCount > sliceCount) || (blockCount % div) == 0); + + return blockSize; +} + +template +static V* getProjectionVectors(const P* geom); + +template<> +SConeProjection* getProjectionVectors(const CConeProjectionGeometry3D* pProjGeom) +{ + return genConeProjections(pProjGeom->getProjectionCount(), + pProjGeom->getDetectorColCount(), + pProjGeom->getDetectorRowCount(), + pProjGeom->getOriginSourceDistance(), + pProjGeom->getOriginDetectorDistance(), + pProjGeom->getDetectorSpacingX(), + pProjGeom->getDetectorSpacingY(), + pProjGeom->getProjectionAngles()); +} + +template<> +SConeProjection* getProjectionVectors(const CConeVecProjectionGeometry3D* pProjGeom) +{ + int nth = pProjGeom->getProjectionCount(); + + SConeProjection* pProjs = new SConeProjection[nth]; + for (int i = 0; i < nth; ++i) + pProjs[i] = pProjGeom->getProjectionVectors()[i]; + + return pProjs; +} + +template<> +SPar3DProjection* getProjectionVectors(const CParallelProjectionGeometry3D* pProjGeom) +{ + return genPar3DProjections(pProjGeom->getProjectionCount(), + pProjGeom->getDetectorColCount(), + pProjGeom->getDetectorRowCount(), + pProjGeom->getDetectorSpacingX(), + pProjGeom->getDetectorSpacingY(), + pProjGeom->getProjectionAngles()); +} + +template<> +SPar3DProjection* getProjectionVectors(const CParallelVecProjectionGeometry3D* pProjGeom) +{ + int nth = pProjGeom->getProjectionCount(); + + SPar3DProjection* pProjs = new SPar3DProjection[nth]; + for (int i = 0; i < nth; ++i) + pProjs[i] = pProjGeom->getProjectionVectors()[i]; + + return pProjs; +} + + +template +static void translateProjectionVectors(V* pProjs, int count, double dv) +{ + for (int i = 0; i < count; ++i) { + pProjs[i].fDetSX += dv * pProjs[i].fDetVX; + pProjs[i].fDetSY += dv * pProjs[i].fDetVY; + pProjs[i].fDetSZ += dv * pProjs[i].fDetVZ; + } +} + + + +static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry3D* pProjGeom, int v, int size) +{ + // First convert to vectors, then translate, then convert into new object + + const CConeProjectionGeometry3D* conegeom = dynamic_cast(pProjGeom); + const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast(pProjGeom); + const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast(pProjGeom); + const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast(pProjGeom); + + if (conegeom || conevec3dgeom) { + SConeProjection* pConeProjs; + if (conegeom) { + pConeProjs = getProjectionVectors(conegeom); + } else { + pConeProjs = getProjectionVectors(conevec3dgeom); + } + + translateProjectionVectors(pConeProjs, pProjGeom->getProjectionCount(), v); + + CProjectionGeometry3D* ret = new CConeVecProjectionGeometry3D(pProjGeom->getProjectionCount(), + size, + pProjGeom->getDetectorColCount(), + pConeProjs); + + + delete[] pConeProjs; + return ret; + } else { + assert(par3dgeom || parvec3dgeom); + SPar3DProjection* pParProjs; + if (par3dgeom) { + pParProjs = getProjectionVectors(par3dgeom); + } else { + pParProjs = getProjectionVectors(parvec3dgeom); + } + + translateProjectionVectors(pParProjs, pProjGeom->getProjectionCount(), v); + + CProjectionGeometry3D* ret = new CParallelVecProjectionGeometry3D(pProjGeom->getProjectionCount(), + size, + pProjGeom->getDetectorColCount(), + pParProjs); + + delete[] pParProjs; + return ret; + } + +} + + + +// split self into sub-parts: +// - each no bigger than maxSize +// - number of sub-parts is divisible by div +// - maybe all approximately the same size? +CCompositeGeometryManager::TPartList CCompositeGeometryManager::CVolumePart::split(size_t maxSize, int div) +{ + TPartList ret; + + if (true) { + // Split in vertical direction only at first, until we figure out + // a model for splitting in other directions + + size_t sliceSize = ((size_t) pGeom->getGridColCount()) * pGeom->getGridRowCount(); + int sliceCount = pGeom->getGridSliceCount(); + size_t blockSize = computeVerticalSplit(maxSize / sliceSize, div, sliceCount); + + int rem = sliceCount % blockSize; + + ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); + + for (int z = -(rem / 2); z < sliceCount; z += blockSize) { + int newsubZ = z; + if (newsubZ < 0) newsubZ = 0; + int endZ = z + blockSize; + if (endZ > sliceCount) endZ = sliceCount; + int size = endZ - newsubZ; + + CVolumePart *sub = new CVolumePart(); + sub->subX = this->subX; + sub->subY = this->subY; + sub->subZ = this->subZ + newsubZ; + + ASTRA_DEBUG("VolumePart split %d %d %d -> %p", sub->subX, sub->subY, sub->subZ, (void*)sub); + + double shift = pGeom->getPixelLengthZ() * newsubZ; + + sub->pData = pData; + sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(), + pGeom->getGridRowCount(), + size, + pGeom->getWindowMinX(), + pGeom->getWindowMinY(), + pGeom->getWindowMinZ() + shift, + pGeom->getWindowMaxX(), + pGeom->getWindowMaxY(), + pGeom->getWindowMinZ() + shift + size * pGeom->getPixelLengthZ()); + + ret.push_back(boost::shared_ptr(sub)); + } + } + + return ret; +} + +CCompositeGeometryManager::CVolumePart* CCompositeGeometryManager::CVolumePart::clone() const +{ + return new CVolumePart(*this); +} + +CCompositeGeometryManager::CProjectionPart::CProjectionPart(const CProjectionPart& other) + : CPart(other) +{ + pGeom = other.pGeom->clone(); +} + +CCompositeGeometryManager::CProjectionPart::~CProjectionPart() +{ + delete pGeom; +} + +void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, size_t &z) +{ + if (!pGeom) { + x = y = z = 0; + return; + } + + x = pGeom->getDetectorColCount(); + y = pGeom->getProjectionCount(); + z = pGeom->getDetectorRowCount(); +} + + +CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::reduce(const CPart *_other) +{ + const CVolumePart *other = dynamic_cast(_other); + assert(other); + + double vmin_g, vmax_g; + + // reduce self to only cover intersection with projection of VolumePart + // (Project corners of volume, take bounding box) + + for (int i = 0; i < pGeom->getProjectionCount(); ++i) { + + double vol_u[8]; + double vol_v[8]; + + double pixx = other->pGeom->getPixelLengthX(); + double pixy = other->pGeom->getPixelLengthY(); + double pixz = other->pGeom->getPixelLengthZ(); + + // TODO: Is 0.5 sufficient? + double xmin = other->pGeom->getWindowMinX() - 0.5 * pixx; + double xmax = other->pGeom->getWindowMaxX() + 0.5 * pixx; + double ymin = other->pGeom->getWindowMinY() - 0.5 * pixy; + double ymax = other->pGeom->getWindowMaxY() + 0.5 * pixy; + double zmin = other->pGeom->getWindowMinZ() - 0.5 * pixz; + double zmax = other->pGeom->getWindowMaxZ() + 0.5 * pixz; + + pGeom->projectPoint(xmin, ymin, zmin, i, vol_u[0], vol_v[0]); + pGeom->projectPoint(xmin, ymin, zmax, i, vol_u[1], vol_v[1]); + pGeom->projectPoint(xmin, ymax, zmin, i, vol_u[2], vol_v[2]); + pGeom->projectPoint(xmin, ymax, zmax, i, vol_u[3], vol_v[3]); + pGeom->projectPoint(xmax, ymin, zmin, i, vol_u[4], vol_v[4]); + pGeom->projectPoint(xmax, ymin, zmax, i, vol_u[5], vol_v[5]); + pGeom->projectPoint(xmax, ymax, zmin, i, vol_u[6], vol_v[6]); + pGeom->projectPoint(xmax, ymax, zmax, i, vol_u[7], vol_v[7]); + + double vmin = vol_v[0]; + double vmax = vol_v[0]; + + for (int j = 1; j < 8; ++j) { + if (vol_v[j] < vmin) + vmin = vol_v[j]; + if (vol_v[j] > vmax) + vmax = vol_v[j]; + } + + if (i == 0 || vmin < vmin_g) + vmin_g = vmin; + if (i == 0 || vmax > vmax_g) + vmax_g = vmax; + } + + // fprintf(stderr, "v extent: %f %f\n", vmin_g, vmax_g); + + int _vmin = (int)floor(vmin_g - 1.0f); + int _vmax = (int)ceil(vmax_g + 1.0f); + if (_vmin < 0) + _vmin = 0; + if (_vmax > pGeom->getDetectorRowCount()) + _vmax = pGeom->getDetectorRowCount(); + + if (_vmin >= _vmax) { + _vmin = _vmax = 0; + } + + CProjectionPart *sub = new CProjectionPart(); + sub->subX = this->subX; + sub->subY = this->subY; + sub->subZ = this->subZ + _vmin; + + sub->pData = pData; + + if (_vmin == _vmax) { + sub->pGeom = 0; + } else { + sub->pGeom = getSubProjectionGeometry(pGeom, _vmin, _vmax - _vmin); + } + + ASTRA_DEBUG("Reduce projection from %d - %d to %d - %d", this->subZ, this->subZ + pGeom->getDetectorRowCount(), this->subZ + _vmin, this->subZ + _vmax); + + return sub; +} + + +CCompositeGeometryManager::TPartList CCompositeGeometryManager::CProjectionPart::split(size_t maxSize, int div) +{ + TPartList ret; + + if (true) { + // Split in vertical direction only at first, until we figure out + // a model for splitting in other directions + + size_t sliceSize = ((size_t) pGeom->getDetectorColCount()) * pGeom->getProjectionCount(); + int sliceCount = pGeom->getDetectorRowCount(); + size_t blockSize = computeVerticalSplit(maxSize / sliceSize, div, sliceCount); + + int rem = sliceCount % blockSize; + + for (int z = -(rem / 2); z < sliceCount; z += blockSize) { + int newsubZ = z; + if (newsubZ < 0) newsubZ = 0; + int endZ = z + blockSize; + if (endZ > sliceCount) endZ = sliceCount; + int size = endZ - newsubZ; + + CProjectionPart *sub = new CProjectionPart(); + sub->subX = this->subX; + sub->subY = this->subY; + sub->subZ = this->subZ + newsubZ; + + ASTRA_DEBUG("ProjectionPart split %d %d %d -> %p", sub->subX, sub->subY, sub->subZ, (void*)sub); + + sub->pData = pData; + + sub->pGeom = getSubProjectionGeometry(pGeom, newsubZ, size); + + ret.push_back(boost::shared_ptr(sub)); + } + } + + return ret; + +} + +CCompositeGeometryManager::CProjectionPart* CCompositeGeometryManager::CProjectionPart::clone() const +{ + return new CProjectionPart(*this); +} + + +bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, + CFloat32ProjectionData3DMemory *pProjData) +{ + ASTRA_DEBUG("CCompositeGeometryManager::doFP"); + // Create single job for FP + // Run result + + CVolumePart *input = new CVolumePart(); + input->pData = pVolData; + input->subX = 0; + input->subY = 0; + input->subZ = 0; + input->pGeom = pVolData->getGeometry()->clone(); + ASTRA_DEBUG("Main FP VolumePart -> %p", (void*)input); + + CProjectionPart *output = new CProjectionPart(); + output->pData = pProjData; + output->subX = 0; + output->subY = 0; + output->subZ = 0; + output->pGeom = pProjData->getGeometry()->clone(); + ASTRA_DEBUG("Main FP ProjectionPart -> %p", (void*)output); + + SJob FP; + FP.pInput = boost::shared_ptr(input); + FP.pOutput = boost::shared_ptr(output); + FP.pProjector = pProjector; + FP.eType = SJob::JOB_FP; + FP.eMode = SJob::MODE_SET; + + TJobList L; + L.push_back(FP); + + return doJobs(L); +} + +bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData, + CFloat32ProjectionData3DMemory *pProjData) +{ + ASTRA_DEBUG("CCompositeGeometryManager::doBP"); + // Create single job for BP + // Run result + + CProjectionPart *input = new CProjectionPart(); + input->pData = pProjData; + input->subX = 0; + input->subY = 0; + input->subZ = 0; + input->pGeom = pProjData->getGeometry()->clone(); + + CVolumePart *output = new CVolumePart(); + output->pData = pVolData; + output->subX = 0; + output->subY = 0; + output->subZ = 0; + output->pGeom = pVolData->getGeometry()->clone(); + + SJob BP; + BP.pInput = boost::shared_ptr(input); + BP.pOutput = boost::shared_ptr(output); + BP.pProjector = pProjector; + BP.eType = SJob::JOB_BP; + BP.eMode = SJob::MODE_SET; + + TJobList L; + L.push_back(BP); + + return doJobs(L); +} + + + +bool CCompositeGeometryManager::doJobs(TJobList &jobs) +{ + ASTRA_DEBUG("CCompositeGeometryManager::doJobs"); + + // Sort job list into job set by output part + TJobSet jobset; + + for (TJobList::iterator i = jobs.begin(); i != jobs.end(); ++i) { + jobset[i->pOutput.get()].push_back(*i); + } + + size_t maxSize = astraCUDA3d::availableGPUMemory(); + if (maxSize == 0) { + ASTRA_WARN("Unable to get available GPU memory. Defaulting to 1GB."); + maxSize = 1024 * 1024 * 1024; + } else { + ASTRA_DEBUG("Detected %lu bytes of GPU memory", maxSize); + } + maxSize = (maxSize * 9) / 10; + + maxSize /= sizeof(float); + int div = 1; + + // TODO: Multi-GPU support + + // Split jobs to fit + TJobSet split; + splitJobs(jobset, maxSize, div, split); + jobset.clear(); + + // Run jobs + + for (TJobSet::iterator iter = split.begin(); iter != split.end(); ++iter) { + + CPart* output = iter->first; + TJobList& L = iter->second; + + assert(!L.empty()); + + bool zero = L.begin()->eMode == SJob::MODE_SET; + + size_t outx, outy, outz; + output->getDims(outx, outy, outz); + + if (L.begin()->eType == SJob::JOB_NOP) { + // just zero output? + if (zero) { + for (size_t z = 0; z < outz; ++z) { + for (size_t y = 0; y < outy; ++y) { + float* ptr = output->pData->getData(); + ptr += (z + output->subX) * (size_t)output->pData->getHeight() * (size_t)output->pData->getWidth(); + ptr += (y + output->subY) * (size_t)output->pData->getWidth(); + ptr += output->subX; + memset(ptr, 0, sizeof(float) * outx); + } + } + } + continue; + } + + + astraCUDA3d::SSubDimensions3D dstdims; + dstdims.nx = output->pData->getWidth(); + dstdims.pitch = dstdims.nx; + dstdims.ny = output->pData->getHeight(); + dstdims.nz = output->pData->getDepth(); + dstdims.subnx = outx; + dstdims.subny = outy; + dstdims.subnz = outz; + ASTRA_DEBUG("dstdims: %d,%d,%d in %d,%d,%d", dstdims.subnx, dstdims.subny, dstdims.subnz, dstdims.nx, dstdims.ny, dstdims.nz); + dstdims.subx = output->subX; + dstdims.suby = output->subY; + dstdims.subz = output->subZ; + float *dst = output->pData->getData(); + + astraCUDA3d::MemHandle3D outputMem = astraCUDA3d::allocateGPUMemory(outx, outy, outz, zero ? astraCUDA3d::INIT_ZERO : astraCUDA3d::INIT_NO); + bool ok = outputMem; + + for (TJobList::iterator i = L.begin(); i != L.end(); ++i) { + SJob &j = *i; + + assert(j.pInput); + + CCudaProjector3D *projector = dynamic_cast(j.pProjector); + Cuda3DProjectionKernel projKernel = ker3d_default; + int detectorSuperSampling = 1; + int voxelSuperSampling = 1; + if (projector) { + projKernel = projector->getProjectionKernel(); + detectorSuperSampling = projector->getDetectorSuperSampling(); + voxelSuperSampling = projector->getVoxelSuperSampling(); + } + + size_t inx, iny, inz; + j.pInput->getDims(inx, iny, inz); + astraCUDA3d::MemHandle3D inputMem = astraCUDA3d::allocateGPUMemory(inx, iny, inz, astraCUDA3d::INIT_NO); + + astraCUDA3d::SSubDimensions3D srcdims; + srcdims.nx = j.pInput->pData->getWidth(); + srcdims.pitch = srcdims.nx; + srcdims.ny = j.pInput->pData->getHeight(); + srcdims.nz = j.pInput->pData->getDepth(); + srcdims.subnx = inx; + srcdims.subny = iny; + srcdims.subnz = inz; + srcdims.subx = j.pInput->subX; + srcdims.suby = j.pInput->subY; + srcdims.subz = j.pInput->subZ; + const float *src = j.pInput->pData->getDataConst(); + + ok = astraCUDA3d::copyToGPUMemory(src, inputMem, srcdims); + if (!ok) ASTRA_ERROR("Error copying input data to GPU"); + + if (j.eType == SJob::JOB_FP) { + assert(dynamic_cast(j.pInput.get())); + assert(dynamic_cast(j.pOutput.get())); + + ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing FP"); + + ok = astraCUDA3d::FP(((CProjectionPart*)j.pOutput.get())->pGeom, outputMem, ((CVolumePart*)j.pInput.get())->pGeom, inputMem, detectorSuperSampling, projKernel); + if (!ok) ASTRA_ERROR("Error performing sub-FP"); + ASTRA_DEBUG("CCompositeGeometryManager::doJobs: FP done"); + } else if (j.eType == SJob::JOB_BP) { + assert(dynamic_cast(j.pOutput.get())); + assert(dynamic_cast(j.pInput.get())); + + ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing BP"); + + ok = astraCUDA3d::BP(((CProjectionPart*)j.pInput.get())->pGeom, inputMem, ((CVolumePart*)j.pOutput.get())->pGeom, outputMem, voxelSuperSampling); + if (!ok) ASTRA_ERROR("Error performing sub-BP"); + ASTRA_DEBUG("CCompositeGeometryManager::doJobs: BP done"); + } else { + assert(false); + } + + ok = astraCUDA3d::freeGPUMemory(inputMem); + if (!ok) ASTRA_ERROR("Error freeing GPU memory"); + + } + + ok = astraCUDA3d::copyFromGPUMemory(dst, outputMem, dstdims); + if (!ok) ASTRA_ERROR("Error copying output data from GPU"); + + ok = astraCUDA3d::freeGPUMemory(outputMem); + if (!ok) ASTRA_ERROR("Error freeing GPU memory"); + } + + return true; +} + + + +} + +#endif diff --git a/src/ConeProjectionGeometry3D.cpp b/src/ConeProjectionGeometry3D.cpp index dd22eba..18f0f8a 100644 --- a/src/ConeProjectionGeometry3D.cpp +++ b/src/ConeProjectionGeometry3D.cpp @@ -29,6 +29,7 @@ $Id$ #include "astra/ConeProjectionGeometry3D.h" #include "astra/Logging.h" +#include "astra/GeometryUtil3D.h" #include #include @@ -230,14 +231,14 @@ CVector3D CConeProjectionGeometry3D::getProjectionDirection(int _iProjectionInde return ret; } -void CConeProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, - int iAngleIndex, - float32 &fU, float32 &fV) const +void CConeProjectionGeometry3D::projectPoint(double fX, double fY, double fZ, + int iAngleIndex, + double &fU, double &fV) const { ASTRA_ASSERT(iAngleIndex >= 0); ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); - float alpha = m_pfProjectionAngles[iAngleIndex]; + double alpha = m_pfProjectionAngles[iAngleIndex]; // Project point onto optical axis @@ -245,14 +246,14 @@ void CConeProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, // 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; + double 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; + double fS = cos(alpha) * fX + sin(alpha) * fY; // Scale fS to detector plane fU = detectorOffsetXToColIndexFloat( (fS * (m_fOriginSourceDistance + m_fOriginDetectorDistance)) / fD ); @@ -261,5 +262,84 @@ void CConeProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, } +void CConeProjectionGeometry3D::backprojectPointX(int iAngleIndex, double fU, double fV, + double fX, double &fY, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SConeProjection *projs = genConeProjections(1, m_iDetectorColCount, m_iDetectorRowCount, + m_fOriginSourceDistance, + m_fOriginDetectorDistance, + m_fDetectorSpacingX, m_fDetectorSpacingY, + &m_pfProjectionAngles[iAngleIndex]); + + SConeProjection &proj = projs[0]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - proj.fSrcX) / (px - proj.fSrcX); + + fY = proj.fSrcY + a * (py - proj.fSrcY); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); + + delete[] projs; +} + +void CConeProjectionGeometry3D::backprojectPointY(int iAngleIndex, double fU, double fV, + double fY, double &fX, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SConeProjection *projs = genConeProjections(1, m_iDetectorColCount, m_iDetectorRowCount, + m_fOriginSourceDistance, + m_fOriginDetectorDistance, + m_fDetectorSpacingX, m_fDetectorSpacingY, + &m_pfProjectionAngles[iAngleIndex]); + + SConeProjection &proj = projs[0]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - proj.fSrcY) / (py - proj.fSrcY); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); + + delete[] projs; +} + +void CConeProjectionGeometry3D::backprojectPointZ(int iAngleIndex, double fU, double fV, + double fZ, double &fX, double &fY) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SConeProjection *projs = genConeProjections(1, m_iDetectorColCount, m_iDetectorRowCount, + m_fOriginSourceDistance, + m_fOriginDetectorDistance, + m_fDetectorSpacingX, m_fDetectorSpacingY, + &m_pfProjectionAngles[iAngleIndex]); + + SConeProjection &proj = projs[0]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - proj.fSrcZ) / (pz - proj.fSrcZ); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fY = proj.fSrcY + a * (py - proj.fSrcY); + + delete[] projs; +} + + } // end namespace astra diff --git a/src/ConeVecProjectionGeometry3D.cpp b/src/ConeVecProjectionGeometry3D.cpp index 47ed630..86e3bd6 100644 --- a/src/ConeVecProjectionGeometry3D.cpp +++ b/src/ConeVecProjectionGeometry3D.cpp @@ -241,9 +241,9 @@ 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, +void CConeVecProjectionGeometry3D::projectPoint(double fX, double fY, double fZ, int iAngleIndex, - float32 &fU, float32 &fV) const + double &fU, double &fV) const { ASTRA_ASSERT(iAngleIndex >= 0); ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); @@ -262,6 +262,60 @@ void CConeVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 } +void CConeVecProjectionGeometry3D::backprojectPointX(int iAngleIndex, double fU, double fV, + double fX, double &fY, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SConeProjection &proj = m_pProjectionAngles[iAngleIndex]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - proj.fSrcX) / (px - proj.fSrcX); + + fY = proj.fSrcY + a * (py - proj.fSrcY); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); +} + +void CConeVecProjectionGeometry3D::backprojectPointY(int iAngleIndex, double fU, double fV, + double fY, double &fX, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SConeProjection &proj = m_pProjectionAngles[iAngleIndex]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - proj.fSrcY) / (py - proj.fSrcY); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); +} + +void CConeVecProjectionGeometry3D::backprojectPointZ(int iAngleIndex, double fU, double fV, + double fZ, double &fX, double &fY) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SConeProjection &proj = m_pProjectionAngles[iAngleIndex]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - proj.fSrcZ) / (pz - proj.fSrcZ); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fY = proj.fSrcY + a * (py - proj.fSrcY); +} + //---------------------------------------------------------------------------------------- bool CConeVecProjectionGeometry3D::_check() diff --git a/src/CudaBackProjectionAlgorithm3D.cpp b/src/CudaBackProjectionAlgorithm3D.cpp index 8cf4c3b..ce8e111 100644 --- a/src/CudaBackProjectionAlgorithm3D.cpp +++ b/src/CudaBackProjectionAlgorithm3D.cpp @@ -37,6 +37,7 @@ $Id$ #include "astra/ParallelProjectionGeometry3D.h" #include "astra/ParallelVecProjectionGeometry3D.h" #include "astra/ConeVecProjectionGeometry3D.h" +#include "astra/CompositeGeometryManager.h" #include "astra/Logging.h" @@ -203,9 +204,16 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations) &volgeom, projgeom, m_iGPUIndex, m_iVoxelSuperSampling); } else { + +#if 1 + CCompositeGeometryManager cgm; + + cgm.doBP(m_pProjector, pReconMem, pSinoMem); +#else astraCudaBP(pReconMem->getData(), pSinoMem->getDataConst(), &volgeom, projgeom, m_iGPUIndex, m_iVoxelSuperSampling); +#endif } } diff --git a/src/CudaForwardProjectionAlgorithm3D.cpp b/src/CudaForwardProjectionAlgorithm3D.cpp index e57e077..209f5a5 100644 --- a/src/CudaForwardProjectionAlgorithm3D.cpp +++ b/src/CudaForwardProjectionAlgorithm3D.cpp @@ -40,6 +40,8 @@ $Id$ #include "astra/ParallelVecProjectionGeometry3D.h" #include "astra/ConeVecProjectionGeometry3D.h" +#include "astra/CompositeGeometryManager.h" + #include "astra/Logging.h" #include "../cuda/3d/astra3d.h" @@ -263,6 +265,12 @@ void CCudaForwardProjectionAlgorithm3D::run(int) // check initialized assert(m_bIsInitialized); +#if 1 + CCompositeGeometryManager cgm; + + cgm.doFP(m_pProjector, m_pVolume, m_pProjections); + +#else const CProjectionGeometry3D* projgeom = m_pProjections->getGeometry(); const CVolumeGeometry3D& volgeom = *m_pVolume->getGeometry(); @@ -294,6 +302,7 @@ void CCudaForwardProjectionAlgorithm3D::run(int) astraCudaFP(m_pVolume->getDataConst(), m_pProjections->getData(), &volgeom, projgeom, m_iGPUIndex, m_iDetectorSuperSampling, projKernel); +#endif } diff --git a/src/GeometryUtil3D.cpp b/src/GeometryUtil3D.cpp index 52dd5a9..c6bfd8b 100644 --- a/src/GeometryUtil3D.cpp +++ b/src/GeometryUtil3D.cpp @@ -28,8 +28,96 @@ $Id$ #include "astra/GeometryUtil3D.h" +#include + namespace astra { + +SConeProjection* genConeProjections(unsigned int iProjAngles, + unsigned int iProjU, + unsigned int iProjV, + double fOriginSourceDistance, + double fOriginDetectorDistance, + double fDetUSize, + double fDetVSize, + const float *pfAngles) +{ + SConeProjection base; + base.fSrcX = 0.0f; + base.fSrcY = -fOriginSourceDistance; + base.fSrcZ = 0.0f; + + base.fDetSX = iProjU * fDetUSize * -0.5f; + base.fDetSY = fOriginDetectorDistance; + base.fDetSZ = iProjV * fDetVSize * -0.5f; + + base.fDetUX = fDetUSize; + base.fDetUY = 0.0f; + base.fDetUZ = 0.0f; + + base.fDetVX = 0.0f; + base.fDetVY = 0.0f; + base.fDetVZ = fDetVSize; + + SConeProjection* p = new SConeProjection[iProjAngles]; + +#define ROTATE0(name,i,alpha) do { p[i].f##name##X = base.f##name##X * cos(alpha) - base.f##name##Y * sin(alpha); p[i].f##name##Y = base.f##name##X * sin(alpha) + base.f##name##Y * cos(alpha); p[i].f##name##Z = base.f##name##Z; } while(0) + + for (unsigned int i = 0; i < iProjAngles; ++i) { + ROTATE0(Src, i, pfAngles[i]); + ROTATE0(DetS, i, pfAngles[i]); + ROTATE0(DetU, i, pfAngles[i]); + ROTATE0(DetV, i, pfAngles[i]); + } + +#undef ROTATE0 + + return p; +} + +SPar3DProjection* genPar3DProjections(unsigned int iProjAngles, + unsigned int iProjU, + unsigned int iProjV, + double fDetUSize, + double fDetVSize, + const float *pfAngles) +{ + SPar3DProjection base; + base.fRayX = 0.0f; + base.fRayY = 1.0f; + base.fRayZ = 0.0f; + + base.fDetSX = iProjU * fDetUSize * -0.5f; + base.fDetSY = 0.0f; + base.fDetSZ = iProjV * fDetVSize * -0.5f; + + base.fDetUX = fDetUSize; + base.fDetUY = 0.0f; + base.fDetUZ = 0.0f; + + base.fDetVX = 0.0f; + base.fDetVY = 0.0f; + base.fDetVZ = fDetVSize; + + SPar3DProjection* p = new SPar3DProjection[iProjAngles]; + +#define ROTATE0(name,i,alpha) do { p[i].f##name##X = base.f##name##X * cos(alpha) - base.f##name##Y * sin(alpha); p[i].f##name##Y = base.f##name##X * sin(alpha) + base.f##name##Y * cos(alpha); p[i].f##name##Z = base.f##name##Z; } while(0) + + for (unsigned int i = 0; i < iProjAngles; ++i) { + ROTATE0(Ray, i, pfAngles[i]); + ROTATE0(DetS, i, pfAngles[i]); + ROTATE0(DetU, i, pfAngles[i]); + ROTATE0(DetV, i, pfAngles[i]); + } + +#undef ROTATE0 + + return p; +} + + + + // (See declaration in header for (mathematical) description of these functions) @@ -72,4 +160,88 @@ void computeBP_UV_Coeffs(const SConeProjection& proj, double &fUX, double &fUY, } +// TODO: Handle cases of rays parallel to coordinate planes + +void backprojectPointX(const SPar3DProjection& proj, double fU, double fV, + double fX, double &fY, double &fZ) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - px) / proj.fRayX; + + fY = py + a * proj.fRayY; + fZ = pz + a * proj.fRayZ; +} + +void backprojectPointY(const SPar3DProjection& proj, double fU, double fV, + double fY, double &fX, double &fZ) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - py) / proj.fRayY; + + fX = px + a * proj.fRayX; + fZ = pz + a * proj.fRayZ; + +} + +void backprojectPointZ(const SPar3DProjection& proj, double fU, double fV, + double fZ, double &fX, double &fY) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - pz) / proj.fRayZ; + + fX = px + a * proj.fRayX; + fY = py + a * proj.fRayY; +} + + + +void backprojectPointX(const SConeProjection& proj, double fU, double fV, + double fX, double &fY, double &fZ) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - proj.fSrcX) / (px - proj.fSrcX); + + fY = proj.fSrcY + a * (py - proj.fSrcY); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); +} + +void backprojectPointY(const SConeProjection& proj, double fU, double fV, + double fY, double &fX, double &fZ) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - proj.fSrcY) / (py - proj.fSrcY); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); +} + +void backprojectPointZ(const SConeProjection& proj, double fU, double fV, + double fZ, double &fX, double &fY) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - proj.fSrcZ) / (pz - proj.fSrcZ); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fY = proj.fSrcY + a * (py - proj.fSrcY); +} + + } diff --git a/src/ParallelProjectionGeometry3D.cpp b/src/ParallelProjectionGeometry3D.cpp index 1c87157..7b64fd9 100644 --- a/src/ParallelProjectionGeometry3D.cpp +++ b/src/ParallelProjectionGeometry3D.cpp @@ -27,8 +27,10 @@ $Id$ */ #include "astra/ParallelProjectionGeometry3D.h" -#include +#include "astra/GeometryUtil3D.h" + +#include #include using namespace std; @@ -185,9 +187,9 @@ CVector3D CParallelProjectionGeometry3D::getProjectionDirection(int _iProjection return CVector3D(fDirX, fDirY, fDirZ); } -void CParallelProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ, +void CParallelProjectionGeometry3D::projectPoint(double fX, double fY, double fZ, int iAngleIndex, - float32 &fU, float32 &fV) const + double &fU, double &fV) const { ASTRA_ASSERT(iAngleIndex >= 0); ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); @@ -214,6 +216,79 @@ CParallelProjectionGeometry2D * CParallelProjectionGeometry3D::createProjectionG return pOutput; } +void CParallelProjectionGeometry3D::backprojectPointX(int iAngleIndex, double fU, double fV, + double fX, double &fY, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SPar3DProjection *projs = genPar3DProjections(1, m_iDetectorColCount, m_iDetectorRowCount, + m_fDetectorSpacingX, m_fDetectorSpacingY, + &m_pfProjectionAngles[iAngleIndex]); + + SPar3DProjection &proj = projs[0]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - px) / proj.fRayX; + + fY = py + a * proj.fRayY; + fZ = pz + a * proj.fRayZ; + + delete[] projs; +} + +void CParallelProjectionGeometry3D::backprojectPointY(int iAngleIndex, double fU, double fV, + double fY, double &fX, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SPar3DProjection *projs = genPar3DProjections(1, m_iDetectorColCount, m_iDetectorRowCount, + m_fDetectorSpacingX, m_fDetectorSpacingY, + &m_pfProjectionAngles[iAngleIndex]); + + SPar3DProjection &proj = projs[0]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - py) / proj.fRayY; + + fX = px + a * proj.fRayX; + fZ = pz + a * proj.fRayZ; + + delete[] projs; +} + +void CParallelProjectionGeometry3D::backprojectPointZ(int iAngleIndex, double fU, double fV, + double fZ, double &fX, double &fY) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SPar3DProjection *projs = genPar3DProjections(1, m_iDetectorColCount, m_iDetectorRowCount, + m_fDetectorSpacingX, m_fDetectorSpacingY, + &m_pfProjectionAngles[iAngleIndex]); + + SPar3DProjection &proj = projs[0]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - pz) / proj.fRayZ; + + fX = px + a * proj.fRayX; + fY = py + a * proj.fRayY; + + delete[] projs; +} + + //---------------------------------------------------------------------------------------- } // end namespace astra diff --git a/src/ParallelVecProjectionGeometry3D.cpp b/src/ParallelVecProjectionGeometry3D.cpp index ffad6d0..d04400b 100644 --- a/src/ParallelVecProjectionGeometry3D.cpp +++ b/src/ParallelVecProjectionGeometry3D.cpp @@ -239,9 +239,9 @@ 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 +void CParallelVecProjectionGeometry3D::projectPoint(double fX, double fY, double fZ, + int iAngleIndex, + double &fU, double &fV) const { ASTRA_ASSERT(iAngleIndex >= 0); ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); @@ -258,6 +258,61 @@ void CParallelVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, floa } +void CParallelVecProjectionGeometry3D::backprojectPointX(int iAngleIndex, double fU, double fV, + double fX, double &fY, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SPar3DProjection &proj = m_pProjectionAngles[iAngleIndex]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - px) / proj.fRayX; + + fY = py + a * proj.fRayY; + fZ = pz + a * proj.fRayZ; +} + +void CParallelVecProjectionGeometry3D::backprojectPointY(int iAngleIndex, double fU, double fV, + double fY, double &fX, double &fZ) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SPar3DProjection &proj = m_pProjectionAngles[iAngleIndex]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - py) / proj.fRayY; + + fX = px + a * proj.fRayX; + fZ = pz + a * proj.fRayZ; +} + +void CParallelVecProjectionGeometry3D::backprojectPointZ(int iAngleIndex, double fU, double fV, + double fZ, double &fX, double &fY) const +{ + ASTRA_ASSERT(iAngleIndex >= 0); + ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount); + + SPar3DProjection &proj = m_pProjectionAngles[iAngleIndex]; + + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - pz) / proj.fRayZ; + + fX = px + a * proj.fRayX; + fY = py + a * proj.fRayY; +} + + //---------------------------------------------------------------------------------------- bool CParallelVecProjectionGeometry3D::_check() -- cgit v1.2.3 From 81e7385c110a6210d0f9bc402df522301ec162f6 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 4 Dec 2015 15:14:19 +0100 Subject: Add utility functions for creating FP/BP JobLists --- src/CompositeGeometryManager.cpp | 113 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 111 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index fc8bc2e..9be4797 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -102,7 +102,6 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div i_out != splitOutput.end(); ++i_out) { boost::shared_ptr outputPart = *i_out; - split[outputPart.get()] = TJobList(); SJob newjob; newjob.pOutput = outputPart; @@ -319,7 +318,7 @@ static size_t computeVerticalSplit(size_t maxBlock, int div, size_t sliceCount) // (This can't make the blocks larger) blockSize = ceildiv(sliceCount, blockCount); - ASTRA_DEBUG("%ld %ld -> %ld * %ld\n", sliceCount, maxBlock, blockCount, blockSize); + ASTRA_DEBUG("%ld %ld -> %ld * %ld", sliceCount, maxBlock, blockCount, blockSize); assert(blockSize <= maxBlock); assert((divCount * divCount > sliceCount) || (blockCount % div) == 0); @@ -725,6 +724,116 @@ bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeDat return doJobs(L); } +bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector& volData, const std::vector& projData) +{ + ASTRA_DEBUG("CCompositeGeometryManager::doFP, multi-volume"); + + std::vector::const_iterator i; + std::vector > inputs; + + for (i = volData.begin(); i != volData.end(); ++i) { + CVolumePart *input = new CVolumePart(); + input->pData = *i; + input->subX = 0; + input->subY = 0; + input->subZ = 0; + input->pGeom = (*i)->getGeometry()->clone(); + + inputs.push_back(boost::shared_ptr(input)); + } + + std::vector::const_iterator j; + std::vector > outputs; + + for (j = projData.begin(); j != projData.end(); ++j) { + CProjectionPart *output = new CProjectionPart(); + output->pData = *j; + output->subX = 0; + output->subY = 0; + output->subZ = 0; + output->pGeom = (*j)->getGeometry()->clone(); + + outputs.push_back(boost::shared_ptr(output)); + } + + std::vector >::iterator i2; + std::vector >::iterator j2; + TJobList L; + + for (i2 = outputs.begin(); i2 != outputs.end(); ++i2) { + SJob FP; + FP.eMode = SJob::MODE_SET; + for (j2 = inputs.begin(); j2 != inputs.end(); ++j2) { + FP.pInput = *j2; + FP.pOutput = *i2; + FP.pProjector = pProjector; + FP.eType = SJob::JOB_FP; + L.push_back(FP); + + // Set first, add rest + FP.eMode = SJob::MODE_ADD; + } + } + + return doJobs(L); +} + +bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, const std::vector& volData, const std::vector& projData) +{ + ASTRA_DEBUG("CCompositeGeometryManager::doBP, multi-volume"); + + + std::vector::const_iterator i; + std::vector > outputs; + + for (i = volData.begin(); i != volData.end(); ++i) { + CVolumePart *output = new CVolumePart(); + output->pData = *i; + output->subX = 0; + output->subY = 0; + output->subZ = 0; + output->pGeom = (*i)->getGeometry()->clone(); + + outputs.push_back(boost::shared_ptr(output)); + } + + std::vector::const_iterator j; + std::vector > inputs; + + for (j = projData.begin(); j != projData.end(); ++j) { + CProjectionPart *input = new CProjectionPart(); + input->pData = *j; + input->subX = 0; + input->subY = 0; + input->subZ = 0; + input->pGeom = (*j)->getGeometry()->clone(); + + inputs.push_back(boost::shared_ptr(input)); + } + + std::vector >::iterator i2; + std::vector >::iterator j2; + TJobList L; + + for (i2 = outputs.begin(); i2 != outputs.end(); ++i2) { + SJob BP; + BP.eMode = SJob::MODE_SET; + for (j2 = inputs.begin(); j2 != inputs.end(); ++j2) { + BP.pInput = *j2; + BP.pOutput = *i2; + BP.pProjector = pProjector; + BP.eType = SJob::JOB_BP; + L.push_back(BP); + + // Set first, add rest + BP.eMode = SJob::MODE_ADD; + } + } + + return doJobs(L); +} + + bool CCompositeGeometryManager::doJobs(TJobList &jobs) -- cgit v1.2.3