From 56809b0359af7e9108adeb1fd21823a225edf6fa Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 20 Jan 2016 18:08:59 +0100 Subject: Remove dependency of libastra on libpython by refactoring PluginAlgorithm --- src/PluginAlgorithm.cpp | 367 +----------------------------------------------- 1 file changed, 1 insertion(+), 366 deletions(-) (limited to 'src') diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp index 9fc511a..1bcfbdb 100644 --- a/src/PluginAlgorithm.cpp +++ b/src/PluginAlgorithm.cpp @@ -26,376 +26,11 @@ along with the ASTRA Toolbox. If not, see . $Id$ */ -#ifdef ASTRA_PYTHON - #include "astra/PluginAlgorithm.h" -#include "astra/Logging.h" -#include "astra/Utilities.h" -#include -#include -#include -#include -#include - -#include -#include "bytesobject.h" namespace astra { +CPluginAlgorithmFactory *CPluginAlgorithmFactory::m_factory = 0; - -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; - 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(){ - 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){ - logPythonError(); - return false; - } - m_bIsInitialized = true; - Py_DECREF(retVal); - return m_bIsInitialized; -} - -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(); - }else{ - Py_DECREF(retVal); - } - PyGILState_Release(state); } -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(){ - if(!Py_IsInitialized()){ - Py_Initialize(); - PyEval_InitThreads(); - } -#ifndef _MSC_VER - if(astra::running_in_matlab) fixLapackLoading(); -#endif - pluginDict = PyDict_New(); - inspect = PyImport_ImportModule("inspect"); - six = PyImport_ImportModule("six"); -} - -CPluginAlgorithmFactory::~CPluginAlgorithmFactory(){ - if(pluginDict!=NULL){ - Py_DECREF(pluginDict); - } - if(inspect!=NULL) Py_DECREF(inspect); - if(six!=NULL) Py_DECREF(six); -} - -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){ - logPythonError(); - return NULL; - } - PyObject *submod = pyclass; - for(unsigned int i=1;i CPluginAlgorithmFactory::getRegisteredMap(){ - std::map ret; - PyObject *key, *value; - Py_ssize_t pos = 0; - while (PyDict_Next(pluginDict, &pos, &key, &value)) { - 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; -} - -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 = ""; - 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); - if(retVal!=NULL){ - 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(); - } - } - if(PyBytes_Check(className)){ - Py_DECREF(pyclass); - } - return ret; -} - -DEFINE_SINGLETON(CPluginAlgorithmFactory); - -#if PY_MAJOR_VERSION >= 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 vec; - boost::split(vec, str, boost::is_any_of(",")); - PyObject *veclist = PyList_New(vec.size()); - for(unsigned int i=0;i nodes = node.getNodes(); - std::list::iterator it = nodes.begin(); - while(it!=nodes.end()){ - XMLNode subnode = *it; - if(subnode.getName()=="Option"){ - 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{ - 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 -- cgit v1.2.3 From 838cfae58d825fb8915dc7d3c974d96e6a4f981c Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Fri, 12 Feb 2016 16:27:08 +0100 Subject: Also split volumes in X/Y directions to respect CUDA limits --- src/CompositeGeometryManager.cpp | 261 +++++++++++++++++++++++++++++++++++---- 1 file changed, 240 insertions(+), 21 deletions(-) (limited to 'src') diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index 96b28e9..1dd12ea 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -51,8 +51,11 @@ along with the ASTRA Toolbox. If not, see . #include #endif + namespace astra { +static const size_t MAX_BLOCK_DIM = 4096; + SGPUParams* CCompositeGeometryManager::s_params = 0; @@ -111,7 +114,20 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div // b. split input part // c. create jobs for new (input,output) subparts - TPartList splitOutput = pOutput->split(maxSize/3, div); + TPartList splitOutput; + pOutput->splitZ(splitOutput, maxSize/3, MAX_BLOCK_DIM, div); + TPartList splitOutput2; + for (TPartList::iterator i_out = splitOutput.begin(); i_out != splitOutput.end(); ++i_out) { + boost::shared_ptr outputPart = *i_out; + outputPart.get()->splitX(splitOutput2, maxSize/3, MAX_BLOCK_DIM, 1); + } + splitOutput.clear(); + for (TPartList::iterator i_out = splitOutput2.begin(); i_out != splitOutput2.end(); ++i_out) { + boost::shared_ptr outputPart = *i_out; + outputPart.get()->splitY(splitOutput, maxSize/3, MAX_BLOCK_DIM, 1); + } + splitOutput2.clear(); + for (TJobList::const_iterator j = L.begin(); j != L.end(); ++j) { @@ -139,8 +155,21 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div size_t remainingSize = ( maxSize - outputPart->getSize() ) / 2; - TPartList splitInput = input->split(remainingSize, 1); + TPartList splitInput; + input->splitZ(splitInput, remainingSize, MAX_BLOCK_DIM, 1); delete input; + TPartList splitInput2; + for (TPartList::iterator i_in = splitInput.begin(); i_in != splitInput.end(); ++i_in) { + boost::shared_ptr inputPart = *i_in; + inputPart.get()->splitX(splitInput2, maxSize/3, MAX_BLOCK_DIM, 1); + } + splitInput.clear(); + for (TPartList::iterator i_in = splitInput2.begin(); i_in != splitInput2.end(); ++i_in) { + boost::shared_ptr inputPart = *i_in; + inputPart.get()->splitY(splitInput, maxSize/3, MAX_BLOCK_DIM, 1); + } + splitInput2.clear(); + ASTRA_DEBUG("Input split into %d parts", splitInput.size()); for (TPartList::iterator i_in = splitInput.begin(); @@ -327,7 +356,7 @@ 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) +static size_t computeLinearSplit(size_t maxBlock, int div, size_t sliceCount) { size_t blockSize = maxBlock; size_t blockCount = ceildiv(sliceCount, blockSize); @@ -410,7 +439,17 @@ SPar3DProjection* getProjectionVectors(const CParallelVecProjectionGeometry3D* p template -static void translateProjectionVectors(V* pProjs, int count, double dv) +static void translateProjectionVectorsU(V* pProjs, int count, double du) +{ + for (int i = 0; i < count; ++i) { + pProjs[i].fDetSX += du * pProjs[i].fDetUX; + pProjs[i].fDetSY += du * pProjs[i].fDetUY; + pProjs[i].fDetSZ += du * pProjs[i].fDetUZ; + } +} + +template +static void translateProjectionVectorsV(V* pProjs, int count, double dv) { for (int i = 0; i < count; ++i) { pProjs[i].fDetSX += dv * pProjs[i].fDetVX; @@ -420,8 +459,58 @@ static void translateProjectionVectors(V* pProjs, int count, double dv) } +static CProjectionGeometry3D* getSubProjectionGeometryU(const CProjectionGeometry3D* pProjGeom, int u, 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); + } + + translateProjectionVectorsU(pConeProjs, pProjGeom->getProjectionCount(), u); + + CProjectionGeometry3D* ret = new CConeVecProjectionGeometry3D(pProjGeom->getProjectionCount(), + pProjGeom->getDetectorRowCount(), + size, + pConeProjs); + + + delete[] pConeProjs; + return ret; + } else { + assert(par3dgeom || parvec3dgeom); + SPar3DProjection* pParProjs; + if (par3dgeom) { + pParProjs = getProjectionVectors(par3dgeom); + } else { + pParProjs = getProjectionVectors(parvec3dgeom); + } + + translateProjectionVectorsU(pParProjs, pProjGeom->getProjectionCount(), u); + + CProjectionGeometry3D* ret = new CParallelVecProjectionGeometry3D(pProjGeom->getProjectionCount(), + pProjGeom->getDetectorRowCount(), + size, + pParProjs); + + delete[] pParProjs; + return ret; + } + +} + + -static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry3D* pProjGeom, int v, int size) +static CProjectionGeometry3D* getSubProjectionGeometryV(const CProjectionGeometry3D* pProjGeom, int v, int size) { // First convert to vectors, then translate, then convert into new object @@ -438,7 +527,7 @@ static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry pConeProjs = getProjectionVectors(conevec3dgeom); } - translateProjectionVectors(pConeProjs, pProjGeom->getProjectionCount(), v); + translateProjectionVectorsV(pConeProjs, pProjGeom->getProjectionCount(), v); CProjectionGeometry3D* ret = new CConeVecProjectionGeometry3D(pProjGeom->getProjectionCount(), size, @@ -457,7 +546,7 @@ static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry pParProjs = getProjectionVectors(parvec3dgeom); } - translateProjectionVectors(pParProjs, pProjGeom->getProjectionCount(), v); + translateProjectionVectorsV(pParProjs, pProjGeom->getProjectionCount(), v); CProjectionGeometry3D* ret = new CParallelVecProjectionGeometry3D(pProjGeom->getProjectionCount(), size, @@ -476,17 +565,110 @@ static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry // - 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) +void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::TPartList& out, size_t maxSize, size_t maxDim, int div) +{ + 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->getGridSliceCount()) * pGeom->getGridRowCount(); + int sliceCount = pGeom->getGridColCount(); + size_t m = std::min(maxSize / sliceSize, maxDim); + size_t blockSize = computeLinearSplit(m, div, sliceCount); + + int rem = sliceCount % blockSize; + + ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); + + for (int x = -(rem / 2); x < sliceCount; x += blockSize) { + int newsubX = x; + if (newsubX < 0) newsubX = 0; + int endX = x + blockSize; + if (endX > sliceCount) endX = sliceCount; + int size = endX - newsubX; + + CVolumePart *sub = new CVolumePart(); + sub->subX = this->subX + newsubX; + sub->subY = this->subY; + sub->subZ = this->subZ; + + ASTRA_DEBUG("VolumePart split %d %d %d -> %p", sub->subX, sub->subY, sub->subZ, (void*)sub); + + double shift = pGeom->getPixelLengthX() * newsubX; + + sub->pData = pData; + sub->pGeom = new CVolumeGeometry3D(size, + pGeom->getGridRowCount(), + pGeom->getGridSliceCount(), + pGeom->getWindowMinX() + shift, + pGeom->getWindowMinY(), + pGeom->getWindowMinZ(), + pGeom->getWindowMinX() + shift + size * pGeom->getPixelLengthX(), + pGeom->getWindowMaxY(), + pGeom->getWindowMaxZ()); + + out.push_back(boost::shared_ptr(sub)); + } + } +} + +void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::TPartList& out, size_t maxSize, size_t maxDim, 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->getGridSliceCount(); + int sliceCount = pGeom->getGridRowCount(); + size_t m = std::min(maxSize / sliceSize, maxDim); + size_t blockSize = computeLinearSplit(m, div, sliceCount); + + int rem = sliceCount % blockSize; + + ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize); + + for (int y = -(rem / 2); y < sliceCount; y += blockSize) { + int newsubY = y; + if (newsubY < 0) newsubY = 0; + int endY = y + blockSize; + if (endY > sliceCount) endY = sliceCount; + int size = endY - newsubY; + + CVolumePart *sub = new CVolumePart(); + sub->subX = this->subX; + sub->subY = this->subY + newsubY; + sub->subZ = this->subZ; + + ASTRA_DEBUG("VolumePart split %d %d %d -> %p", sub->subX, sub->subY, sub->subZ, (void*)sub); + + double shift = pGeom->getPixelLengthY() * newsubY; + + sub->pData = pData; + sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(), + size, + pGeom->getGridSliceCount(), + pGeom->getWindowMinX(), + pGeom->getWindowMinY() + shift, + pGeom->getWindowMinZ(), + pGeom->getWindowMaxX(), + pGeom->getWindowMinY() + shift + size * pGeom->getPixelLengthY(), + pGeom->getWindowMaxZ()); + out.push_back(boost::shared_ptr(sub)); + } + } +} + +void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::TPartList& out, size_t maxSize, size_t maxDim, int div) +{ 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); + size_t m = std::min(maxSize / sliceSize, maxDim); + size_t blockSize = computeLinearSplit(m, div, sliceCount); int rem = sliceCount % blockSize; @@ -519,11 +701,9 @@ CCompositeGeometryManager::TPartList CCompositeGeometryManager::CVolumePart::spl pGeom->getWindowMaxY(), pGeom->getWindowMinZ() + shift + size * pGeom->getPixelLengthZ()); - ret.push_back(boost::shared_ptr(sub)); + out.push_back(boost::shared_ptr(sub)); } } - - return ret; } CCompositeGeometryManager::CVolumePart* CCompositeGeometryManager::CVolumePart::clone() const @@ -630,7 +810,7 @@ CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::re if (_vmin == _vmax) { sub->pGeom = 0; } else { - sub->pGeom = getSubProjectionGeometry(pGeom, _vmin, _vmax - _vmin); + sub->pGeom = getSubProjectionGeometryV(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); @@ -639,17 +819,58 @@ CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::re } -CCompositeGeometryManager::TPartList CCompositeGeometryManager::CProjectionPart::split(size_t maxSize, int div) +void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManager::TPartList &out, size_t maxSize, size_t maxDim, int div) +{ + 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->getDetectorRowCount()) * pGeom->getProjectionCount(); + int sliceCount = pGeom->getDetectorColCount(); + size_t m = std::min(maxSize / sliceSize, maxDim); + size_t blockSize = computeLinearSplit(m, div, sliceCount); + + int rem = sliceCount % blockSize; + + for (int x = -(rem / 2); x < sliceCount; x += blockSize) { + int newsubX = x; + if (newsubX < 0) newsubX = 0; + int endX = x + blockSize; + if (endX > sliceCount) endX = sliceCount; + int size = endX - newsubX; + + CProjectionPart *sub = new CProjectionPart(); + sub->subX = this->subX + newsubX; + sub->subY = this->subY; + sub->subZ = this->subZ; + + ASTRA_DEBUG("ProjectionPart split %d %d %d -> %p", sub->subX, sub->subY, sub->subZ, (void*)sub); + + sub->pData = pData; + + sub->pGeom = getSubProjectionGeometryU(pGeom, newsubX, size); + + out.push_back(boost::shared_ptr(sub)); + } + } +} + +void CCompositeGeometryManager::CProjectionPart::splitY(CCompositeGeometryManager::TPartList &out, size_t maxSize, size_t maxDim, int div) { - TPartList ret; + // TODO + out.push_back(boost::shared_ptr(clone())); +} +void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManager::TPartList &out, size_t maxSize, size_t maxDim, int div) +{ 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); + size_t m = std::min(maxSize / sliceSize, maxDim); + size_t blockSize = computeLinearSplit(m, div, sliceCount); int rem = sliceCount % blockSize; @@ -669,14 +890,12 @@ CCompositeGeometryManager::TPartList CCompositeGeometryManager::CProjectionPart: sub->pData = pData; - sub->pGeom = getSubProjectionGeometry(pGeom, newsubZ, size); + sub->pGeom = getSubProjectionGeometryV(pGeom, newsubZ, size); - ret.push_back(boost::shared_ptr(sub)); + out.push_back(boost::shared_ptr(sub)); } } - return ret; - } CCompositeGeometryManager::CProjectionPart* CCompositeGeometryManager::CProjectionPart::clone() const -- cgit v1.2.3 From e9fad320817cd8ab84f7ef81940fda63f975551e Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 15 Feb 2016 13:37:22 +0100 Subject: Skip some unnecessary splitting --- src/CompositeGeometryManager.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index 1dd12ea..1991731 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -45,6 +45,7 @@ along with the ASTRA Toolbox. If not, see . #include #include +#include #ifndef USE_PTHREADS #include @@ -115,18 +116,20 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div // c. create jobs for new (input,output) subparts TPartList splitOutput; - pOutput->splitZ(splitOutput, maxSize/3, MAX_BLOCK_DIM, div); + pOutput->splitZ(splitOutput, maxSize/3, SIZE_MAX, div); +#if 0 TPartList splitOutput2; for (TPartList::iterator i_out = splitOutput.begin(); i_out != splitOutput.end(); ++i_out) { boost::shared_ptr outputPart = *i_out; - outputPart.get()->splitX(splitOutput2, maxSize/3, MAX_BLOCK_DIM, 1); + outputPart.get()->splitX(splitOutput2, SIZE_MAX, SIZE_MAX, 1); } splitOutput.clear(); for (TPartList::iterator i_out = splitOutput2.begin(); i_out != splitOutput2.end(); ++i_out) { boost::shared_ptr outputPart = *i_out; - outputPart.get()->splitY(splitOutput, maxSize/3, MAX_BLOCK_DIM, 1); + outputPart.get()->splitY(splitOutput, SIZE_MAX, SIZE_MAX, 1); } splitOutput2.clear(); +#endif for (TJobList::const_iterator j = L.begin(); j != L.end(); ++j) @@ -161,12 +164,12 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div TPartList splitInput2; for (TPartList::iterator i_in = splitInput.begin(); i_in != splitInput.end(); ++i_in) { boost::shared_ptr inputPart = *i_in; - inputPart.get()->splitX(splitInput2, maxSize/3, MAX_BLOCK_DIM, 1); + inputPart.get()->splitX(splitInput2, SIZE_MAX, MAX_BLOCK_DIM, 1); } splitInput.clear(); for (TPartList::iterator i_in = splitInput2.begin(); i_in != splitInput2.end(); ++i_in) { boost::shared_ptr inputPart = *i_in; - inputPart.get()->splitY(splitInput, maxSize/3, MAX_BLOCK_DIM, 1); + inputPart.get()->splitY(splitInput, SIZE_MAX, MAX_BLOCK_DIM, 1); } splitInput2.clear(); -- cgit v1.2.3 From 447e7acfb0c220f66d5fe25f31b25c989d4ec1d7 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 15 Feb 2016 15:44:49 +0100 Subject: Avoid (unlikely) integer overflow --- src/CompositeGeometryManager.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index 1991731..cafc452 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -362,7 +362,11 @@ static size_t ceildiv(size_t a, size_t b) { static size_t computeLinearSplit(size_t maxBlock, int div, size_t sliceCount) { size_t blockSize = maxBlock; - size_t blockCount = ceildiv(sliceCount, blockSize); + size_t blockCount; + if (sliceCount <= blockSize) + blockCount = 1; + else + blockCount = ceildiv(sliceCount, blockSize); // Increase number of blocks to be divisible by div size_t divCount = div * ceildiv(blockCount, div); -- cgit v1.2.3 From bc2e4018054f494fcba01e6a27a63e151bf1e9a4 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 16 Feb 2016 15:15:13 +0100 Subject: Refactor AstraObjectManager to add an AstraIndexManager The new AstraIndexManager can be used to obtain information about objects without knowing their type. --- src/AstraObjectManager.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) (limited to 'src') diff --git a/src/AstraObjectManager.cpp b/src/AstraObjectManager.cpp index c49f273..46eae4b 100644 --- a/src/AstraObjectManager.cpp +++ b/src/AstraObjectManager.cpp @@ -31,13 +31,13 @@ $Id$ namespace astra { -int CAstraIndexManager::m_iPreviousIndex = 0; - -DEFINE_SINGLETON(CAstraObjectManager); -DEFINE_SINGLETON(CAstraObjectManager); -DEFINE_SINGLETON(CAstraObjectManager); -DEFINE_SINGLETON(CAstraObjectManager); -DEFINE_SINGLETON(CAstraObjectManager); -DEFINE_SINGLETON(CAstraObjectManager); +DEFINE_SINGLETON(CProjector2DManager); +DEFINE_SINGLETON(CProjector3DManager); +DEFINE_SINGLETON(CData2DManager); +DEFINE_SINGLETON(CData3DManager); +DEFINE_SINGLETON(CAlgorithmManager); +DEFINE_SINGLETON(CMatrixManager); + +DEFINE_SINGLETON(CAstraIndexManager); } // end namespace -- cgit v1.2.3 From 3743fdc534b39958c105f4124ad1130d3e8b042a Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 16 Feb 2016 17:53:24 +0100 Subject: Query max texture size instead of hardcoding it --- src/CompositeGeometryManager.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) (limited to 'src') diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index cafc452..c9cbaaa 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -55,9 +55,6 @@ along with the ASTRA Toolbox. If not, see . namespace astra { -static const size_t MAX_BLOCK_DIM = 4096; - - SGPUParams* CCompositeGeometryManager::s_params = 0; CCompositeGeometryManager::CCompositeGeometryManager() @@ -102,6 +99,9 @@ CCompositeGeometryManager::CCompositeGeometryManager() bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div, TJobSet &split) { + int maxBlockDim = astraCUDA3d::maxBlockDimension(); + ASTRA_DEBUG("Found max block dim %d", maxBlockDim); + split.clear(); for (TJobSet::const_iterator i = jobs.begin(); i != jobs.end(); ++i) @@ -159,17 +159,17 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div size_t remainingSize = ( maxSize - outputPart->getSize() ) / 2; TPartList splitInput; - input->splitZ(splitInput, remainingSize, MAX_BLOCK_DIM, 1); + input->splitZ(splitInput, remainingSize, maxBlockDim, 1); delete input; TPartList splitInput2; for (TPartList::iterator i_in = splitInput.begin(); i_in != splitInput.end(); ++i_in) { boost::shared_ptr inputPart = *i_in; - inputPart.get()->splitX(splitInput2, SIZE_MAX, MAX_BLOCK_DIM, 1); + inputPart.get()->splitX(splitInput2, SIZE_MAX, maxBlockDim, 1); } splitInput.clear(); for (TPartList::iterator i_in = splitInput2.begin(); i_in != splitInput2.end(); ++i_in) { boost::shared_ptr inputPart = *i_in; - inputPart.get()->splitY(splitInput, SIZE_MAX, MAX_BLOCK_DIM, 1); + inputPart.get()->splitY(splitInput, SIZE_MAX, maxBlockDim, 1); } splitInput2.clear(); -- cgit v1.2.3 From 7b8a508f0bc7a8a02766b15fa094dfd18c1b0525 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 1 Mar 2016 15:05:53 +0100 Subject: Fix build --- src/CompositeGeometryManager.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) (limited to 'src') diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index c9cbaaa..084ba8c 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -45,7 +45,7 @@ along with the ASTRA Toolbox. If not, see . #include #include -#include +#include #ifndef USE_PTHREADS #include @@ -116,17 +116,17 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div // c. create jobs for new (input,output) subparts TPartList splitOutput; - pOutput->splitZ(splitOutput, maxSize/3, SIZE_MAX, div); + pOutput->splitZ(splitOutput, maxSize/3, UINT_MAX, div); #if 0 TPartList splitOutput2; for (TPartList::iterator i_out = splitOutput.begin(); i_out != splitOutput.end(); ++i_out) { boost::shared_ptr outputPart = *i_out; - outputPart.get()->splitX(splitOutput2, SIZE_MAX, SIZE_MAX, 1); + outputPart.get()->splitX(splitOutput2, UINT_MAX, UINT_MAX, 1); } splitOutput.clear(); for (TPartList::iterator i_out = splitOutput2.begin(); i_out != splitOutput2.end(); ++i_out) { boost::shared_ptr outputPart = *i_out; - outputPart.get()->splitY(splitOutput, SIZE_MAX, SIZE_MAX, 1); + outputPart.get()->splitY(splitOutput, UINT_MAX, UINT_MAX, 1); } splitOutput2.clear(); #endif @@ -164,12 +164,12 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div TPartList splitInput2; for (TPartList::iterator i_in = splitInput.begin(); i_in != splitInput.end(); ++i_in) { boost::shared_ptr inputPart = *i_in; - inputPart.get()->splitX(splitInput2, SIZE_MAX, maxBlockDim, 1); + inputPart.get()->splitX(splitInput2, UINT_MAX, maxBlockDim, 1); } splitInput.clear(); for (TPartList::iterator i_in = splitInput2.begin(); i_in != splitInput2.end(); ++i_in) { boost::shared_ptr inputPart = *i_in; - inputPart.get()->splitY(splitInput, SIZE_MAX, maxBlockDim, 1); + inputPart.get()->splitY(splitInput, UINT_MAX, maxBlockDim, 1); } splitInput2.clear(); -- cgit v1.2.3 From 14bef5ea534e4aa4e6d0819e728d0a8d2b0b7925 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 1 Mar 2016 15:05:53 +0100 Subject: Fix build --- src/CompositeGeometryManager.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) (limited to 'src') diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index c9cbaaa..084ba8c 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -45,7 +45,7 @@ along with the ASTRA Toolbox. If not, see . #include #include -#include +#include #ifndef USE_PTHREADS #include @@ -116,17 +116,17 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div // c. create jobs for new (input,output) subparts TPartList splitOutput; - pOutput->splitZ(splitOutput, maxSize/3, SIZE_MAX, div); + pOutput->splitZ(splitOutput, maxSize/3, UINT_MAX, div); #if 0 TPartList splitOutput2; for (TPartList::iterator i_out = splitOutput.begin(); i_out != splitOutput.end(); ++i_out) { boost::shared_ptr outputPart = *i_out; - outputPart.get()->splitX(splitOutput2, SIZE_MAX, SIZE_MAX, 1); + outputPart.get()->splitX(splitOutput2, UINT_MAX, UINT_MAX, 1); } splitOutput.clear(); for (TPartList::iterator i_out = splitOutput2.begin(); i_out != splitOutput2.end(); ++i_out) { boost::shared_ptr outputPart = *i_out; - outputPart.get()->splitY(splitOutput, SIZE_MAX, SIZE_MAX, 1); + outputPart.get()->splitY(splitOutput, UINT_MAX, UINT_MAX, 1); } splitOutput2.clear(); #endif @@ -164,12 +164,12 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div TPartList splitInput2; for (TPartList::iterator i_in = splitInput.begin(); i_in != splitInput.end(); ++i_in) { boost::shared_ptr inputPart = *i_in; - inputPart.get()->splitX(splitInput2, SIZE_MAX, maxBlockDim, 1); + inputPart.get()->splitX(splitInput2, UINT_MAX, maxBlockDim, 1); } splitInput.clear(); for (TPartList::iterator i_in = splitInput2.begin(); i_in != splitInput2.end(); ++i_in) { boost::shared_ptr inputPart = *i_in; - inputPart.get()->splitY(splitInput, SIZE_MAX, maxBlockDim, 1); + inputPart.get()->splitY(splitInput, UINT_MAX, maxBlockDim, 1); } splitInput2.clear(); -- cgit v1.2.3 From 495903529d473a9968c1333d5a515e3b94732f0b Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 23 Mar 2016 15:29:43 +0100 Subject: Move CUDA algorithm initialization to its own function --- src/CudaReconstructionAlgorithm2D.cpp | 22 +++++++++++++++------- src/CudaSirtAlgorithm.cpp | 35 +++++++++++------------------------ 2 files changed, 26 insertions(+), 31 deletions(-) (limited to 'src') diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp index 5a1910c..2798434 100644 --- a/src/CudaReconstructionAlgorithm2D.cpp +++ b/src/CudaReconstructionAlgorithm2D.cpp @@ -328,6 +328,20 @@ bool CCudaReconstructionAlgorithm2D::setupGeometry() return true; } +//---------------------------------------------------------------------------------------- + +void CCudaReconstructionAlgorithm2D::initCUDAAlgorithm() +{ + bool ok; + + ok = setupGeometry(); + ASTRA_ASSERT(ok); + + ok = m_pAlgo->allocateBuffers(); + ASTRA_ASSERT(ok); +} + + //---------------------------------------------------------------------------------------- // Iterate void CCudaReconstructionAlgorithm2D::run(int _iNrIterations) @@ -339,13 +353,7 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations) const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); if (!m_bAlgoInit) { - - ok = setupGeometry(); - ASTRA_ASSERT(ok); - - ok = m_pAlgo->allocateBuffers(); - ASTRA_ASSERT(ok); - + initCUDAAlgorithm(); m_bAlgoInit = true; } diff --git a/src/CudaSirtAlgorithm.cpp b/src/CudaSirtAlgorithm.cpp index 33e381a..7beb30e 100644 --- a/src/CudaSirtAlgorithm.cpp +++ b/src/CudaSirtAlgorithm.cpp @@ -113,36 +113,23 @@ bool CCudaSirtAlgorithm::initialize(CProjector2D* _pProjector, } //---------------------------------------------------------------------------------------- -// Iterate -void CCudaSirtAlgorithm::run(int _iNrIterations) -{ - // check initialized - ASTRA_ASSERT(m_bIsInitialized); - if (!m_bAlgoInit) { - // We only override the initialisation step to copy the min/max masks +void CCudaSirtAlgorithm::initCUDAAlgorithm() +{ + CCudaReconstructionAlgorithm2D::initCUDAAlgorithm(); - bool ok = setupGeometry(); - ASTRA_ASSERT(ok); + astraCUDA::SIRT* pSirt = dynamic_cast(m_pAlgo); - ok = m_pAlgo->allocateBuffers(); + if (m_pMinMask || m_pMaxMask) { + const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); + const float *pfMinMaskData = 0; + const float *pfMaxMaskData = 0; + if (m_pMinMask) pfMinMaskData = m_pMinMask->getDataConst(); + if (m_pMaxMask) pfMaxMaskData = m_pMaxMask->getDataConst(); + bool ok = pSirt->uploadMinMaxMasks(pfMinMaskData, pfMaxMaskData, volgeom.getGridColCount()); ASTRA_ASSERT(ok); - - if (m_pMinMask || m_pMaxMask) { - const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); - astraCUDA::SIRT* pSirt = dynamic_cast(m_pAlgo); - const float *pfMinMaskData = 0; - const float *pfMaxMaskData = 0; - if (m_pMinMask) pfMinMaskData = m_pMinMask->getDataConst(); - if (m_pMaxMask) pfMaxMaskData = m_pMaxMask->getDataConst(); - ok = pSirt->uploadMinMaxMasks(pfMinMaskData, pfMaxMaskData, volgeom.getGridColCount()); - ASTRA_ASSERT(ok); - } - - m_bAlgoInit = true; } - CCudaReconstructionAlgorithm2D::run(_iNrIterations); } -- cgit v1.2.3 From f03ceb16d2dbde0c43e8c90683c5feafe01e5356 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 23 Mar 2016 15:30:47 +0100 Subject: Rename ART lambda option to Relaxation --- src/ArtAlgorithm.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/ArtAlgorithm.cpp b/src/ArtAlgorithm.cpp index b59bd93..526c263 100644 --- a/src/ArtAlgorithm.cpp +++ b/src/ArtAlgorithm.cpp @@ -156,8 +156,12 @@ bool CArtAlgorithm::initialize(const Config& _cfg) return false; } + // "Lambda" is replaced by the more descriptive "Relaxation" m_fLambda = _cfg.self.getOptionNumerical("Lambda", 1.0f); - CC.markOptionParsed("Lambda"); + m_fLambda = _cfg.self.getOptionNumerical("Relaxation", m_fLambda); + if (!_cfg.self.hasOption("Relaxation")) + CC.markOptionParsed("Lambda"); + CC.markOptionParsed("Relaxation"); // success m_bIsInitialized = _check(); @@ -232,7 +236,7 @@ map CArtAlgorithm::getInformation() { map res; res["RayOrder"] = getInformation("RayOrder"); - res["Lambda"] = getInformation("Lambda"); + res["Relaxation"] = getInformation("Relaxation"); return mergeMap(CReconstructionAlgorithm2D::getInformation(), res); }; @@ -240,7 +244,7 @@ map CArtAlgorithm::getInformation() // Information - Specific boost::any CArtAlgorithm::getInformation(std::string _sIdentifier) { - if (_sIdentifier == "Lambda") { return m_fLambda; } + if (_sIdentifier == "Relaxation") { return m_fLambda; } if (_sIdentifier == "RayOrder") { vector res; for (int i = 0; i < m_iRayCount; i++) { -- cgit v1.2.3 From 5edb35edc2c721b458334a65512b534912c2c542 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 23 Mar 2016 15:30:56 +0100 Subject: Add relaxation parameters to SIRT, SART --- src/CudaSartAlgorithm.cpp | 17 ++++++++++++++++- src/CudaSirtAlgorithm.cpp | 6 ++++++ src/SartAlgorithm.cpp | 8 +++++++- src/SirtAlgorithm.cpp | 11 +++++++++-- 4 files changed, 38 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/CudaSartAlgorithm.cpp b/src/CudaSartAlgorithm.cpp index d202847..bf97224 100644 --- a/src/CudaSartAlgorithm.cpp +++ b/src/CudaSartAlgorithm.cpp @@ -107,7 +107,8 @@ bool CCudaSartAlgorithm::initialize(const Config& _cfg) CC.markOptionParsed("ProjectionOrderList"); } - + m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f); + CC.markOptionParsed("Relaxation"); return true; } @@ -123,12 +124,26 @@ bool CCudaSartAlgorithm::initialize(CProjector2D* _pProjector, if (!m_bIsInitialized) return false; + m_fLambda = 1.0f; + m_pAlgo = new astraCUDA::SART(); m_bAlgoInit = false; return true; } +//---------------------------------------------------------------------------------------- + +void CCudaSartAlgorithm::initCUDAAlgorithm() +{ + CCudaReconstructionAlgorithm2D::initCUDAAlgorithm(); + + astraCUDA::SART* pSart = dynamic_cast(m_pAlgo); + + pSart->setRelaxation(m_fLambda); +} + + } // namespace astra diff --git a/src/CudaSirtAlgorithm.cpp b/src/CudaSirtAlgorithm.cpp index 7beb30e..c8dc677 100644 --- a/src/CudaSirtAlgorithm.cpp +++ b/src/CudaSirtAlgorithm.cpp @@ -50,6 +50,8 @@ CCudaSirtAlgorithm::CCudaSirtAlgorithm() m_pMinMask = 0; m_pMaxMask = 0; + + m_fLambda = 1.0f; } //---------------------------------------------------------------------------------------- @@ -86,6 +88,8 @@ bool CCudaSirtAlgorithm::initialize(const Config& _cfg) } CC.markOptionParsed("MaxMaskId"); + m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f); + CC.markOptionParsed("Relaxation"); m_pAlgo = new astraCUDA::SIRT(); m_bAlgoInit = false; @@ -108,6 +112,7 @@ bool CCudaSirtAlgorithm::initialize(CProjector2D* _pProjector, m_pAlgo = new astraCUDA::SIRT(); m_bAlgoInit = false; + m_fLambda = 1.0f; return true; } @@ -130,6 +135,7 @@ void CCudaSirtAlgorithm::initCUDAAlgorithm() ASTRA_ASSERT(ok); } + pSirt->setRelaxation(m_fLambda); } diff --git a/src/SartAlgorithm.cpp b/src/SartAlgorithm.cpp index 9346160..403f851 100644 --- a/src/SartAlgorithm.cpp +++ b/src/SartAlgorithm.cpp @@ -151,6 +151,9 @@ bool CSartAlgorithm::initialize(const Config& _cfg) CC.markOptionParsed("ProjectionOrderList"); } + m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f); + CC.markOptionParsed("Relaxation"); + // create data objects m_pTotalRayLength = new CFloat32ProjectionData2D(m_pProjector->getProjectionGeometry()); m_pTotalPixelWeight = new CFloat32VolumeData2D(m_pProjector->getVolumeGeometry()); @@ -246,6 +249,7 @@ map CSartAlgorithm::getInformation() { map res; res["ProjectionOrder"] = getInformation("ProjectionOrder"); + res["Relaxation"] = getInformation("Relaxation"); return mergeMap(CReconstructionAlgorithm2D::getInformation(), res); }; @@ -253,6 +257,8 @@ map CSartAlgorithm::getInformation() // Information - Specific boost::any CSartAlgorithm::getInformation(std::string _sIdentifier) { + if (_sIdentifier == "Relaxation") + return m_fLambda; if (_sIdentifier == "ProjectionOrder") { vector res; for (int i = 0; i < m_iProjectionCount; i++) { @@ -286,7 +292,7 @@ void CSartAlgorithm::run(int _iNrIterations) m_pProjector, SinogramMaskPolicy(m_pSinogramMask), // sinogram mask ReconstructionMaskPolicy(m_pReconstructionMask), // reconstruction mask - SIRTBPPolicy(m_pReconstruction, m_pDiffSinogram, m_pTotalPixelWeight, m_pTotalRayLength), // SIRT backprojection + SIRTBPPolicy(m_pReconstruction, m_pDiffSinogram, m_pTotalPixelWeight, m_pTotalRayLength, m_fLambda), // SIRT backprojection m_bUseSinogramMask, m_bUseReconstructionMask, true // options on/off ); diff --git a/src/SirtAlgorithm.cpp b/src/SirtAlgorithm.cpp index d9f3a65..ff25648 100644 --- a/src/SirtAlgorithm.cpp +++ b/src/SirtAlgorithm.cpp @@ -76,6 +76,7 @@ void CSirtAlgorithm::_clear() m_pDiffSinogram = NULL; m_pTmpVolume = NULL; + m_fLambda = 1.0f; m_iIterationCount = 0; } @@ -91,6 +92,7 @@ void CSirtAlgorithm::clear() ASTRA_DELETE(m_pDiffSinogram); ASTRA_DELETE(m_pTmpVolume); + m_fLambda = 1.0f; m_iIterationCount = 0; } @@ -128,6 +130,9 @@ bool CSirtAlgorithm::initialize(const Config& _cfg) return false; } + m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f); + CC.markOptionParsed("Relaxation"); + // init data objects and data projectors _init(); @@ -152,6 +157,8 @@ bool CSirtAlgorithm::initialize(CProjector2D* _pProjector, m_pSinogram = _pSinogram; m_pReconstruction = _pReconstruction; + m_fLambda = 1.0f; + // init data objects and data projectors _init(); @@ -248,7 +255,7 @@ void CSirtAlgorithm::run(int _iNrIterations) x = 1.0f / x; else x = 0.0f; - pfT[i] = x; + pfT[i] = m_fLambda * x; } pfT = m_pTotalRayLength->getData(); for (int i = 0; i < m_pTotalRayLength->getSize(); ++i) { @@ -296,7 +303,7 @@ void CSirtAlgorithm::run(int _iNrIterations) m_pTmpVolume->setData(0.0f); pBackProjector->project(); - // divide by pixel weights + // multiply with relaxation factor divided by pixel weights (*m_pTmpVolume) *= (*m_pTotalPixelWeight); (*m_pReconstruction) += (*m_pTmpVolume); -- cgit v1.2.3 From 16430239d04ff738a21146c410918c285552543f Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 23 Mar 2016 15:50:24 +0100 Subject: Add relaxation parameters to SIRT3D --- src/CudaSirtAlgorithm3D.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) (limited to 'src') diff --git a/src/CudaSirtAlgorithm3D.cpp b/src/CudaSirtAlgorithm3D.cpp index 605c470..c819f8e 100644 --- a/src/CudaSirtAlgorithm3D.cpp +++ b/src/CudaSirtAlgorithm3D.cpp @@ -56,6 +56,7 @@ CCudaSirtAlgorithm3D::CCudaSirtAlgorithm3D() m_iGPUIndex = -1; m_iVoxelSuperSampling = 1; m_iDetectorSuperSampling = 1; + m_fLambda = 1.0f; } //---------------------------------------------------------------------------------------- @@ -128,6 +129,8 @@ bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg) return false; } + m_fLambda = _cfg.self.getOptionNumerical("Relaxation"); + initializeFromProjector(); // Deprecated options @@ -135,6 +138,7 @@ bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg) m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling); m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex); m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); + CC.markOptionParsed("VoxelSuperSampling"); CC.markOptionParsed("DetectorSuperSampling"); CC.markOptionParsed("GPUIndex"); @@ -164,6 +168,8 @@ bool CCudaSirtAlgorithm3D::initialize(CProjector3D* _pProjector, clear(); } + m_fLambda = 1.0f; + // required classes m_pProjector = _pProjector; m_pSinogram = _pSinogram; @@ -224,6 +230,8 @@ void CCudaSirtAlgorithm3D::run(int _iNrIterations) ASTRA_ASSERT(ok); + m_pSirt->setRelaxation(m_fLambda); + m_bAstraSIRTInit = true; } -- cgit v1.2.3 From f9cc36d3507f7cde4d20165836d65a584ced720f Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 23 Mar 2016 18:21:42 +0100 Subject: Fix accumulating multiple raylengths in SART Thanks to @mohamedadaly for noticing. --- src/SartAlgorithm.cpp | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/SartAlgorithm.cpp b/src/SartAlgorithm.cpp index 9346160..f80df61 100644 --- a/src/SartAlgorithm.cpp +++ b/src/SartAlgorithm.cpp @@ -272,9 +272,8 @@ void CSartAlgorithm::run(int _iNrIterations) m_bShouldAbort = false; - int iIteration = 0; - // data projectors + CDataProjectorInterface* pFirstForwardProjector; CDataProjectorInterface* pForwardProjector; CDataProjectorInterface* pBackProjector; @@ -292,7 +291,7 @@ void CSartAlgorithm::run(int _iNrIterations) // first time forward projection data projector, // also computes total pixel weight and total ray length - pForwardProjector = dispatchDataProjector( + pFirstForwardProjector = dispatchDataProjector( m_pProjector, SinogramMaskPolicy(m_pSinogramMask), // sinogram mask ReconstructionMaskPolicy(m_pReconstructionMask), // reconstruction mask @@ -303,16 +302,30 @@ void CSartAlgorithm::run(int _iNrIterations) m_bUseSinogramMask, m_bUseReconstructionMask, true // options on/off ); + // forward projection data projector + pForwardProjector = dispatchDataProjector( + m_pProjector, + SinogramMaskPolicy(m_pSinogramMask), // sinogram mask + ReconstructionMaskPolicy(m_pReconstructionMask), // reconstruction mask + CombinePolicy( // 2 basic operations + DiffFPPolicy(m_pReconstruction, m_pDiffSinogram, m_pSinogram), // forward projection with difference calculation + TotalPixelWeightPolicy(m_pTotalPixelWeight)), // calculate the total pixel weights + m_bUseSinogramMask, m_bUseReconstructionMask, true // options on/off + ); + // iteration loop - for (; iIteration < _iNrIterations && !m_bShouldAbort; ++iIteration) { + for (int iIteration = 0; iIteration < _iNrIterations && !m_bShouldAbort; ++iIteration) { int iProjection = m_piProjectionOrder[m_iIterationCount % m_iProjectionCount]; // forward projection and difference calculation m_pTotalPixelWeight->setData(0.0f); - pForwardProjector->projectSingleProjection(iProjection); + if (iIteration < m_iProjectionCount) + pFirstForwardProjector->projectSingleProjection(iProjection); + else + pForwardProjector->projectSingleProjection(iProjection); // backprojection pBackProjector->projectSingleProjection(iProjection); // update iteration count @@ -325,6 +338,7 @@ void CSartAlgorithm::run(int _iNrIterations) } + ASTRA_DELETE(pFirstForwardProjector); ASTRA_DELETE(pForwardProjector); ASTRA_DELETE(pBackProjector); -- cgit v1.2.3 From 5c6a9523523c1680e898b4cf897531cd75435f22 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 24 Mar 2016 17:34:34 +0100 Subject: Fix cppcheck warnings --- src/CudaDataOperationAlgorithm.cpp | 2 +- src/Float32VolumeData3DMemory.cpp | 1 - src/XMLNode.cpp | 4 ++-- 3 files changed, 3 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/CudaDataOperationAlgorithm.cpp b/src/CudaDataOperationAlgorithm.cpp index 15886a4..82b676b 100644 --- a/src/CudaDataOperationAlgorithm.cpp +++ b/src/CudaDataOperationAlgorithm.cpp @@ -76,7 +76,7 @@ bool CCudaDataOperationAlgorithm::initialize(const Config& _cfg) node = _cfg.self.getSingleNode("DataId"); ASTRA_CONFIG_CHECK(node, "CCudaDataOperationAlgorithm", "No DataId tag specified."); vector data = node.getContentArray(); - for (vector::iterator it = data.begin(); it != data.end(); it++){ + for (vector::iterator it = data.begin(); it != data.end(); ++it){ int id = StringUtil::stringToInt(*it); m_pData.push_back(dynamic_cast(CData2DManager::getSingleton().get(id))); } diff --git a/src/Float32VolumeData3DMemory.cpp b/src/Float32VolumeData3DMemory.cpp index af45cb9..14adb1a 100644 --- a/src/Float32VolumeData3DMemory.cpp +++ b/src/Float32VolumeData3DMemory.cpp @@ -136,7 +136,6 @@ CFloat32VolumeData2D * CFloat32VolumeData3DMemory::fetchSliceZ(int _iSliceIndex) CFloat32VolumeData2D* res = new CFloat32VolumeData2D(&volGeom); // copy data - int iSliceCount = m_pGeometry->getGridSliceCount(); float * pfTargetData = res->getData(); for(int iRowIndex = 0; iRowIndex < iRowCount; iRowIndex++) { diff --git a/src/XMLNode.cpp b/src/XMLNode.cpp index 40a9b22..cf268c2 100644 --- a/src/XMLNode.cpp +++ b/src/XMLNode.cpp @@ -158,7 +158,7 @@ vector XMLNode::getContentArray() const vector res(iSize); // loop all list item nodes list nodes = getNodes("ListItem"); - for (list::iterator it = nodes.begin(); it != nodes.end(); it++) { + for (list::iterator it = nodes.begin(); it != nodes.end(); ++it) { int iIndex = it->getAttributeNumerical("index"); string sValue = it->getAttribute("value"); ASTRA_ASSERT(iIndex < iSize); @@ -290,7 +290,7 @@ vector XMLNode::getOptionNumericalArray(string _sKey) const if (!hasOption(_sKey)) return vector(); list nodes = getNodes("Option"); - for (list::iterator it = nodes.begin(); it != nodes.end(); it++) { + for (list::iterator it = nodes.begin(); it != nodes.end(); ++it) { if (it->getAttribute("key") == _sKey) { vector vals = it->getContentNumericalArray(); return vals; -- cgit v1.2.3 From ab583834eec75d23e16b8f205ab65788d8df6ffe Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 24 Mar 2016 17:34:50 +0100 Subject: Fix memory leak --- src/FilteredBackProjectionAlgorithm.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) (limited to 'src') diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp index c195578..ccbfec6 100644 --- a/src/FilteredBackProjectionAlgorithm.cpp +++ b/src/FilteredBackProjectionAlgorithm.cpp @@ -117,12 +117,10 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) int angleCount = projectionIndex.size(); int detectorCount = m_pProjector->getProjectionGeometry()->getDetectorCount(); + // TODO: There is no need to allocate this. Better just + // create the CFloat32ProjectionData2D object directly, and use its + // memory. float32 * sinogramData2D = new float32[angleCount* detectorCount]; - float32 ** sinogramData = new float32*[angleCount]; - for (int i = 0; i < angleCount; i++) - { - sinogramData[i] = &(sinogramData2D[i * detectorCount]); - } float32 * projectionAngles = new float32[angleCount]; float32 detectorWidth = m_pProjector->getProjectionGeometry()->getDetectorWidth(); @@ -130,6 +128,8 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) for (int i = 0; i < angleCount; i ++) { if (projectionIndex[i] > m_pProjector->getProjectionGeometry()->getProjectionAngleCount() -1 ) { + delete[] sinogramData2D; + delete[] projectionAngles; ASTRA_ERROR("Invalid Projection Index"); return false; } else { @@ -139,7 +139,6 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) { sinogramData2D[i*detectorCount+ iDetector] = m_pSinogram->getData2D()[orgIndex][iDetector]; } -// sinogramData[i] = m_pSinogram->getSingleProjectionData(projectionIndex[i]); projectionAngles[i] = m_pProjector->getProjectionGeometry()->getProjectionAngle((int)projectionIndex[i] ); } @@ -148,6 +147,9 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) CParallelProjectionGeometry2D * pg = new CParallelProjectionGeometry2D(angleCount, detectorCount,detectorWidth,projectionAngles); m_pProjector = new CParallelBeamLineKernelProjector2D(pg,m_pReconstruction->getGeometry()); m_pSinogram = new CFloat32ProjectionData2D(pg, sinogramData2D); + + delete[] sinogramData2D; + delete[] projectionAngles; } // TODO: check that the angles are linearly spaced between 0 and pi -- cgit v1.2.3 From 547def0ea6e3eab07b7e4c48cee6d6a81f6155e1 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 18 Apr 2016 11:43:48 +0200 Subject: Fix stringToInt parsing doubles --- src/Utilities.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/Utilities.cpp b/src/Utilities.cpp index 4b80503..c9740bf 100644 --- a/src/Utilities.cpp +++ b/src/Utilities.cpp @@ -42,7 +42,7 @@ namespace StringUtil { int stringToInt(const std::string& s) { - double i; + int i; std::istringstream iss(s); iss.imbue(std::locale::classic()); iss >> i; -- cgit v1.2.3