diff options
73 files changed, 1693 insertions, 874 deletions
diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..4a179f1 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,45 @@ +language: python + +python: + - "2.7" + - "3.5" + +os: + - linux + +sudo: false + +addons: + apt: + packages: + - libboost-all-dev + - nvidia-common + - nvidia-current + - nvidia-cuda-toolkit + - nvidia-cuda-dev +env: + - CUDA=yes + - CUDA=no + +before_install: + - if [[ "$TRAVIS_PYTHON_VERSION" == "2.7" ]]; then + wget https://repo.continuum.io/miniconda/Miniconda-latest-Linux-x86_64.sh -O miniconda.sh; + else + wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh; + fi + - bash miniconda.sh -b -p $HOME/miniconda + - export PATH="$HOME/miniconda/bin:$PATH" + - conda config --set always_yes yes --set changeps1 no + - conda update conda + +install: + - conda install python=$TRAVIS_PYTHON_VERSION six numpy scipy cython + - conda info -a + - cd build/linux + - ./autogen.sh + - if [ $CUDA == yes ]; then ./configure --prefix=$HOME/astra --with-python --with-cuda; else ./configure --prefix=$HOME/astra --with-python --without-cuda; fi + - make -j 4 + - make install + +script: + - LD_LIBRARY_PATH=$HOME/astra/lib/:$LD_LIBRARY_PATH PYTHONPATH=$HOME/astra/python/:$PYTHONPATH python -c "import astra" diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in index 3018674..a199bf6 100644 --- a/build/linux/Makefile.in +++ b/build/linux/Makefile.in @@ -62,8 +62,6 @@ PYCPPFLAGS := $(CPPFLAGS) PYCPPFLAGS += -I../include PYLDFLAGS = $(LDFLAGS) PYLDFLAGS += -L$(abs_top_builddir)/.libs -LIBS += -l$(PYLIBVER) -LDFLAGS += -L$(PYLIBDIR) endif # This is below where PYCPPFLAGS copies CPPFLAGS. The python code is built @@ -97,6 +95,11 @@ ifeq ($(cuda),yes) MEXFLAGS += -DASTRA_CUDA endif +ifeq ($(python),yes) +MEXPYLDFLAGS='$$LDFLAGS $(LDFLAGS) -L$(PYLIBDIR)' +MEXPYLIBS=$(MEXLIBS) -l$(PYLIBVER) +endif + endif LIBDIR=/usr/local/lib @@ -147,6 +150,7 @@ BASE_OBJECTS=\ src/ParallelProjectionGeometry3D.lo \ src/ParallelVecProjectionGeometry3D.lo \ src/PlatformDepSystemCode.lo \ + src/PluginAlgorithm.lo \ src/ProjectionGeometry2D.lo \ src/ProjectionGeometry3D.lo \ src/Projector2D.lo \ @@ -252,7 +256,6 @@ MATLAB_MEX=\ matlab/mex/astra_mex_direct_c.$(MEXSUFFIX) ifeq ($(python),yes) -ALL_OBJECTS+=src/PluginAlgorithm.lo MATLAB_MEX+=matlab/mex/astra_mex_plugin_c.$(MEXSUFFIX) endif @@ -267,6 +270,11 @@ mex: $(MATLAB_MEX) %.$(MEXSUFFIX): %.o $(MATLAB_CXX_OBJECTS) libastra.la $(MEX) LDFLAGS=$(MEXLDFLAGS) $(MEXFLAGS) $(LIBS) $(MEXLIBS) -lastra -output $* $*.o $(MATLAB_CXX_OBJECTS) + +ifeq ($(python),yes) +matlab/mex/astra_mex_plugin_c.$(MEXSUFFIX): matlab/mex/astra_mex_plugin_c.o $(MATLAB_CXX_OBJECTS) libastra.la + $(MEX) LDFLAGS=$(MEXPYLDFLAGS) $(MEXFLAGS) $(LIBS) $(MEXPYLIBS) -lastra -output matlab/mex/astra_mex_plugin_c $< $(MATLAB_CXX_OBJECTS) +endif endif ifeq ($(python),yes) @@ -307,8 +315,10 @@ ifeq ($(cuda),yes) ifeq ($(gen_static_libs),yes) @$(NVCC) $(NVCCFLAGS) -c $(<) -o $*.o >/dev/null 2>&1 endif - @# Generate a .d file, with target name $*.lo - @$(NVCC) $(NVCCFLAGS) -M $(<) -MT $(*F).lo -odir $(*D) -o $(*D)/$(DEPDIR)/$(*F).d + @# Generate a .d file, and change the target name in it from .o to .lo + @$(NVCC) $(NVCCFLAGS) -M $(<) -odir $(*D) -o $(*D)/$(DEPDIR)/$(*F).d2 + @sed '1s/\.o :/.lo :/' < $(*D)/$(DEPDIR)/$(*F).d2 > $(*D)/$(DEPDIR)/$(*F).d + @rm -f $(*D)/$(DEPDIR)/$(*F).d2 @# Generate empty targets for all dependencies listed in the .d file. @# This mimics gcc's -MP option. @for x in `cat $(*D)/$(DEPDIR)/$(*F).d`; do if test a$$x != a: -a a$$x != a\\; then echo -e "\n$$x:\n" >> $(*D)/$(DEPDIR)/$(*F).d; fi; done @@ -411,7 +421,7 @@ $(srcdir)/configure: $(srcdir)/configure.ac @echo "configure.ac has been changed. Regenerating configure script" cd $(srcdir) && $(SHELL) ./autogen.sh -.PHONY: all mex test clean distclean install install-libraries +.PHONY: all mex test clean distclean install install-libraries py python-root-install install-python # don't remove intermediate files: .SECONDARY: diff --git a/build/msvc/gen.py b/build/msvc/gen.py index 999710f..cc69a62 100644 --- a/build/msvc/gen.py +++ b/build/msvc/gen.py @@ -1,6 +1,8 @@ from __future__ import print_function import sys import os +import codecs +import six vcppguid = "8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942" # C++ project siguid = "2150E333-8FDC-42A3-9474-1A3956D46DE8" # project group @@ -428,7 +430,10 @@ P_astra["files"].sort() projects = [ P_astra, F_astra_mex, P0, P1, P2, P3, P4, P5, P6, P7, P8 ] -bom = "\xef\xbb\xbf" +if six.PY2: + bom = "\xef\xbb\xbf" +else: + bom = codecs.BOM_UTF8.decode("utf-8") class Configuration: def __init__(self, debug, cuda, x64): diff --git a/cuda/2d/sart.cu b/cuda/2d/sart.cu index e5cb5bb..c8608a3 100644 --- a/cuda/2d/sart.cu +++ b/cuda/2d/sart.cu @@ -71,6 +71,8 @@ SART::SART() : ReconAlgo() projectionCount = 0; iteration = 0; customOrder = false; + + fRelaxation = 1.0f; } @@ -98,6 +100,7 @@ void SART::reset() projectionCount = 0; iteration = 0; customOrder = false; + fRelaxation = 1.0f; ReconAlgo::reset(); } @@ -200,10 +203,10 @@ bool SART::iterate(unsigned int iterations) // BP, mask, and add back // TODO: Try putting the masking directly in the BP zeroVolumeData(D_tmpData, tmpPitch, dims); - callBP_SART(D_tmpData, tmpPitch, D_projData, projPitch, angle, 1.0f); + callBP_SART(D_tmpData, tmpPitch, D_projData, projPitch, angle, fRelaxation); processVol<opAddMul>(D_volumeData, D_maskData, D_tmpData, volumePitch, dims); } else { - callBP_SART(D_volumeData, volumePitch, D_projData, projPitch, angle, 1.0f); + callBP_SART(D_volumeData, volumePitch, D_projData, projPitch, angle, fRelaxation); } if (useMinConstraint) diff --git a/cuda/2d/sart.h b/cuda/2d/sart.h index 7dcd641..eff9ecf 100644 --- a/cuda/2d/sart.h +++ b/cuda/2d/sart.h @@ -50,6 +50,8 @@ public: virtual float computeDiffNorm(); + void setRelaxation(float r) { fRelaxation = r; } + protected: void reset(); bool precomputeWeights(); @@ -78,6 +80,8 @@ protected: // Geometry-specific precomputed data float* D_lineWeight; unsigned int linePitch; + + float fRelaxation; }; } diff --git a/cuda/2d/sirt.cu b/cuda/2d/sirt.cu index 162ee77..4baaccb 100644 --- a/cuda/2d/sirt.cu +++ b/cuda/2d/sirt.cu @@ -50,6 +50,8 @@ SIRT::SIRT() : ReconAlgo() D_minMaskData = 0; D_maxMaskData = 0; + fRelaxation = 1.0f; + freeMinMaxMasks = false; } @@ -83,6 +85,8 @@ void SIRT::reset() useVolumeMask = false; useSinogramMask = false; + fRelaxation = 1.0f; + ReconAlgo::reset(); } @@ -139,6 +143,9 @@ bool SIRT::precomputeWeights() processVol<opMul>(D_pixelWeight, D_maskData, pixelPitch, dims); } + // Also fold the relaxation factor into pixel weights + processVol<opMul>(D_pixelWeight, fRelaxation, pixelPitch, dims); + return true; } @@ -253,6 +260,7 @@ bool SIRT::iterate(unsigned int iterations) callBP(D_tmpData, tmpPitch, D_projData, projPitch, 1.0f); + // pixel weights also contain the volume mask and relaxation factor processVol<opAddMul>(D_volumeData, D_pixelWeight, D_tmpData, volumePitch, dims); if (useMinConstraint) diff --git a/cuda/2d/sirt.h b/cuda/2d/sirt.h index 21094a1..bc913f4 100644 --- a/cuda/2d/sirt.h +++ b/cuda/2d/sirt.h @@ -53,6 +53,8 @@ public: bool uploadMinMaxMasks(const float* minMaskData, const float* maxMaskData, unsigned int pitch); + void setRelaxation(float r) { fRelaxation = r; } + virtual bool iterate(unsigned int iterations); virtual float computeDiffNorm(); @@ -81,6 +83,8 @@ protected: unsigned int minMaskPitch; float* D_maxMaskData; unsigned int maxMaskPitch; + + float fRelaxation; }; bool doSIRT(float* D_volumeData, unsigned int volumePitch, diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu index 42cece2..dd6d551 100644 --- a/cuda/3d/astra3d.cu +++ b/cuda/3d/astra3d.cu @@ -263,6 +263,7 @@ public: float* angles; float fOriginSourceDistance; float fOriginDetectorDistance; + float fRelaxation; SConeProjection* projs; SPar3DProjection* parprojs; @@ -302,6 +303,8 @@ AstraSIRT3d::AstraSIRT3d() pData->projs = 0; pData->parprojs = 0; + pData->fRelaxation = 1.0f; + pData->initialized = false; pData->setStartReconstruction = false; @@ -380,6 +383,14 @@ bool AstraSIRT3d::enableSuperSampling(unsigned int iVoxelSuperSampling, return true; } +void AstraSIRT3d::setRelaxation(float r) +{ + if (pData->initialized) + return; + + pData->fRelaxation = r; +} + bool AstraSIRT3d::enableVolumeMask() { if (pData->initialized) @@ -439,6 +450,8 @@ bool AstraSIRT3d::init() if (!ok) return false; + pData->sirt.setRelaxation(pData->fRelaxation); + pData->D_volumeData = allocateVolumeData(pData->dims); ok = pData->D_volumeData.ptr; if (!ok) diff --git a/cuda/3d/astra3d.h b/cuda/3d/astra3d.h index 1f1cee2..3345ab8 100644 --- a/cuda/3d/astra3d.h +++ b/cuda/3d/astra3d.h @@ -67,6 +67,8 @@ public: bool enableSuperSampling(unsigned int iVoxelSuperSampling, unsigned int iDetectorSuperSampling); + void setRelaxation(float r); + // Enable volume/sinogram masks // // This may optionally be called before init(). diff --git a/cuda/3d/mem3d.cu b/cuda/3d/mem3d.cu index d09c203..9a239da 100644 --- a/cuda/3d/mem3d.cu +++ b/cuda/3d/mem3d.cu @@ -62,6 +62,25 @@ size_t availableGPUMemory() return free; } +int maxBlockDimension() +{ + int dev; + cudaError_t err = cudaGetDevice(&dev); + if (err != cudaSuccess) { + ASTRA_WARN("Error querying device"); + return 0; + } + + cudaDeviceProp props; + err = cudaGetDeviceProperties(&props, dev); + if (err != cudaSuccess) { + ASTRA_WARN("Error querying device %d properties", dev); + return 0; + } + + return std::min(props.maxTexture3D[0], std::min(props.maxTexture3D[1], props.maxTexture3D[2])); +} + MemHandle3D allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, Mem3DZeroMode zero) { SMemHandle3D_internal hnd; diff --git a/cuda/3d/mem3d.h b/cuda/3d/mem3d.h index acb72cb..6fff80b 100644 --- a/cuda/3d/mem3d.h +++ b/cuda/3d/mem3d.h @@ -78,6 +78,7 @@ enum Mem3DZeroMode { }; size_t availableGPUMemory(); +int maxBlockDimension(); MemHandle3D allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, Mem3DZeroMode zero); diff --git a/cuda/3d/sirt3d.cu b/cuda/3d/sirt3d.cu index 4cb35c3..ba6a159 100644 --- a/cuda/3d/sirt3d.cu +++ b/cuda/3d/sirt3d.cu @@ -59,6 +59,8 @@ SIRT::SIRT() : ReconAlgo3D() useMinConstraint = false; useMaxConstraint = false; + + fRelaxation = 1.0f; } @@ -89,6 +91,8 @@ void SIRT::reset() useVolumeMask = false; useSinogramMask = false; + fRelaxation = 1.0f; + ReconAlgo3D::reset(); } @@ -196,6 +200,8 @@ bool SIRT::precomputeWeights() // scale pixel weights with mask to zero out masked pixels processVol3D<opMul>(D_pixelWeight, D_maskData, dims); } + processVol3D<opMul>(D_pixelWeight, fRelaxation, dims); + return true; } @@ -307,7 +313,7 @@ bool SIRT::iterate(unsigned int iterations) } #endif - + // pixel weights also contain the volume mask and relaxation factor processVol3D<opAddMul>(D_volumeData, D_tmpData, D_pixelWeight, dims); if (useMinConstraint) diff --git a/cuda/3d/sirt3d.h b/cuda/3d/sirt3d.h index bb3864a..5e93deb 100644 --- a/cuda/3d/sirt3d.h +++ b/cuda/3d/sirt3d.h @@ -48,6 +48,9 @@ public: // init should be called after setting all geometry bool init(); + // Set relaxation factor. This may be called after init and before iterate. + void setRelaxation(float r) { fRelaxation = r; } + // setVolumeMask should be called after init and before iterate, // but only if enableVolumeMask was called before init. // It may be called again after iterate. @@ -91,6 +94,8 @@ protected: float fMinConstraint; float fMaxConstraint; + float fRelaxation; + cudaPitchedPtr D_maskData; cudaPitchedPtr D_smaskData; diff --git a/include/astra/Algorithm.h b/include/astra/Algorithm.h index 049417c..18c465f 100644 --- a/include/astra/Algorithm.h +++ b/include/astra/Algorithm.h @@ -81,7 +81,7 @@ public: * * @return initialized */ - bool isInitialized(); + bool isInitialized() const; /** get a description of the class * @@ -133,7 +133,7 @@ private: // inline functions inline std::string CAlgorithm::description() const { return "Algorithm"; }; -inline bool CAlgorithm::isInitialized() { return m_bIsInitialized; } +inline bool CAlgorithm::isInitialized() const { return m_bIsInitialized; } } // end namespace diff --git a/include/astra/ArtAlgorithm.h b/include/astra/ArtAlgorithm.h index d232b95..1ad9f3f 100644 --- a/include/astra/ArtAlgorithm.h +++ b/include/astra/ArtAlgorithm.h @@ -59,7 +59,7 @@ namespace astra { * \astra_xml_item_option{MinConstraintValue, float, 0, Minimum constraint value.} * \astra_xml_item_option{UseMaxConstraint, bool, false, Use maximum value constraint.} * \astra_xml_item_option{MaxConstraintValue, float, 255, Maximum constraint value.} - * \astra_xml_item_option{Lamda, float, 1, The relaxation factor.} + * \astra_xml_item_option{Relaxation, float, 1, The relaxation factor.} * \astra_xml_item_option{RayOrder, string, "sequential", the order in which the rays are updated. 'sequential' or 'custom'} * \astra_xml_item_option{RayOrderList, n by 2 vector of float, not used, if RayOrder='custom': use this ray order. Each row consist of a projection id and detector id.} * @@ -73,7 +73,7 @@ namespace astra { * cfg.option.UseMinConstraint = 'yes';\n * cfg.option.UseMaxConstraint = 'yes';\n * cfg.option.MaxConstraintValue = 1024;\n - * cfg.option.Lamda = 0.7;\n + * cfg.option.Relaxation = 0.7;\n * cfg.option.RayOrder = 'custom';\n * cfg.option.RayOrderList = [0\,0; 0\,2; 1\,0];\n * alg_id = astra_mex_algorithm('create'\, cfg);\n diff --git a/include/astra/AstraObjectFactory.h b/include/astra/AstraObjectFactory.h index 325989e..6af9cd8 100644 --- a/include/astra/AstraObjectFactory.h +++ b/include/astra/AstraObjectFactory.h @@ -155,8 +155,11 @@ class _AstraExport CAlgorithmFactory : public CAstraObjectFactory<CAlgorithm, Al template <> inline CAlgorithm* CAstraObjectFactory<CAlgorithm, AlgorithmTypeList>::findPlugin(std::string _sType) { - CPluginAlgorithmFactory *fac = CPluginAlgorithmFactory::getSingletonPtr(); - return fac->getPlugin(_sType); + CPluginAlgorithmFactory *fac = CPluginAlgorithmFactory::getFactory(); + if (fac) + return fac->getPlugin(_sType); + else + return 0; } #endif diff --git a/include/astra/AstraObjectManager.h b/include/astra/AstraObjectManager.h index 895f955..ad89c2a 100644 --- a/include/astra/AstraObjectManager.h +++ b/include/astra/AstraObjectManager.h @@ -52,17 +52,49 @@ namespace astra { * among all ObjectManagers. */ +class CAstraObjectManagerBase { +public: + virtual std::string getInfo(int index) const =0; + virtual void remove(int index) =0; + virtual std::string getType() const =0; +}; -class CAstraIndexManager { -protected: - /** The index of the previously stored data object. + +class CAstraIndexManager : public Singleton<CAstraIndexManager> { +public: + CAstraIndexManager() : m_iLastIndex(0) { } + + int store(CAstraObjectManagerBase* m) { + m_table[++m_iLastIndex] = m; + return m_iLastIndex; + } + + CAstraObjectManagerBase* get(int index) const { + std::map<int, CAstraObjectManagerBase*>::const_iterator i; + i = m_table.find(index); + if (i != m_table.end()) + return i->second; + else + return 0; + } + + void remove(int index) { + std::map<int, CAstraObjectManagerBase*>::iterator i; + i = m_table.find(index); + if (i != m_table.end()) + m_table.erase(i); + } + +private: + /** The index last handed out */ - static int m_iPreviousIndex; + int m_iLastIndex; + std::map<int, CAstraObjectManagerBase*> m_table; }; template <typename T> -class CAstraObjectManager : public Singleton<CAstraObjectManager<T> >, CAstraIndexManager { +class CAstraObjectManager : public CAstraObjectManagerBase { public: @@ -117,7 +149,11 @@ public: */ void clear(); - /** Get info. + /** Get info of object. + */ + std::string getInfo(int index) const; + + /** Get list with info of all managed objects. */ std::string info(); @@ -149,9 +185,9 @@ CAstraObjectManager<T>::~CAstraObjectManager() template <typename T> int CAstraObjectManager<T>::store(T* _pDataObject) { - m_iPreviousIndex++; - m_mIndexToObject[m_iPreviousIndex] = _pDataObject; - return m_iPreviousIndex; + int iIndex = CAstraIndexManager::getSingleton().store(this); + m_mIndexToObject[iIndex] = _pDataObject; + return iIndex; } //---------------------------------------------------------------------------------------- @@ -180,15 +216,16 @@ T* CAstraObjectManager<T>::get(int _iIndex) const template <typename T> void CAstraObjectManager<T>::remove(int _iIndex) { - if (!hasIndex(_iIndex)) { - return; - } // find data typename map<int,T*>::iterator it = m_mIndexToObject.find(_iIndex); + if (it == m_mIndexToObject.end()) + return; // delete data delete (*it).second; // delete from map - m_mIndexToObject.erase(it); + m_mIndexToObject.erase(it); + + CAstraIndexManager::getSingleton().remove(_iIndex); } //---------------------------------------------------------------------------------------- @@ -220,19 +257,29 @@ void CAstraObjectManager<T>::clear() //---------------------------------------------------------------------------------------- // Print info to string template <typename T> +std::string CAstraObjectManager<T>::getInfo(int index) const { + typename map<int,T*>::const_iterator it = m_mIndexToObject.find(index); + if (it == m_mIndexToObject.end()) + return ""; + const T* pObject = it->second; + std::stringstream res; + res << index << " \t"; + if (pObject->isInitialized()) { + res << "v "; + } else { + res << "x "; + } + res << pObject->description(); + return res.str(); +} + +template <typename T> std::string CAstraObjectManager<T>::info() { std::stringstream res; res << "id init description" << std::endl; res << "-----------------------------------------" << std::endl; - for (typename map<int,T*>::iterator it = m_mIndexToObject.begin(); it != m_mIndexToObject.end(); it++) { - res << (*it).first << " \t"; - T* pObject = m_mIndexToObject[(*it).first]; - if (pObject->isInitialized()) { - res << "v "; - } else { - res << "x "; - } - res << pObject->description() << endl; + for (typename map<int,T*>::const_iterator it = m_mIndexToObject.begin(); it != m_mIndexToObject.end(); it++) { + res << getInfo(it->first) << endl; } res << "-----------------------------------------" << std::endl; return res.str(); @@ -247,42 +294,60 @@ std::string CAstraObjectManager<T>::info() { * assigned to each data object by which it can be accessed in the future. * Indices are always >= 1. */ -class _AstraExport CProjector2DManager : public CAstraObjectManager<CProjector2D>{}; +class _AstraExport CProjector2DManager : public Singleton<CProjector2DManager>, public CAstraObjectManager<CProjector2D> +{ + virtual std::string getType() const { return "projector2d"; } +}; /** * This class contains functionality to store 3D projector objects. A unique index handle will be * assigned to each data object by which it can be accessed in the future. * Indices are always >= 1. */ -class _AstraExport CProjector3DManager : public CAstraObjectManager<CProjector3D>{}; +class _AstraExport CProjector3DManager : public Singleton<CProjector3DManager>, public CAstraObjectManager<CProjector3D> +{ + virtual std::string getType() const { return "projector3d"; } +}; /** * This class contains functionality to store 2D data objects. A unique index handle will be * assigned to each data object by which it can be accessed in the future. * Indices are always >= 1. */ -class _AstraExport CData2DManager : public CAstraObjectManager<CFloat32Data2D>{}; +class _AstraExport CData2DManager : public Singleton<CData2DManager>, public CAstraObjectManager<CFloat32Data2D> +{ + virtual std::string getType() const { return "data2d"; } +}; /** * This class contains functionality to store 3D data objects. A unique index handle will be * assigned to each data object by which it can be accessed in the future. * Indices are always >= 1. */ -class _AstraExport CData3DManager : public CAstraObjectManager<CFloat32Data3D>{}; +class _AstraExport CData3DManager : public Singleton<CData3DManager>, public CAstraObjectManager<CFloat32Data3D> +{ + virtual std::string getType() const { return "data3d"; } +}; /** * This class contains functionality to store algorithm objects. A unique index handle will be * assigned to each data object by which it can be accessed in the future. * Indices are always >= 1. */ -class _AstraExport CAlgorithmManager : public CAstraObjectManager<CAlgorithm>{}; +class _AstraExport CAlgorithmManager : public Singleton<CAlgorithmManager>, public CAstraObjectManager<CAlgorithm> +{ + virtual std::string getType() const { return "algorithm"; } +}; /** * This class contains functionality to store matrix objects. A unique index handle will be * assigned to each data object by which it can be accessed in the future. * Indices are always >= 1. */ -class _AstraExport CMatrixManager : public CAstraObjectManager<CSparseMatrix>{}; +class _AstraExport CMatrixManager : public Singleton<CMatrixManager>, public CAstraObjectManager<CSparseMatrix> +{ + virtual std::string getType() const { return "matrix"; } +}; } // end namespace diff --git a/include/astra/CompositeGeometryManager.h b/include/astra/CompositeGeometryManager.h index 4338994..18dd72f 100644 --- a/include/astra/CompositeGeometryManager.h +++ b/include/astra/CompositeGeometryManager.h @@ -79,7 +79,9 @@ public: bool uploadToGPU(); bool downloadFromGPU(/*mode?*/); - virtual TPartList split(size_t maxSize, int div) = 0; + virtual void splitX(TPartList& out, size_t maxSize, size_t maxDim, int div) = 0; + virtual void splitY(TPartList& out, size_t maxSize, size_t maxDim, int div) = 0; + virtual void splitZ(TPartList& out, size_t maxSize, size_t maxDim, int div) = 0; virtual CPart* reduce(const CPart *other) = 0; virtual void getDims(size_t &x, size_t &y, size_t &z) = 0; size_t getSize(); @@ -93,7 +95,9 @@ public: CVolumeGeometry3D* pGeom; - virtual TPartList split(size_t maxSize, int div); + virtual void splitX(TPartList& out, size_t maxSize, size_t maxDim, int div); + virtual void splitY(TPartList& out, size_t maxSize, size_t maxDim, int div); + virtual void splitZ(TPartList& out, size_t maxSize, size_t maxDim, int div); virtual CPart* reduce(const CPart *other); virtual void getDims(size_t &x, size_t &y, size_t &z); @@ -107,7 +111,9 @@ public: CProjectionGeometry3D* pGeom; - virtual TPartList split(size_t maxSize, int div); + virtual void splitX(TPartList& out, size_t maxSize, size_t maxDim, int div); + virtual void splitY(TPartList& out, size_t maxSize, size_t maxDim, int div); + virtual void splitZ(TPartList& out, size_t maxSize, size_t maxDim, int div); virtual CPart* reduce(const CPart *other); virtual void getDims(size_t &x, size_t &y, size_t &z); diff --git a/include/astra/CudaReconstructionAlgorithm2D.h b/include/astra/CudaReconstructionAlgorithm2D.h index dc93a1a..bb5f2a7 100644 --- a/include/astra/CudaReconstructionAlgorithm2D.h +++ b/include/astra/CudaReconstructionAlgorithm2D.h @@ -141,6 +141,9 @@ protected: */ bool setupGeometry(); + /** Initialize CUDA algorithm. For internal use only. + */ + virtual void initCUDAAlgorithm(); /** The internally used CUDA algorithm object */ diff --git a/include/astra/CudaSartAlgorithm.h b/include/astra/CudaSartAlgorithm.h index c22dc4f..e5e18f0 100644 --- a/include/astra/CudaSartAlgorithm.h +++ b/include/astra/CudaSartAlgorithm.h @@ -46,6 +46,7 @@ namespace astra { * \astra_xml_item{ProjectionDataId, integer, Identifier of a projection data object as it is stored in the DataManager.} * \astra_xml_item{ReconstructionDataId, integer, Identifier of a volume data object as it is stored in the DataManager.} * \astra_xml_item_option{ReconstructionMaskId, integer, not used, Identifier of a volume data object that acts as a reconstruction mask. 0 = reconstruct on this pixel. 1 = don't reconstruct on this pixel.} + * \astra_xml_item_option{Relaxation, float, 1, The relaxation factor.} * * \par MATLAB example * \astra_code{ @@ -53,6 +54,7 @@ namespace astra { * cfg.ProjectionDataId = sino_id;\n * cfg.ReconstructionDataId = recon_id;\n * cfg.option.ReconstructionMaskId = mask_id;\n + * cfg.option.Relaxation = 1.0;\n * alg_id = astra_mex_algorithm('create'\, cfg);\n * astra_mex_algorithm('iterate'\, alg_id\, 10);\n * astra_mex_algorithm('delete'\, alg_id);\n @@ -97,6 +99,14 @@ public: * @return description string */ virtual std::string description() const; + +protected: + + /** Relaxation factor + */ + float m_fLambda; + + virtual void initCUDAAlgorithm(); }; // inline functions diff --git a/include/astra/CudaSirtAlgorithm.h b/include/astra/CudaSirtAlgorithm.h index 929ac30..3cd8acc 100644 --- a/include/astra/CudaSirtAlgorithm.h +++ b/include/astra/CudaSirtAlgorithm.h @@ -53,6 +53,7 @@ namespace astra { * \astra_xml_item{ProjectionDataId, integer, Identifier of a projection data object as it is stored in the DataManager.} * \astra_xml_item{ReconstructionDataId, integer, Identifier of a volume data object as it is stored in the DataManager.} * \astra_xml_item_option{ReconstructionMaskId, integer, not used, Identifier of a volume data object that acts as a reconstruction mask. 0 = reconstruct on this pixel. 1 = don't reconstruct on this pixel.} + * \astra_xml_item_option{Relaxation, float, 1, Relaxation parameter.} * * \par MATLAB example * \astra_code{ @@ -62,6 +63,7 @@ namespace astra { * cfg.ProjectionDataId = sino_id;\n * cfg.ReconstructionDataId = recon_id;\n * cfg.option.ReconstructionMaskId = mask_id;\n + * cfg.option.Relaxation = 1.0;\n * alg_id = astra_mex_algorithm('create'\, cfg);\n * astra_mex_algorithm('iterate'\, alg_id\, 10);\n * astra_mex_algorithm('delete'\, alg_id);\n @@ -93,8 +95,6 @@ public: */ virtual bool initialize(const Config& _cfg); - virtual void run(int _iNrIterations); - /** Initialize class. * * @param _pProjector Projector Object. (Optional) @@ -114,6 +114,12 @@ public: protected: CFloat32VolumeData2D* m_pMinMask; CFloat32VolumeData2D* m_pMaxMask; + + /** Relaxation factor + */ + float m_fLambda; + + virtual void initCUDAAlgorithm(); }; // inline functions diff --git a/include/astra/CudaSirtAlgorithm3D.h b/include/astra/CudaSirtAlgorithm3D.h index 379720e..60191cd 100644 --- a/include/astra/CudaSirtAlgorithm3D.h +++ b/include/astra/CudaSirtAlgorithm3D.h @@ -50,7 +50,7 @@ class AstraSIRT3d; * * The update step of pixel \f$v_j\f$ for iteration \f$k\f$ is given by: * \f[ - * v_j^{(k+1)} = v_j^{(k)} + \alpha \sum_{i=1}^{M} \left( \frac{w_{ij}\left( p_i - \sum_{r=1}^{N} w_{ir}v_r^{(k)}\right)}{\sum_{k=1}^{N} w_{ik}} \right) \frac{1}{\sum_{l=1}^{M}w_{lj}} + * v_j^{(k+1)} = v_j^{(k)} + \lambda \sum_{i=1}^{M} \left( \frac{w_{ij}\left( p_i - \sum_{r=1}^{N} w_{ir}v_r^{(k)}\right)}{\sum_{k=1}^{N} w_{ik}} \right) \frac{1}{\sum_{l=1}^{M}w_{lj}} * \f] * * \par XML Configuration @@ -175,6 +175,7 @@ protected: bool m_bAstraSIRTInit; int m_iDetectorSuperSampling; int m_iVoxelSuperSampling; + float m_fLambda; void initializeFromProjector(); }; diff --git a/include/astra/DataProjectorPolicies.h b/include/astra/DataProjectorPolicies.h index c258f5c..acfb36f 100644 --- a/include/astra/DataProjectorPolicies.h +++ b/include/astra/DataProjectorPolicies.h @@ -319,10 +319,12 @@ class SIRTBPPolicy { CFloat32ProjectionData2D* m_pTotalRayLength; CFloat32VolumeData2D* m_pTotalPixelWeight; + float m_fRelaxation; + public: FORCEINLINE SIRTBPPolicy(); - FORCEINLINE SIRTBPPolicy(CFloat32VolumeData2D* _pReconstruction, CFloat32ProjectionData2D* _pSinogram, CFloat32VolumeData2D* _pTotalPixelWeight, CFloat32ProjectionData2D* _pTotalRayLength); + FORCEINLINE SIRTBPPolicy(CFloat32VolumeData2D* _pReconstruction, CFloat32ProjectionData2D* _pSinogram, CFloat32VolumeData2D* _pTotalPixelWeight, CFloat32ProjectionData2D* _pTotalRayLength, float _fRelaxation); FORCEINLINE ~SIRTBPPolicy(); FORCEINLINE bool rayPrior(int _iRayIndex); diff --git a/include/astra/DataProjectorPolicies.inl b/include/astra/DataProjectorPolicies.inl index 0c0eddd..f300761 100644 --- a/include/astra/DataProjectorPolicies.inl +++ b/include/astra/DataProjectorPolicies.inl @@ -712,12 +712,14 @@ SIRTBPPolicy::SIRTBPPolicy() SIRTBPPolicy::SIRTBPPolicy(CFloat32VolumeData2D* _pReconstruction, CFloat32ProjectionData2D* _pSinogram, CFloat32VolumeData2D* _pTotalPixelWeight, - CFloat32ProjectionData2D* _pTotalRayLength) + CFloat32ProjectionData2D* _pTotalRayLength, + float _fRelaxation) { m_pReconstruction = _pReconstruction; m_pSinogram = _pSinogram; m_pTotalPixelWeight = _pTotalPixelWeight; m_pTotalRayLength = _pTotalRayLength; + m_fRelaxation = _fRelaxation; } //---------------------------------------------------------------------------------------- SIRTBPPolicy::~SIRTBPPolicy() @@ -739,7 +741,7 @@ void SIRTBPPolicy::addWeight(int _iRayIndex, int _iVolumeIndex, float32 _fWeight { float32 fGammaBeta = m_pTotalPixelWeight->getData()[_iVolumeIndex] * m_pTotalRayLength->getData()[_iRayIndex]; if ((fGammaBeta > 0.001f) || (fGammaBeta < -0.001f)) { - m_pReconstruction->getData()[_iVolumeIndex] += _fWeight * m_pSinogram->getData()[_iRayIndex] / fGammaBeta; + m_pReconstruction->getData()[_iVolumeIndex] += _fWeight * m_fRelaxation * m_pSinogram->getData()[_iRayIndex] / fGammaBeta; } } //---------------------------------------------------------------------------------------- diff --git a/include/astra/PluginAlgorithm.h b/include/astra/PluginAlgorithm.h index 667e813..cbd80fc 100644 --- a/include/astra/PluginAlgorithm.h +++ b/include/astra/PluginAlgorithm.h @@ -29,62 +29,37 @@ $Id$ #ifndef _INC_ASTRA_PLUGINALGORITHM #define _INC_ASTRA_PLUGINALGORITHM -#ifdef ASTRA_PYTHON - -#include "astra/Algorithm.h" -#include "astra/Singleton.h" -#include "astra/XMLDocument.h" -#include "astra/XMLNode.h" - -// Slightly hackish forward declaration of PyObject -struct _object; -typedef _object PyObject; +#include "astra/Globals.h" +#include <map> +#include <string> namespace astra { -class _AstraExport CPluginAlgorithm : public CAlgorithm { - -public: - - CPluginAlgorithm(PyObject* pyclass); - ~CPluginAlgorithm(); - - bool initialize(const Config& _cfg); - void run(int _iNrIterations); - -private: - PyObject * instance; -}; +class CAlgorithm; -class _AstraExport CPluginAlgorithmFactory : public Singleton<CPluginAlgorithmFactory> { +class _AstraExport CPluginAlgorithmFactory { public: + CPluginAlgorithmFactory() { } + virtual ~CPluginAlgorithmFactory() { } - CPluginAlgorithmFactory(); - ~CPluginAlgorithmFactory(); + virtual CAlgorithm * getPlugin(const std::string &name) = 0; - CPluginAlgorithm * getPlugin(std::string name); + virtual bool registerPlugin(std::string name, std::string className) = 0; + virtual bool registerPlugin(std::string className) = 0; - bool registerPlugin(std::string name, std::string className); - bool registerPlugin(std::string className); - bool registerPluginClass(std::string name, PyObject * className); - bool registerPluginClass(PyObject * className); + virtual std::map<std::string, std::string> getRegisteredMap() = 0; - PyObject * getRegistered(); - std::map<std::string, std::string> getRegisteredMap(); - - std::string getHelp(std::string name); + virtual std::string getHelp(const std::string &name) = 0; + + static void registerFactory(CPluginAlgorithmFactory *factory) { m_factory = factory; } + static CPluginAlgorithmFactory* getFactory() { return m_factory; } private: - PyObject * pluginDict; - PyObject *inspect, *six; + static CPluginAlgorithmFactory *m_factory; }; -PyObject* XMLNode2dict(XMLNode node); - } #endif - -#endif diff --git a/include/astra/Projector2D.h b/include/astra/Projector2D.h index a1ea0ce..c7a899d 100644 --- a/include/astra/Projector2D.h +++ b/include/astra/Projector2D.h @@ -174,7 +174,7 @@ public: * * @return initialized successfully */ - bool isInitialized(); + bool isInitialized() const; /** get a description of the class * @@ -191,7 +191,7 @@ private: }; // inline functions -inline bool CProjector2D::isInitialized() { return m_bIsInitialized; } +inline bool CProjector2D::isInitialized() const { return m_bIsInitialized; } inline CProjectionGeometry2D* CProjector2D::getProjectionGeometry() { return m_pProjectionGeometry; } inline CVolumeGeometry2D* CProjector2D::getVolumeGeometry() { return m_pVolumeGeometry; } diff --git a/include/astra/Projector3D.h b/include/astra/Projector3D.h index 1ef1ba7..88ca2be 100644 --- a/include/astra/Projector3D.h +++ b/include/astra/Projector3D.h @@ -153,7 +153,7 @@ public: * * @return initialized successfully */ - bool isInitialized(); + bool isInitialized() const; /** get a description of the class * @@ -174,7 +174,7 @@ private: }; // inline functions -inline bool CProjector3D::isInitialized() { return m_bIsInitialized; } +inline bool CProjector3D::isInitialized() const { return m_bIsInitialized; } inline CProjectionGeometry3D* CProjector3D::getProjectionGeometry() { return m_pProjectionGeometry; } inline CVolumeGeometry3D* CProjector3D::getVolumeGeometry() { return m_pVolumeGeometry; } diff --git a/include/astra/SartAlgorithm.h b/include/astra/SartAlgorithm.h index eb4c61e..cdae029 100644 --- a/include/astra/SartAlgorithm.h +++ b/include/astra/SartAlgorithm.h @@ -49,7 +49,7 @@ namespace astra { * * The update step of pixel \f$v_j\f$ for projection \f$phi\f$ and iteration \f$k\f$ is given by: * \f[ - * v_j^{(k+1)} = v_j^{(k)} + \frac{\sum_{p_i \in P_\phi} \left( \lambda \frac{p_i - \sum_{r=1}^{N} w_{ir}v_r^{(k)}} {\sum_{r=1}^{N}w_{ir} } \right)} {\sum_{p_i \in P_\phi}w_{ij}} + * v_j^{(k+1)} = v_j^{(k)} + \lambda \frac{\sum_{p_i \in P_\phi} \left( \frac{p_i - \sum_{r=1}^{N} w_{ir}v_r^{(k)}} {\sum_{r=1}^{N}w_{ir} } \right)} {\sum_{p_i \in P_\phi}w_{ij}} * \f] * * \par XML Configuration @@ -64,6 +64,7 @@ namespace astra { * \astra_xml_item_option{MaxConstraintValue, float, 255, Maximum constraint value.} * \astra_xml_item_option{ProjectionOrder, string, "sequential", the order in which the projections are updated. 'sequential', 'random' or 'custom'} * \astra_xml_item_option{ProjectionOrderList, vector of float, not used, if ProjectionOrder='custom': use this order.} + * \astra_xml_item_option{Relaxation, float, 1, The relaxation parameter.} * * \par MATLAB example * \astra_code{ @@ -76,7 +77,8 @@ namespace astra { * cfg.option.UseMaxConstraint = 'yes';\n * cfg.option.MaxConstraintValue = 1024;\n * cfg.option.ProjectionOrder = 'custom';\n -* cfg.option.ProjectionOrderList = randperm(100);\n + * cfg.option.ProjectionOrderList = randperm(100);\n + * cfg.option.Relaxation = 1.0;\n * alg_id = astra_mex_algorithm('create'\, cfg);\n * astra_mex_algorithm('iterate'\, alg_id\, 10);\n * astra_mex_algorithm('delete'\, alg_id);\n @@ -215,6 +217,8 @@ protected: //< Current index in the projection order array. int m_iCurrentProjection; + //< Relaxation parameter + float m_fLambda; }; // inline functions diff --git a/include/astra/SirtAlgorithm.h b/include/astra/SirtAlgorithm.h index 05b3fa9..8044d09 100644 --- a/include/astra/SirtAlgorithm.h +++ b/include/astra/SirtAlgorithm.h @@ -49,7 +49,7 @@ namespace astra { * * The update step of pixel \f$v_j\f$ for iteration \f$k\f$ is given by: * \f[ - * v_j^{(k+1)} = v_j^{(k)} + \alpha \sum_{i=1}^{M} \left( \frac{w_{ij}\left( p_i - \sum_{r=1}^{N} w_{ir}v_r^{(k)}\right)}{\sum_{k=1}^{N} w_{ik}} \right) \frac{1}{\sum_{l=1}^{M}w_{lj}} + * v_j^{(k+1)} = v_j^{(k)} + \lambda \sum_{i=1}^{M} \left( \frac{w_{ij}\left( p_i - \sum_{r=1}^{N} w_{ir}v_r^{(k)}\right)}{\sum_{k=1}^{N} w_{ik}} \right) \frac{1}{\sum_{l=1}^{M}w_{lj}} * \f] * * \par XML Configuration @@ -62,6 +62,7 @@ namespace astra { * \astra_xml_item_option{MinConstraintValue, float, 0, Minimum constraint value.} * \astra_xml_item_option{UseMaxConstraint, bool, false, Use maximum value constraint.} * \astra_xml_item_option{MaxConstraintValue, float, 255, Maximum constraint value.} + * \astra_xml_item_option{Relaxation, float, 1, The relaxation factor.} * * \par XML Example * \astra_code{ @@ -74,6 +75,7 @@ namespace astra { * <Option key="UseMinConstraint" value="yes"/>\n * <Option key="UseMaxConstraint" value="yes"/>\n * <Option key="MaxConstraintValue" value="1024"/>\n + * <Option key="Relaxation" value="1"/>\n * </Algorithm> * } * @@ -88,6 +90,7 @@ namespace astra { * cfg.option.UseMinConstraint = 'yes';\n * cfg.option.UseMaxConstraint = 'yes';\n * cfg.option.MaxConstraintValue = 1024;\n + * cfg.option.Relaxation = 1.0;\n * alg_id = astra_mex_algorithm('create'\, cfg);\n * astra_mex_algorithm('iterate'\, alg_id\, 10);\n * astra_mex_algorithm('delete'\, alg_id);\n @@ -136,6 +139,10 @@ protected: */ int m_iIterationCount; + /** Relaxation parameter + */ + float m_fLambda; + public: // type of the algorithm, needed to register with CAlgorithmFactory diff --git a/include/astra/Utilities.h b/include/astra/Utilities.h index 3ae0e6c..8d7c44d 100644 --- a/include/astra/Utilities.h +++ b/include/astra/Utilities.h @@ -50,40 +50,40 @@ public: //< Parse string as int. //< Throw exception on failure. -int stringToInt(const std::string& s); +_AstraExport int stringToInt(const std::string& s); //< Parse string as float. //< Throw exception on failure. -float stringToFloat(const std::string& s); +_AstraExport float stringToFloat(const std::string& s); //< Parse string as double. //< Throw exception on failure. -double stringToDouble(const std::string& s); +_AstraExport double stringToDouble(const std::string& s); template<typename T> -T stringTo(const std::string& s); +_AstraExport T stringTo(const std::string& s); //< Parse comma/semicolon-separated string as float vector. //< Throw exception on failure. -std::vector<float> stringToFloatVector(const std::string& s); +_AstraExport std::vector<float> stringToFloatVector(const std::string& s); //< Parse comma/semicolon-separated string as double vector. //< Throw exception on failure. -std::vector<double> stringToDoubleVector(const std::string& s); +_AstraExport std::vector<double> stringToDoubleVector(const std::string& s); template<typename T> -std::vector<T> stringToVector(const std::string& s); +_AstraExport std::vector<T> stringToVector(const std::string& s); //< Generate string from float. -std::string floatToString(float f); +_AstraExport std::string floatToString(float f); //< Generate string from double. -std::string doubleToString(double f); +_AstraExport std::string doubleToString(double f); template<typename T> -std::string toString(T f); +_AstraExport std::string toString(T f); } diff --git a/matlab/mex/astra_mex_c.cpp b/matlab/mex/astra_mex_c.cpp index fdf4f33..f499528 100644 --- a/matlab/mex/astra_mex_c.cpp +++ b/matlab/mex/astra_mex_c.cpp @@ -36,10 +36,14 @@ $Id$ #include "mexInitFunctions.h" #include "astra/Globals.h" +#include "astra/AstraObjectManager.h" + #ifdef ASTRA_CUDA #include "../cuda/2d/darthelper.h" #include "astra/CompositeGeometryManager.h" #endif + + using namespace std; using namespace astra; @@ -144,10 +148,51 @@ void astra_mex_version(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[ //----------------------------------------------------------------------------------------- +void astra_mex_info(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + if (nrhs < 2) { + mexErrMsgTxt("Usage: astra_mex('info', index/indices);\n"); + return; + } + + for (int i = 1; i < nrhs; i++) { + int iDataID = (int)(mxGetScalar(prhs[i])); + CAstraObjectManagerBase *ptr; + ptr = CAstraIndexManager::getSingleton().get(iDataID); + if (ptr) { + mexPrintf("%s\t%s\n", ptr->getType().c_str(), ptr->getInfo(iDataID).c_str()); + } + } + +} + +//----------------------------------------------------------------------------------------- + +void astra_mex_delete(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + if (nrhs < 2) { + mexErrMsgTxt("Usage: astra_mex('delete', index/indices);\n"); + return; + } + + for (int i = 1; i < nrhs; i++) { + int iDataID = (int)(mxGetScalar(prhs[i])); + CAstraObjectManagerBase *ptr; + ptr = CAstraIndexManager::getSingleton().get(iDataID); + if (ptr) + ptr->remove(iDataID); + } + +} + + + +//----------------------------------------------------------------------------------------- + static void printHelp() { mexPrintf("Please specify a mode of operation.\n"); - mexPrintf(" Valid modes: version, use_cuda, credits\n"); + mexPrintf(" Valid modes: version, use_cuda, credits, set_gpu_index, info, delete\n"); } //----------------------------------------------------------------------------------------- @@ -178,6 +223,10 @@ void mexFunction(int nlhs, mxArray* plhs[], astra_mex_credits(nlhs, plhs, nrhs, prhs); } else if (sMode == std::string("set_gpu_index")) { astra_mex_set_gpu_index(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("info")) { + astra_mex_info(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("delete")) { + astra_mex_delete(nlhs, plhs, nrhs, prhs); } else { printHelp(); } diff --git a/matlab/mex/astra_mex_plugin_c.cpp b/matlab/mex/astra_mex_plugin_c.cpp index 177fcf4..4ed534e 100644 --- a/matlab/mex/astra_mex_plugin_c.cpp +++ b/matlab/mex/astra_mex_plugin_c.cpp @@ -37,9 +37,63 @@ $Id$ #include "astra/PluginAlgorithm.h" +#include <Python.h> + using namespace std; using namespace astra; +static 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); // RTLD_NOW|RTLD_DEEPBIND + if (retVal != NULL) { + PyObject *lapack = PyImport_ImportModule("numpy.linalg.lapack_lite"); + if (lapack != NULL) { + Py_DECREF(lapack); + } + PyObject *retVal2 = PyObject_CallMethod(sys, "setdlopenflags", "O",curFlags); + if (retVal2 != NULL) { + Py_DECREF(retVal2); + } + Py_DECREF(retVal); + } + Py_DECREF(curFlags); + } + Py_DECREF(sys); + } +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_plugin('init'); + * + * Initialize plugin support by initializing python and importing astra + */ +void astra_mex_plugin_init() +{ + if(!Py_IsInitialized()){ + Py_Initialize(); + PyEval_InitThreads(); + } + +#ifndef _MSC_VER + fixLapackLoading(); +#endif + + // Importing astra may be overkill, since we only need to initialize + // PythonPluginAlgorithmFactory from astra.plugin_c. + PyObject *mod = PyImport_ImportModule("astra"); + Py_XDECREF(mod); +} + //----------------------------------------------------------------------------------------- /** astra_mex_plugin('get_registered'); @@ -48,7 +102,11 @@ using namespace astra; */ void astra_mex_plugin_get_registered(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { - astra::CPluginAlgorithmFactory *fact = astra::CPluginAlgorithmFactory::getSingletonPtr(); + astra::CPluginAlgorithmFactory *fact = astra::CPluginAlgorithmFactory::getFactory(); + if (!fact) { + mexPrintf("Plugin support not initialized."); + return; + } std::map<std::string, std::string> mp = fact->getRegisteredMap(); for(std::map<std::string,std::string>::iterator it=mp.begin();it!=mp.end();it++){ mexPrintf("%s: %s\n",it->first.c_str(), it->second.c_str()); @@ -62,9 +120,13 @@ void astra_mex_plugin_get_registered(int nlhs, mxArray* plhs[], int nrhs, const */ void astra_mex_plugin_register(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + astra::CPluginAlgorithmFactory *fact = astra::CPluginAlgorithmFactory::getFactory(); + if (!fact) { + mexPrintf("Plugin support not initialized."); + return; + } if (2 <= nrhs) { string class_name = mexToString(prhs[1]); - astra::CPluginAlgorithmFactory *fact = astra::CPluginAlgorithmFactory::getSingletonPtr(); fact->registerPlugin(class_name); }else{ mexPrintf("astra_mex_plugin('register', class_name);\n"); @@ -78,9 +140,13 @@ void astra_mex_plugin_register(int nlhs, mxArray* plhs[], int nrhs, const mxArra */ void astra_mex_plugin_get_help(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { + astra::CPluginAlgorithmFactory *fact = astra::CPluginAlgorithmFactory::getFactory(); + if (!fact) { + mexPrintf("Plugin support not initialized."); + return; + } if (2 <= nrhs) { string name = mexToString(prhs[1]); - astra::CPluginAlgorithmFactory *fact = astra::CPluginAlgorithmFactory::getSingletonPtr(); mexPrintf((fact->getHelp(name)+"\n").c_str()); }else{ mexPrintf("astra_mex_plugin('get_help', name);\n"); @@ -116,12 +182,14 @@ void mexFunction(int nlhs, mxArray* plhs[], initASTRAMex(); // SWITCH (MODE) - if (sMode == std::string("get_registered")) { - astra_mex_plugin_get_registered(nlhs, plhs, nrhs, prhs); - }else if (sMode == std::string("get_help")) { - astra_mex_plugin_get_help(nlhs, plhs, nrhs, prhs); - }else if (sMode == std::string("register")) { - astra_mex_plugin_register(nlhs, plhs, nrhs, prhs); + if (sMode == "init") { + astra_mex_plugin_init(); + } else if (sMode == std::string("get_registered")) { + astra_mex_plugin_get_registered(nlhs, plhs, nrhs, prhs); + }else if (sMode == std::string("get_help")) { + astra_mex_plugin_get_help(nlhs, plhs, nrhs, prhs); + }else if (sMode == std::string("register")) { + astra_mex_plugin_register(nlhs, plhs, nrhs, prhs); } else { printHelp(); } diff --git a/matlab/mex/mexInitFunctions.cpp b/matlab/mex/mexInitFunctions.cpp index bd3df2c..7245af2 100644 --- a/matlab/mex/mexInitFunctions.cpp +++ b/matlab/mex/mexInitFunctions.cpp @@ -23,5 +23,13 @@ void initASTRAMex(){ if(!astra::CLogger::setCallbackScreen(&logCallBack)){ mexErrMsgTxt("Error initializing mex functions."); } + mexIsInitialized=true; + + + // If we have support for plugins, initialize them. + // (NB: Call this after setting mexIsInitialized, to avoid recursively + // calling initASTRAMex) + mexEvalString("if exist('astra_mex_plugin_c') == 3; astra_mex_plugin_c('init'); end"); + } diff --git a/matlab/tools/opTomo.m b/matlab/tools/opTomo.m index 71dfb1e..04b3634 100644 --- a/matlab/tools/opTomo.m +++ b/matlab/tools/opTomo.m @@ -44,11 +44,9 @@ classdef opTomo < opSpot vol_id fp_alg_id bp_alg_id + proj_id % ASTRA IDs handle astra_handle - % geometries - proj_geom; - vol_geom; end % properties properties ( SetAccess = private, GetAccess = public ) @@ -139,6 +137,17 @@ classdef opTomo < opSpot error(['Only type ' 39 'cuda' 39 ' is supported ' ... 'for 3D geometries.']) end + + % setup projector + cfg = astra_struct('cuda3d'); + cfg.ProjectionGeometry = proj_geom; + cfg.VolumeGeometry = vol_geom; + cfg.option.GPUindex = gpu_index; + + % create projector + op.proj_id = astra_mex_projector3d('create', cfg); + % create handle to ASTRA object, for cleaning up + op.astra_handle = opTomo_helper_handle(op.proj_id); % create a function handle op.funHandle = @opTomo_intrnl3D; @@ -148,8 +157,6 @@ classdef opTomo < opSpot % pass object properties op.proj_size = proj_size; op.vol_size = vol_size; - op.proj_geom = proj_geom; - op.vol_geom = vol_geom; op.cflag = false; op.sweepflag = false; @@ -169,10 +176,12 @@ classdef opTomo < opSpot if issparse(x) x = full(x); end - - % convert input to single - if isa(x, 'single') == false + + if isa(x, 'double') + isdouble = true; x = single(x); + else + isdouble = false; end % the multiplication @@ -180,6 +189,10 @@ classdef opTomo < opSpot % make sure output is column vector y = y(:); + + if isdouble + y = double(y); + end end % multiply @@ -194,7 +207,7 @@ classdef opTomo < opSpot function y = opTomo_intrnl2D(op,x,mode) if mode == 1 - % X is passed as a vector, reshape it into an image. + % x is passed as a vector, reshape it into an image. x = reshape(x, op.vol_size); % Matlab data copied to ASTRA data @@ -206,7 +219,7 @@ classdef opTomo < opSpot % retrieve Matlab array y = astra_mex_data2d('get_single', op.sino_id); else - % X is passed as a vector, reshape it into a sinogram. + % x is passed as a vector, reshape it into a sinogram. x = reshape(x, op.proj_size); % Matlab data copied to ASTRA data @@ -218,6 +231,7 @@ classdef opTomo < opSpot % retrieve Matlab array y = astra_mex_data2d('get_single', op.vol_id); end + end % opTomo_intrnl2D @@ -225,55 +239,16 @@ classdef opTomo < opSpot function y = opTomo_intrnl3D(op,x,mode) if mode == 1 - % X is passed as a vector, reshape it into an image + % x is passed as a vector, reshape it into an image x = reshape(x, op.vol_size); - % initialize output - y = zeros(op.proj_size, 'single'); - - % link matlab array to ASTRA - vol_id = astra_mex_data3d_c('link', '-vol', op.vol_geom, x, 0); - sino_id = astra_mex_data3d_c('link', '-sino', op.proj_geom, y, 1); - - % initialize fp algorithm - cfg = astra_struct('FP3D_CUDA'); - cfg.ProjectionDataId = sino_id; - cfg.VolumeDataId = vol_id; - - alg_id = astra_mex_algorithm('create', cfg); - % forward projection - astra_mex_algorithm('iterate', alg_id); - - % cleanup - astra_mex_data3d('delete', vol_id); - astra_mex_data3d('delete', sino_id); - astra_mex_algorithm('delete', alg_id); + y = astra_mex_direct('FP3D', op.proj_id, x); else - % X is passed as a vector, reshape it into projection data + % x is passed as a vector, reshape it into projection data x = reshape(x, op.proj_size); - - % initialize output - y = zeros(op.vol_size,'single'); - - % link matlab array to ASTRA - vol_id = astra_mex_data3d_c('link', '-vol', op.vol_geom, y, 1); - sino_id = astra_mex_data3d_c('link', '-sino', op.proj_geom, x, 0); - % initialize bp algorithm - cfg = astra_struct('BP3D_CUDA'); - cfg.ProjectionDataId = sino_id; - cfg.ReconstructionDataId = vol_id; - - alg_id = astra_mex_algorithm('create', cfg); - - % backprojection - astra_mex_algorithm('iterate', alg_id); - - % cleanup - astra_mex_data3d('delete', vol_id); - astra_mex_data3d('delete', sino_id); - astra_mex_algorithm('delete', alg_id); + y = astra_mex_direct('BP3D', op.proj_id, x); end end % opTomo_intrnl3D diff --git a/python/astra/PyIndexManager.pxd b/python/astra/PyIndexManager.pxd new file mode 100644 index 0000000..c1ad502 --- /dev/null +++ b/python/astra/PyIndexManager.pxd @@ -0,0 +1,40 @@ +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam +# +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox +# +# This file is part of the ASTRA Toolbox. +# +# +# The ASTRA Toolbox is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# The ASTRA Toolbox is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +# ----------------------------------------------------------------------- + +from libcpp.string cimport string + +from .PyIncludes cimport * + +cdef extern from "astra/AstraObjectManager.h" namespace "astra": + cdef cppclass CAstraObjectManagerBase: + string getInfo(int) + void remove(int) + string getType() + cdef cppclass CAstraIndexManager: + CAstraObjectManagerBase* get(int) + +cdef extern from "astra/AstraObjectManager.h" namespace "astra::CAstraIndexManager": + cdef CAstraIndexManager* getSingletonPtr() + diff --git a/python/astra/algorithm_c.pyx b/python/astra/algorithm_c.pyx index 3231c1f..4e96578 100644 --- a/python/astra/algorithm_c.pyx +++ b/python/astra/algorithm_c.pyx @@ -1,29 +1,28 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox # +# This file is part of the ASTRA Toolbox. # -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). # -#The Python interface to 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 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 Python interface to 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. +# 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 Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +# ----------------------------------------------------------------------- # -#----------------------------------------------------------------------- - # distutils: language = c++ # distutils: libraries = astra diff --git a/python/astra/astra.py b/python/astra/astra.py index 9328b6b..61c26ee 100644 --- a/python/astra/astra.py +++ b/python/astra/astra.py @@ -56,3 +56,23 @@ def set_gpu_index(idx, memory=0): :type idx: :class:`int` """ a.set_gpu_index(idx, memory) + +def delete(ids): + """Delete an astra object. + + :param ids: ID or list of ID's to delete. + :type ids: :class:`int` or :class:`list` + + """ + return a.delete(ids) + +def info(ids): + """Print info about an astra object. + + :param ids: ID or list of ID's to show. + :type ids: :class:`int` or :class:`list` + + """ + return a.info(ids) + + diff --git a/python/astra/astra_c.pyx b/python/astra/astra_c.pyx index 2a9c816..8e30e69 100644 --- a/python/astra/astra_c.pyx +++ b/python/astra/astra_c.pyx @@ -1,28 +1,28 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox # +# This file is part of the ASTRA Toolbox. # -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). # -#The Python interface to 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 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 Python interface to 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. +# 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 Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +# ----------------------------------------------------------------------- # -#----------------------------------------------------------------------- # distutils: language = c++ # distutils: libraries = astra @@ -33,6 +33,9 @@ from .utils import wrap_from_bytes from libcpp.string cimport string from libcpp.vector cimport vector from libcpp cimport bool +cimport PyIndexManager +from .PyIndexManager cimport CAstraObjectManagerBase + cdef extern from "astra/Globals.h" namespace "astra": int getVersion() string getVersionString() @@ -51,6 +54,7 @@ cdef extern from "astra/CompositeGeometryManager.h" namespace "astra": cdef extern from "astra/CompositeGeometryManager.h" namespace "astra::CCompositeGeometryManager": void setGlobalGPUParams(SGPUParams&) + def credits(): six.print_("""The ASTRA Toolbox has been developed at the University of Antwerp and CWI, Amsterdam by * Prof. dr. Joost Batenburg @@ -77,12 +81,12 @@ def version(printToScreen=False): else: return getVersion() -def set_gpu_index(idx, memory=0): - import types +IF HAVE_CUDA==True: + def set_gpu_index(idx, memory=0): import collections cdef SGPUParams params if use_cuda()==True: - if not isinstance(idx, collections.Iterable) or isinstance(idx, types.StringTypes): + if not isinstance(idx, collections.Iterable) or isinstance(idx, six.string_types + (six.text_type,six.binary_type)): idx = (idx,) params.memory = memory params.GPUIndices = idx @@ -90,3 +94,27 @@ def set_gpu_index(idx, memory=0): ret = setGPUIndex(params.GPUIndices[0]) if not ret: six.print_("Failed to set GPU " + str(params.GPUIndices[0])) +ELSE: + def set_gpu_index(idx, memory=0): + raise NotImplementedError("CUDA support is not enabled in ASTRA") + +def delete(ids): + import collections + cdef CAstraObjectManagerBase* ptr + if not isinstance(ids, collections.Iterable) or isinstance(ids, six.string_types + (six.text_type,six.binary_type)): + ids = (ids,) + for i in ids: + ptr = PyIndexManager.getSingletonPtr().get(i) + if ptr: + ptr.remove(i) + +def info(ids): + import collections + cdef CAstraObjectManagerBase* ptr + if not isinstance(ids, collections.Iterable) or isinstance(ids, six.string_types + (six.text_type,six.binary_type)): + ids = (ids,) + for i in ids: + ptr = PyIndexManager.getSingletonPtr().get(i) + if ptr: + s = ptr.getType() + six.b("\t") + ptr.getInfo(i) + six.print_(wrap_from_bytes(s)) diff --git a/python/astra/data2d_c.pyx b/python/astra/data2d_c.pyx index 801fd8e..65e80ce 100644 --- a/python/astra/data2d_c.pyx +++ b/python/astra/data2d_c.pyx @@ -1,28 +1,28 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox # +# This file is part of the ASTRA Toolbox. # -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). # -#The Python interface to 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 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 Python interface to 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. +# 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 Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +# ----------------------------------------------------------------------- # -#----------------------------------------------------------------------- # distutils: language = c++ # distutils: libraries = astra diff --git a/python/astra/data3d.py b/python/astra/data3d.py index e5ef6b0..f143659 100644 --- a/python/astra/data3d.py +++ b/python/astra/data3d.py @@ -89,7 +89,7 @@ def get_single(i): :returns: :class:`numpy.ndarray` -- The object data. """ - return g.get_single(i) + return d.get_single(i) def store(i,data): """Fill existing 3D object with data. diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx index 3b27ab7..207d9a5 100644 --- a/python/astra/data3d_c.pyx +++ b/python/astra/data3d_c.pyx @@ -1,28 +1,28 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox # +# This file is part of the ASTRA Toolbox. # -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). # -#The Python interface to 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 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 Python interface to 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. +# 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 Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +# ----------------------------------------------------------------------- # -#----------------------------------------------------------------------- # distutils: language = c++ # distutils: libraries = astra diff --git a/python/astra/experimental.pyx b/python/astra/experimental.pyx index aafc002..9bb73a2 100644 --- a/python/astra/experimental.pyx +++ b/python/astra/experimental.pyx @@ -1,6 +1,6 @@ -#----------------------------------------------------------------------- -# Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp -# 2014-2015, CWI, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # # Contact: astra@uantwerpen.be # Website: http://sf.net/projects/astra-toolbox @@ -21,8 +21,8 @@ # You should have received a copy of the GNU General Public License # along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. # -#----------------------------------------------------------------------- - +# ----------------------------------------------------------------------- +# # distutils: language = c++ # distutils: libraries = astra diff --git a/python/astra/extrautils.pyx b/python/astra/extrautils.pyx index 998f5cb..5bc315f 100644 --- a/python/astra/extrautils.pyx +++ b/python/astra/extrautils.pyx @@ -1,28 +1,27 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox # +# This file is part of the ASTRA Toolbox. # -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). # -#The Python interface to 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 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 Python interface to 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. +# 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 Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. # -#----------------------------------------------------------------------- +# ----------------------------------------------------------------------- def clipCircle(img): cdef int i,j diff --git a/python/astra/functions.py b/python/astra/functions.py index e38b5bc..3f4aa82 100644 --- a/python/astra/functions.py +++ b/python/astra/functions.py @@ -115,7 +115,7 @@ def add_noise_to_sino(sinogram_in, I0, seed=None): sinogram_out = -max_sinogramRaw * np.log(sinogramCT_D) if not isinstance(sinogram_in, np.ndarray): - at.data2d.store(sinogram_in, sinogram_out) + data2d.store(sinogram_in, sinogram_out) if not seed==None: np.random.set_state(curstate) diff --git a/python/astra/log_c.pyx b/python/astra/log_c.pyx index 55c63e6..0d187e9 100644 --- a/python/astra/log_c.pyx +++ b/python/astra/log_c.pyx @@ -1,28 +1,28 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox # +# This file is part of the ASTRA Toolbox. # -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). # -#The Python interface to 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 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 Python interface to 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. +# 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 Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +# ----------------------------------------------------------------------- # -#----------------------------------------------------------------------- # distutils: language = c++ # distutils: libraries = astra diff --git a/python/astra/matrix_c.pyx b/python/astra/matrix_c.pyx index d099a75..f5c0938 100644 --- a/python/astra/matrix_c.pyx +++ b/python/astra/matrix_c.pyx @@ -1,28 +1,28 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox # +# This file is part of the ASTRA Toolbox. # -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). # -#The Python interface to 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 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 Python interface to 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. +# 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 Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +# ----------------------------------------------------------------------- # -#----------------------------------------------------------------------- # distutils: language = c++ # distutils: libraries = astra diff --git a/python/astra/optomo.py b/python/astra/optomo.py index 4a64150..dd10713 100644 --- a/python/astra/optomo.py +++ b/python/astra/optomo.py @@ -160,7 +160,7 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): return self._matvec(v) return scipy.sparse.linalg.LinearOperator.__mul__(self, v) - def reconstruct(self, method, s, iterations=1, extraOptions = {}): + def reconstruct(self, method, s, iterations=1, extraOptions = None): """Reconstruct an object. :param method: Method to use for reconstruction. @@ -172,6 +172,8 @@ class OpTomo(scipy.sparse.linalg.LinearOperator): :param extraOptions: Extra options to use during reconstruction (i.e. for cfg['option']). :type extraOptions: :class:`dict` """ + if extraOptions is None: + extraOptions={} s = self.__checkArray(s, self.sshape) sid = self.data_mod.link('-sino',self.pg,s) v = np.zeros(self.vshape,dtype=np.float32) diff --git a/python/astra/plugin_c.pyx b/python/astra/plugin_c.pyx index 8d6816b..ee04853 100644 --- a/python/astra/plugin_c.pyx +++ b/python/astra/plugin_c.pyx @@ -1,28 +1,28 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox # +# This file is part of the ASTRA Toolbox. # -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). # -#The Python interface to 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 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 Python interface to 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. +# 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 Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +# ----------------------------------------------------------------------- # -#----------------------------------------------------------------------- # distutils: language = c++ # distutils: libraries = astra @@ -32,21 +32,25 @@ import inspect from libcpp.string cimport string from libcpp cimport bool -cdef CPluginAlgorithmFactory *fact = getSingletonPtr() +cdef CPythonPluginAlgorithmFactory *fact = getSingletonPtr() from . import utils -cdef extern from "astra/PluginAlgorithm.h" namespace "astra": - cdef cppclass CPluginAlgorithmFactory: +cdef extern from "src/PythonPluginAlgorithm.h" namespace "astra": + cdef cppclass CPythonPluginAlgorithmFactory: bool registerPlugin(string className) bool registerPlugin(string name, string className) bool registerPluginClass(object className) bool registerPluginClass(string name, object className) object getRegistered() - string getHelp(string name) + string getHelp(string &name) + +cdef extern from "src/PythonPluginAlgorithm.h" namespace "astra::CPythonPluginAlgorithmFactory": + cdef CPythonPluginAlgorithmFactory* getSingletonPtr() cdef extern from "astra/PluginAlgorithm.h" namespace "astra::CPluginAlgorithmFactory": - cdef CPluginAlgorithmFactory* getSingletonPtr() + # NB: Using wrong pointer type here for convenience + cdef void registerFactory(CPythonPluginAlgorithmFactory *) def register(className, name=None): if inspect.isclass(className): @@ -65,3 +69,6 @@ def get_registered(): def get_help(name): return utils.wrap_from_bytes(fact.getHelp(six.b(name))) + +# Register python plugin factory with astra +registerFactory(fact) diff --git a/python/astra/projector3d_c.pyx b/python/astra/projector3d_c.pyx index aec9cde..e355e38 100644 --- a/python/astra/projector3d_c.pyx +++ b/python/astra/projector3d_c.pyx @@ -1,28 +1,28 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox # +# This file is part of the ASTRA Toolbox. # -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). # -#The Python interface to 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 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 Python interface to 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. +# 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 Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +# ----------------------------------------------------------------------- # -#----------------------------------------------------------------------- # distutils: language = c++ # distutils: libraries = astra diff --git a/python/astra/projector_c.pyx b/python/astra/projector_c.pyx index 77c64a4..53d38c3 100644 --- a/python/astra/projector_c.pyx +++ b/python/astra/projector_c.pyx @@ -1,28 +1,28 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox # +# This file is part of the ASTRA Toolbox. # -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). # -#The Python interface to 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 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 Python interface to 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. +# 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 Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +# ----------------------------------------------------------------------- # -#----------------------------------------------------------------------- # distutils: language = c++ # distutils: libraries = astra diff --git a/python/astra/src/PythonPluginAlgorithm.cpp b/python/astra/src/PythonPluginAlgorithm.cpp new file mode 100644 index 0000000..617c0f4 --- /dev/null +++ b/python/astra/src/PythonPluginAlgorithm.cpp @@ -0,0 +1,372 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp + 2014-2016, CWI, Amsterdam + +Contact: astra@uantwerpen.be +Website: http://sf.net/projects/astra-toolbox + +This file is part of the ASTRA Toolbox. + + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. + +----------------------------------------------------------------------- +*/ + +#ifdef ASTRA_PYTHON + +#include "PythonPluginAlgorithm.h" + +#include "astra/Logging.h" +#include "astra/Utilities.h" +#include <boost/algorithm/string.hpp> +#include <boost/algorithm/string/split.hpp> +#include <iostream> +#include <fstream> +#include <string> + +#include <Python.h> +#include "bytesobject.h" + +namespace astra { + + + +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); +} + +CPythonPluginAlgorithmFactory::CPythonPluginAlgorithmFactory(){ + if(!Py_IsInitialized()){ + Py_Initialize(); + PyEval_InitThreads(); + } + pluginDict = PyDict_New(); + inspect = PyImport_ImportModule("inspect"); + six = PyImport_ImportModule("six"); +} + +CPythonPluginAlgorithmFactory::~CPythonPluginAlgorithmFactory(){ + 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<std::string> 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<items.size();i++){ + submod = PyObject_GetAttrString(submod,items[i].c_str()); + Py_DECREF(pyclass); + pyclass = submod; + if(pyclass==NULL){ + logPythonError(); + return NULL; + } + } + return pyclass; +} + +bool CPythonPluginAlgorithmFactory::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 CPythonPluginAlgorithmFactory::registerPlugin(std::string className){ + PyObject *pyclass = getClassFromString(className); + if(pyclass==NULL) return false; + bool ret = registerPluginClass(pyclass); + Py_DECREF(pyclass); + return ret; +} + +bool CPythonPluginAlgorithmFactory::registerPluginClass(std::string name, PyObject * className){ + PyDict_SetItemString(pluginDict, name.c_str(), className); + return true; +} + +bool CPythonPluginAlgorithmFactory::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; +} + +CAlgorithm * CPythonPluginAlgorithmFactory::getPlugin(const std::string &name){ + PyObject *className = PyDict_GetItemString(pluginDict, name.c_str()); + if(className==NULL) return NULL; + CPluginAlgorithm *alg = NULL; + if(PyBytes_Check(className)){ + std::string str = std::string(PyBytes_AsString(className)); + PyObject *pyclass = getClassFromString(str); + if(pyclass!=NULL){ + alg = new CPluginAlgorithm(pyclass); + Py_DECREF(pyclass); + } + }else{ + alg = new CPluginAlgorithm(className); + } + return alg; +} + +PyObject * CPythonPluginAlgorithmFactory::getRegistered(){ + Py_INCREF(pluginDict); + return pluginDict; +} + +std::map<std::string, std::string> CPythonPluginAlgorithmFactory::getRegisteredMap(){ + std::map<std::string, std::string> 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 CPythonPluginAlgorithmFactory::getHelp(const 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(CPythonPluginAlgorithmFactory); + +#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<std::string> rows, row; + boost::split(rows, str, boost::is_any_of(";")); + PyObject *mat = PyList_New(rows.size()); + for(unsigned int i=0; i<rows.size(); i++){ + boost::split(row, rows[i], boost::is_any_of(",")); + PyObject *rowlist = PyList_New(row.size()); + for(unsigned int j=0;j<row.size();j++){ + PyList_SetItem(rowlist, j, PyFloat_FromDouble(StringUtil::stringToDouble(row[j]))); + } + PyList_SetItem(mat, i, rowlist); + } + return mat; + } + if(str.find(",")!=std::string::npos){ + std::vector<std::string> vec; + boost::split(vec, str, boost::is_any_of(",")); + PyObject *veclist = PyList_New(vec.size()); + for(unsigned int i=0;i<vec.size();i++){ + PyList_SetItem(veclist, i, PyFloat_FromDouble(StringUtil::stringToDouble(vec[i]))); + } + return veclist; + } + try{ + return PyLong_FromLong(StringUtil::stringToInt(str)); + }catch(const StringUtil::bad_cast &){ + try{ + return PyFloat_FromDouble(StringUtil::stringToDouble(str)); + }catch(const StringUtil::bad_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<XMLNode> nodes = node.getNodes(); + std::list<XMLNode>::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 diff --git a/python/astra/src/PythonPluginAlgorithm.h b/python/astra/src/PythonPluginAlgorithm.h new file mode 100644 index 0000000..ea4c6fb --- /dev/null +++ b/python/astra/src/PythonPluginAlgorithm.h @@ -0,0 +1,88 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp + 2014-2016, CWI, Amsterdam + +Contact: astra@uantwerpen.be +Website: http://sf.net/projects/astra-toolbox + +This file is part of the ASTRA Toolbox. + + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. + +----------------------------------------------------------------------- +*/ + +#ifndef _INC_PYTHONPLUGINALGORITHM +#define _INC_PYTHONPLUGINALGORITHM + +#ifdef ASTRA_PYTHON + +#include "astra/Algorithm.h" +#include "astra/Singleton.h" +#include "astra/XMLDocument.h" +#include "astra/XMLNode.h" +#include "astra/PluginAlgorithm.h" + +#include <Python.h> + +namespace astra { +class CPluginAlgorithm : public CAlgorithm { + +public: + + CPluginAlgorithm(PyObject* pyclass); + ~CPluginAlgorithm(); + + bool initialize(const Config& _cfg); + void run(int _iNrIterations); + +private: + PyObject * instance; + +}; + +class CPythonPluginAlgorithmFactory : public CPluginAlgorithmFactory, public Singleton<CPythonPluginAlgorithmFactory> { + +public: + + CPythonPluginAlgorithmFactory(); + virtual ~CPythonPluginAlgorithmFactory(); + + virtual CAlgorithm * getPlugin(const std::string &name); + + virtual bool registerPlugin(std::string name, std::string className); + virtual bool registerPlugin(std::string className); + bool registerPluginClass(std::string name, PyObject * className); + bool registerPluginClass(PyObject * className); + + PyObject * getRegistered(); + virtual std::map<std::string, std::string> getRegisteredMap(); + + virtual std::string getHelp(const std::string &name); + +private: + PyObject * pluginDict; + PyObject *inspect, *six; +}; + +PyObject* XMLNode2dict(XMLNode node); + +} + + +#endif + +#endif diff --git a/python/astra/utils.pyx b/python/astra/utils.pyx index 07727ce..34d1902 100644 --- a/python/astra/utils.pyx +++ b/python/astra/utils.pyx @@ -1,28 +1,28 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox # +# This file is part of the ASTRA Toolbox. # -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). # -#The Python interface to 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 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 Python interface to 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. +# 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 Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# +# ----------------------------------------------------------------------- # -#----------------------------------------------------------------------- # distutils: language = c++ # distutils: libraries = astra @@ -32,7 +32,7 @@ import six if six.PY3: import builtins else: - import __builtin__ + import __builtin__ as builtins from libcpp.string cimport string from libcpp.vector cimport vector from libcpp.list cimport list @@ -95,14 +95,14 @@ cdef void readDict(XMLNode root, _dc): dc = convert_item(_dc) for item in dc: val = dc[item] - if isinstance(val, __builtins__.list) or isinstance(val, tuple): + if isinstance(val, builtins.list) or isinstance(val, tuple): val = np.array(val,dtype=np.float64) if isinstance(val, np.ndarray): if val.size == 0: break listbase = root.addChildNode(item) contig_data = np.ascontiguousarray(val,dtype=np.float64) - data = <double*>np.PyArray_DATA(contig_data) + data = <double*>np.PyArray_DATA(contig_data) if val.ndim == 2: listbase.setContent(data, val.shape[1], val.shape[0], False) elif val.ndim == 1: @@ -119,6 +119,8 @@ cdef void readDict(XMLNode root, _dc): if item == six.b('type'): root.addAttribute(< string > six.b('type'), <string> wrap_to_bytes(val)) else: + if isinstance(val, builtins.bool): + val = int(val) itm = root.addChildNode(item, wrap_to_bytes(val)) cdef void readOptions(XMLNode node, dc): @@ -131,7 +133,7 @@ cdef void readOptions(XMLNode node, dc): val = dc[item] if node.hasOption(item): raise Exception('Duplicate Option: %s' % item) - if isinstance(val, __builtins__.list) or isinstance(val, tuple): + if isinstance(val, builtins.list) or isinstance(val, tuple): val = np.array(val,dtype=np.float64) if isinstance(val, np.ndarray): if val.size == 0: @@ -147,6 +149,8 @@ cdef void readOptions(XMLNode node, dc): else: raise Exception("Only 1 or 2 dimensions are allowed") else: + if isinstance(val, builtins.bool): + val = int(val) node.addOption(item, wrap_to_bytes(val)) cdef configToDict(Config *cfg): diff --git a/python/builder.py b/python/builder.py index 018b26b..dcd62d8 100644 --- a/python/builder.py +++ b/python/builder.py @@ -1,28 +1,28 @@ -#----------------------------------------------------------------------- -#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam +# ----------------------------------------------------------------------- +# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp +# 2013-2016, CWI, Amsterdam # -#Author: Daniel M. Pelt -#Contact: D.M.Pelt@cwi.nl -#Website: http://dmpelt.github.io/pyastratoolbox/ +# Contact: astra@uantwerpen.be +# Website: http://sf.net/projects/astra-toolbox # +# This file is part of the ASTRA Toolbox. # -#This file is part of the Python interface to the -#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). # -#The Python interface to 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 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 Python interface to 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. +# 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 Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. +# You should have received a copy of the GNU General Public License +# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. # #----------------------------------------------------------------------- + import sys import os import numpy as np @@ -70,12 +70,16 @@ ext_modules = [ ] ext_modules = cythonize("astra/*.pyx", language_level=2) cmdclass = { 'build_ext': build_ext } +for m in ext_modules: + if m.name == 'astra.plugin_c': + m.sources.append('astra/src/PythonPluginAlgorithm.cpp') + setup (name = 'PyASTRAToolbox', version = '1.7.1', description = 'Python interface to the ASTRA-Toolbox', author='D.M. Pelt', author_email='D.M.Pelt@cwi.nl', - url='http://dmpelt.github.io/pyastratoolbox/', + url='http://sf.net/projects/astra-toolbox', #ext_package='astra', #ext_modules = cythonize(Extension("astra/*.pyx",extra_compile_args=extra_compile_args,extra_linker_args=extra_compile_args)), license='GPLv3', diff --git a/python/conda/build.sh b/python/conda/build.sh index 814ea7e..13ae3f8 100644 --- a/python/conda/build.sh +++ b/python/conda/build.sh @@ -5,12 +5,4 @@ if [ $MAKEOPTS == '<UNDEFINED>' ] then MAKEOPTS="" fi -make $MAKEOPTS install-libraries -make $MAKEOPTS python-root-install -LIBPATH=lib -if [ $ARCH == 64 ] - then - LIBPATH+=64 -fi -cp -P $CUDA_ROOT/$LIBPATH/libcudart.so.* $PREFIX/lib -cp -P $CUDA_ROOT/$LIBPATH/libcufft.so.* $PREFIX/lib +make $MAKEOPTS python-root-install
\ No newline at end of file diff --git a/python/conda/libastra/build.sh b/python/conda/libastra/build.sh new file mode 100644 index 0000000..e1d9700 --- /dev/null +++ b/python/conda/libastra/build.sh @@ -0,0 +1,15 @@ +cd build/linux +./autogen.sh +./configure --with-cuda=$CUDA_ROOT --prefix=$PREFIX +if [ $MAKEOPTS == '<UNDEFINED>' ] + then + MAKEOPTS="" +fi +make $MAKEOPTS install-libraries +LIBPATH=lib +if [ $ARCH == 64 ] + then + LIBPATH+=64 +fi +cp -P $CUDA_ROOT/$LIBPATH/libcudart.so.* $PREFIX/lib +cp -P $CUDA_ROOT/$LIBPATH/libcufft.so.* $PREFIX/lib diff --git a/python/conda/libastra/meta.yaml b/python/conda/libastra/meta.yaml new file mode 100644 index 0000000..73fa0d7 --- /dev/null +++ b/python/conda/libastra/meta.yaml @@ -0,0 +1,22 @@ +package: + name: libastra + version: '1.8b' + +source: + git_url: https://github.com/astra-toolbox/astra-toolbox.git + #git_tag: v1.7.1 # Change to 1.8 after release + +build: + number: 0 + script_env: + - CUDA_ROOT + - MAKEOPTS + +about: + home: http://www.astra-toolbox.com + license: GPLv3 + summary: 'The ASTRA Toolbox is a Python toolbox of high-performance GPU primitives for 2D and 3D tomography.' + +# See +# http://docs.continuum.io/conda/build.html for +# more information about meta.yaml diff --git a/python/conda/meta.yaml b/python/conda/meta.yaml index 7e4679b..e6a7f52 100644 --- a/python/conda/meta.yaml +++ b/python/conda/meta.yaml @@ -1,10 +1,10 @@ package: name: astra-toolbox - version: '1.7.1' + version: '1.8b' source: git_url: https://github.com/astra-toolbox/astra-toolbox.git - git_tag: v1.7.1 + #git_tag: v1.7.1 # Change to 1.8 after release build: number: 0 @@ -29,10 +29,11 @@ requirements: - numpy - scipy - six + - libastra ==1.8b about: - home: http://sourceforge.net/p/astra-toolbox/wiki/Home/ + home: http://www.astra-toolbox.com license: GPLv3 summary: 'The ASTRA Toolbox is a Python toolbox of high-performance GPU primitives for 2D and 3D tomography.' diff --git a/src/ArtAlgorithm.cpp b/src/ArtAlgorithm.cpp index b59bd93..526c263 100644 --- a/src/ArtAlgorithm.cpp +++ b/src/ArtAlgorithm.cpp @@ -156,8 +156,12 @@ bool CArtAlgorithm::initialize(const Config& _cfg) return false; } + // "Lambda" is replaced by the more descriptive "Relaxation" m_fLambda = _cfg.self.getOptionNumerical("Lambda", 1.0f); - CC.markOptionParsed("Lambda"); + m_fLambda = _cfg.self.getOptionNumerical("Relaxation", m_fLambda); + if (!_cfg.self.hasOption("Relaxation")) + CC.markOptionParsed("Lambda"); + CC.markOptionParsed("Relaxation"); // success m_bIsInitialized = _check(); @@ -232,7 +236,7 @@ map<string,boost::any> CArtAlgorithm::getInformation() { map<string, boost::any> res; res["RayOrder"] = getInformation("RayOrder"); - res["Lambda"] = getInformation("Lambda"); + res["Relaxation"] = getInformation("Relaxation"); return mergeMap<string,boost::any>(CReconstructionAlgorithm2D::getInformation(), res); }; @@ -240,7 +244,7 @@ map<string,boost::any> CArtAlgorithm::getInformation() // Information - Specific boost::any CArtAlgorithm::getInformation(std::string _sIdentifier) { - if (_sIdentifier == "Lambda") { return m_fLambda; } + if (_sIdentifier == "Relaxation") { return m_fLambda; } if (_sIdentifier == "RayOrder") { vector<float32> res; for (int i = 0; i < m_iRayCount; i++) { diff --git a/src/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<CProjector2D>); -DEFINE_SINGLETON(CAstraObjectManager<CProjector3D>); -DEFINE_SINGLETON(CAstraObjectManager<CFloat32Data2D>); -DEFINE_SINGLETON(CAstraObjectManager<CFloat32Data3D>); -DEFINE_SINGLETON(CAstraObjectManager<CAlgorithm>); -DEFINE_SINGLETON(CAstraObjectManager<CSparseMatrix>); +DEFINE_SINGLETON(CProjector2DManager); +DEFINE_SINGLETON(CProjector3DManager); +DEFINE_SINGLETON(CData2DManager); +DEFINE_SINGLETON(CData3DManager); +DEFINE_SINGLETON(CAlgorithmManager); +DEFINE_SINGLETON(CMatrixManager); + +DEFINE_SINGLETON(CAstraIndexManager); } // end namespace diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp index 96b28e9..084ba8c 100644 --- a/src/CompositeGeometryManager.cpp +++ b/src/CompositeGeometryManager.cpp @@ -45,14 +45,15 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. #include <cstring> #include <sstream> +#include <climits> #ifndef USE_PTHREADS #include <boost/thread/mutex.hpp> #include <boost/thread.hpp> #endif -namespace astra { +namespace astra { SGPUParams* CCompositeGeometryManager::s_params = 0; @@ -98,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) @@ -111,7 +115,22 @@ 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, UINT_MAX, div); +#if 0 + TPartList splitOutput2; + for (TPartList::iterator i_out = splitOutput.begin(); i_out != splitOutput.end(); ++i_out) { + boost::shared_ptr<CPart> outputPart = *i_out; + 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<CPart> outputPart = *i_out; + outputPart.get()->splitY(splitOutput, UINT_MAX, UINT_MAX, 1); + } + splitOutput2.clear(); +#endif + for (TJobList::const_iterator j = L.begin(); j != L.end(); ++j) { @@ -139,8 +158,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, maxBlockDim, 1); delete input; + TPartList splitInput2; + for (TPartList::iterator i_in = splitInput.begin(); i_in != splitInput.end(); ++i_in) { + boost::shared_ptr<CPart> inputPart = *i_in; + 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<CPart> inputPart = *i_in; + inputPart.get()->splitY(splitInput, UINT_MAX, maxBlockDim, 1); + } + splitInput2.clear(); + ASTRA_DEBUG("Input split into %d parts", splitInput.size()); for (TPartList::iterator i_in = splitInput.begin(); @@ -327,10 +359,14 @@ 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); + 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); @@ -410,7 +446,17 @@ SPar3DProjection* getProjectionVectors(const CParallelVecProjectionGeometry3D* p template<class V> -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<class V> +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 +466,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 -static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry3D* pProjGeom, int v, int size) + const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(pProjGeom); + const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(pProjGeom); + const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(pProjGeom); + const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(pProjGeom); + + if (conegeom || conevec3dgeom) { + SConeProjection* pConeProjs; + if (conegeom) { + pConeProjs = getProjectionVectors<SConeProjection>(conegeom); + } else { + pConeProjs = getProjectionVectors<SConeProjection>(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<SPar3DProjection>(par3dgeom); + } else { + pParProjs = getProjectionVectors<SPar3DProjection>(parvec3dgeom); + } + + translateProjectionVectorsU(pParProjs, pProjGeom->getProjectionCount(), u); + + CProjectionGeometry3D* ret = new CParallelVecProjectionGeometry3D(pProjGeom->getProjectionCount(), + pProjGeom->getDetectorRowCount(), + size, + pParProjs); + + delete[] pParProjs; + return ret; + } + +} + + + +static CProjectionGeometry3D* getSubProjectionGeometryV(const CProjectionGeometry3D* pProjGeom, int v, int size) { // First convert to vectors, then translate, then convert into new object @@ -438,7 +534,7 @@ static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry pConeProjs = getProjectionVectors<SConeProjection>(conevec3dgeom); } - translateProjectionVectors(pConeProjs, pProjGeom->getProjectionCount(), v); + translateProjectionVectorsV(pConeProjs, pProjGeom->getProjectionCount(), v); CProjectionGeometry3D* ret = new CConeVecProjectionGeometry3D(pProjGeom->getProjectionCount(), size, @@ -457,7 +553,7 @@ static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry pParProjs = getProjectionVectors<SPar3DProjection>(parvec3dgeom); } - translateProjectionVectors(pParProjs, pProjGeom->getProjectionCount(), v); + translateProjectionVectorsV(pParProjs, pProjGeom->getProjectionCount(), v); CProjectionGeometry3D* ret = new CParallelVecProjectionGeometry3D(pProjGeom->getProjectionCount(), size, @@ -476,17 +572,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<CPart>(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<CPart>(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 +708,9 @@ CCompositeGeometryManager::TPartList CCompositeGeometryManager::CVolumePart::spl pGeom->getWindowMaxY(), pGeom->getWindowMinZ() + shift + size * pGeom->getPixelLengthZ()); - ret.push_back(boost::shared_ptr<CPart>(sub)); + out.push_back(boost::shared_ptr<CPart>(sub)); } } - - return ret; } CCompositeGeometryManager::CVolumePart* CCompositeGeometryManager::CVolumePart::clone() const @@ -630,7 +817,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 +826,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<CPart>(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<CPart>(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 +897,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<CPart>(sub)); + out.push_back(boost::shared_ptr<CPart>(sub)); } } - return ret; - } CCompositeGeometryManager::CProjectionPart* CCompositeGeometryManager::CProjectionPart::clone() const diff --git a/src/CudaDataOperationAlgorithm.cpp b/src/CudaDataOperationAlgorithm.cpp index 15886a4..82b676b 100644 --- a/src/CudaDataOperationAlgorithm.cpp +++ b/src/CudaDataOperationAlgorithm.cpp @@ -76,7 +76,7 @@ bool CCudaDataOperationAlgorithm::initialize(const Config& _cfg) node = _cfg.self.getSingleNode("DataId"); ASTRA_CONFIG_CHECK(node, "CCudaDataOperationAlgorithm", "No DataId tag specified."); vector<string> data = node.getContentArray(); - for (vector<string>::iterator it = data.begin(); it != data.end(); it++){ + for (vector<string>::iterator it = data.begin(); it != data.end(); ++it){ int id = StringUtil::stringToInt(*it); m_pData.push_back(dynamic_cast<CFloat32Data2D*>(CData2DManager::getSingleton().get(id))); } diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp index 5a1910c..2798434 100644 --- a/src/CudaReconstructionAlgorithm2D.cpp +++ b/src/CudaReconstructionAlgorithm2D.cpp @@ -329,6 +329,20 @@ bool CCudaReconstructionAlgorithm2D::setupGeometry() } //---------------------------------------------------------------------------------------- + +void CCudaReconstructionAlgorithm2D::initCUDAAlgorithm() +{ + bool ok; + + ok = setupGeometry(); + ASTRA_ASSERT(ok); + + ok = m_pAlgo->allocateBuffers(); + ASTRA_ASSERT(ok); +} + + +//---------------------------------------------------------------------------------------- // Iterate void CCudaReconstructionAlgorithm2D::run(int _iNrIterations) { @@ -339,13 +353,7 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations) const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); if (!m_bAlgoInit) { - - ok = setupGeometry(); - ASTRA_ASSERT(ok); - - ok = m_pAlgo->allocateBuffers(); - ASTRA_ASSERT(ok); - + initCUDAAlgorithm(); m_bAlgoInit = true; } diff --git a/src/CudaSartAlgorithm.cpp b/src/CudaSartAlgorithm.cpp index d202847..bf97224 100644 --- a/src/CudaSartAlgorithm.cpp +++ b/src/CudaSartAlgorithm.cpp @@ -107,7 +107,8 @@ bool CCudaSartAlgorithm::initialize(const Config& _cfg) CC.markOptionParsed("ProjectionOrderList"); } - + m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f); + CC.markOptionParsed("Relaxation"); return true; } @@ -123,12 +124,26 @@ bool CCudaSartAlgorithm::initialize(CProjector2D* _pProjector, if (!m_bIsInitialized) return false; + m_fLambda = 1.0f; + m_pAlgo = new astraCUDA::SART(); m_bAlgoInit = false; return true; } +//---------------------------------------------------------------------------------------- + +void CCudaSartAlgorithm::initCUDAAlgorithm() +{ + CCudaReconstructionAlgorithm2D::initCUDAAlgorithm(); + + astraCUDA::SART* pSart = dynamic_cast<astraCUDA::SART*>(m_pAlgo); + + pSart->setRelaxation(m_fLambda); +} + + } // namespace astra diff --git a/src/CudaSirtAlgorithm.cpp b/src/CudaSirtAlgorithm.cpp index 33e381a..c8dc677 100644 --- a/src/CudaSirtAlgorithm.cpp +++ b/src/CudaSirtAlgorithm.cpp @@ -50,6 +50,8 @@ CCudaSirtAlgorithm::CCudaSirtAlgorithm() m_pMinMask = 0; m_pMaxMask = 0; + + m_fLambda = 1.0f; } //---------------------------------------------------------------------------------------- @@ -86,6 +88,8 @@ bool CCudaSirtAlgorithm::initialize(const Config& _cfg) } CC.markOptionParsed("MaxMaskId"); + m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f); + CC.markOptionParsed("Relaxation"); m_pAlgo = new astraCUDA::SIRT(); m_bAlgoInit = false; @@ -108,41 +112,30 @@ bool CCudaSirtAlgorithm::initialize(CProjector2D* _pProjector, m_pAlgo = new astraCUDA::SIRT(); m_bAlgoInit = false; + m_fLambda = 1.0f; return true; } //---------------------------------------------------------------------------------------- -// Iterate -void CCudaSirtAlgorithm::run(int _iNrIterations) + +void CCudaSirtAlgorithm::initCUDAAlgorithm() { - // check initialized - ASTRA_ASSERT(m_bIsInitialized); + CCudaReconstructionAlgorithm2D::initCUDAAlgorithm(); - if (!m_bAlgoInit) { - // We only override the initialisation step to copy the min/max masks + astraCUDA::SIRT* pSirt = dynamic_cast<astraCUDA::SIRT*>(m_pAlgo); - bool ok = setupGeometry(); + if (m_pMinMask || m_pMaxMask) { + const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); + const float *pfMinMaskData = 0; + const float *pfMaxMaskData = 0; + if (m_pMinMask) pfMinMaskData = m_pMinMask->getDataConst(); + if (m_pMaxMask) pfMaxMaskData = m_pMaxMask->getDataConst(); + bool ok = pSirt->uploadMinMaxMasks(pfMinMaskData, pfMaxMaskData, volgeom.getGridColCount()); ASTRA_ASSERT(ok); - - ok = m_pAlgo->allocateBuffers(); - ASTRA_ASSERT(ok); - - if (m_pMinMask || m_pMaxMask) { - const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); - astraCUDA::SIRT* pSirt = dynamic_cast<astraCUDA::SIRT*>(m_pAlgo); - const float *pfMinMaskData = 0; - const float *pfMaxMaskData = 0; - if (m_pMinMask) pfMinMaskData = m_pMinMask->getDataConst(); - if (m_pMaxMask) pfMaxMaskData = m_pMaxMask->getDataConst(); - ok = pSirt->uploadMinMaxMasks(pfMinMaskData, pfMaxMaskData, volgeom.getGridColCount()); - ASTRA_ASSERT(ok); - } - - m_bAlgoInit = true; } - CCudaReconstructionAlgorithm2D::run(_iNrIterations); + pSirt->setRelaxation(m_fLambda); } diff --git a/src/CudaSirtAlgorithm3D.cpp b/src/CudaSirtAlgorithm3D.cpp index 605c470..c819f8e 100644 --- a/src/CudaSirtAlgorithm3D.cpp +++ b/src/CudaSirtAlgorithm3D.cpp @@ -56,6 +56,7 @@ CCudaSirtAlgorithm3D::CCudaSirtAlgorithm3D() m_iGPUIndex = -1; m_iVoxelSuperSampling = 1; m_iDetectorSuperSampling = 1; + m_fLambda = 1.0f; } //---------------------------------------------------------------------------------------- @@ -128,6 +129,8 @@ bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg) return false; } + m_fLambda = _cfg.self.getOptionNumerical("Relaxation"); + initializeFromProjector(); // Deprecated options @@ -135,6 +138,7 @@ bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg) m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling); m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex); m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex); + CC.markOptionParsed("VoxelSuperSampling"); CC.markOptionParsed("DetectorSuperSampling"); CC.markOptionParsed("GPUIndex"); @@ -164,6 +168,8 @@ bool CCudaSirtAlgorithm3D::initialize(CProjector3D* _pProjector, clear(); } + m_fLambda = 1.0f; + // required classes m_pProjector = _pProjector; m_pSinogram = _pSinogram; @@ -224,6 +230,8 @@ void CCudaSirtAlgorithm3D::run(int _iNrIterations) ASTRA_ASSERT(ok); + m_pSirt->setRelaxation(m_fLambda); + m_bAstraSIRTInit = true; } diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp index c195578..ccbfec6 100644 --- a/src/FilteredBackProjectionAlgorithm.cpp +++ b/src/FilteredBackProjectionAlgorithm.cpp @@ -117,12 +117,10 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) int angleCount = projectionIndex.size(); int detectorCount = m_pProjector->getProjectionGeometry()->getDetectorCount(); + // TODO: There is no need to allocate this. Better just + // create the CFloat32ProjectionData2D object directly, and use its + // memory. float32 * sinogramData2D = new float32[angleCount* detectorCount]; - float32 ** sinogramData = new float32*[angleCount]; - for (int i = 0; i < angleCount; i++) - { - sinogramData[i] = &(sinogramData2D[i * detectorCount]); - } float32 * projectionAngles = new float32[angleCount]; float32 detectorWidth = m_pProjector->getProjectionGeometry()->getDetectorWidth(); @@ -130,6 +128,8 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) for (int i = 0; i < angleCount; i ++) { if (projectionIndex[i] > m_pProjector->getProjectionGeometry()->getProjectionAngleCount() -1 ) { + delete[] sinogramData2D; + delete[] projectionAngles; ASTRA_ERROR("Invalid Projection Index"); return false; } else { @@ -139,7 +139,6 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) { sinogramData2D[i*detectorCount+ iDetector] = m_pSinogram->getData2D()[orgIndex][iDetector]; } -// sinogramData[i] = m_pSinogram->getSingleProjectionData(projectionIndex[i]); projectionAngles[i] = m_pProjector->getProjectionGeometry()->getProjectionAngle((int)projectionIndex[i] ); } @@ -148,6 +147,9 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg) CParallelProjectionGeometry2D * pg = new CParallelProjectionGeometry2D(angleCount, detectorCount,detectorWidth,projectionAngles); m_pProjector = new CParallelBeamLineKernelProjector2D(pg,m_pReconstruction->getGeometry()); m_pSinogram = new CFloat32ProjectionData2D(pg, sinogramData2D); + + delete[] sinogramData2D; + delete[] projectionAngles; } // TODO: check that the angles are linearly spaced between 0 and pi diff --git a/src/Float32VolumeData3DMemory.cpp b/src/Float32VolumeData3DMemory.cpp index af45cb9..14adb1a 100644 --- a/src/Float32VolumeData3DMemory.cpp +++ b/src/Float32VolumeData3DMemory.cpp @@ -136,7 +136,6 @@ CFloat32VolumeData2D * CFloat32VolumeData3DMemory::fetchSliceZ(int _iSliceIndex) CFloat32VolumeData2D* res = new CFloat32VolumeData2D(&volGeom); // copy data - int iSliceCount = m_pGeometry->getGridSliceCount(); float * pfTargetData = res->getData(); for(int iRowIndex = 0; iRowIndex < iRowCount; iRowIndex++) { diff --git a/src/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 <http://www.gnu.org/licenses/>. $Id$ */ -#ifdef ASTRA_PYTHON - #include "astra/PluginAlgorithm.h" -#include "astra/Logging.h" -#include "astra/Utilities.h" -#include <boost/algorithm/string.hpp> -#include <boost/algorithm/string/split.hpp> -#include <iostream> -#include <fstream> -#include <string> - -#include <Python.h> -#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<std::string> 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<items.size();i++){ - submod = PyObject_GetAttrString(submod,items[i].c_str()); - Py_DECREF(pyclass); - pyclass = submod; - if(pyclass==NULL){ - logPythonError(); - return NULL; - } - } - 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; - CPluginAlgorithm *alg = NULL; - if(PyBytes_Check(className)){ - std::string str = std::string(PyBytes_AsString(className)); - PyObject *pyclass = getClassFromString(str); - if(pyclass!=NULL){ - alg = new CPluginAlgorithm(pyclass); - Py_DECREF(pyclass); - } - }else{ - alg = new CPluginAlgorithm(className); - } - return alg; -} - -PyObject * CPluginAlgorithmFactory::getRegistered(){ - Py_INCREF(pluginDict); - return pluginDict; -} - -std::map<std::string, std::string> CPluginAlgorithmFactory::getRegisteredMap(){ - std::map<std::string, std::string> 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<std::string> rows, row; - boost::split(rows, str, boost::is_any_of(";")); - PyObject *mat = PyList_New(rows.size()); - for(unsigned int i=0; i<rows.size(); i++){ - boost::split(row, rows[i], boost::is_any_of(",")); - PyObject *rowlist = PyList_New(row.size()); - for(unsigned int j=0;j<row.size();j++){ - PyList_SetItem(rowlist, j, PyFloat_FromDouble(StringUtil::stringToDouble(row[j]))); - } - PyList_SetItem(mat, i, rowlist); - } - return mat; - } - if(str.find(",")!=std::string::npos){ - std::vector<std::string> vec; - boost::split(vec, str, boost::is_any_of(",")); - PyObject *veclist = PyList_New(vec.size()); - for(unsigned int i=0;i<vec.size();i++){ - PyList_SetItem(veclist, i, PyFloat_FromDouble(StringUtil::stringToDouble(vec[i]))); - } - return veclist; - } - try{ - return PyLong_FromLong(StringUtil::stringToInt(str)); - }catch(const StringUtil::bad_cast &){ - try{ - return PyFloat_FromDouble(StringUtil::stringToDouble(str)); - }catch(const StringUtil::bad_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<XMLNode> nodes = node.getNodes(); - std::list<XMLNode>::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 diff --git a/src/SartAlgorithm.cpp b/src/SartAlgorithm.cpp index 9346160..8df3342 100644 --- a/src/SartAlgorithm.cpp +++ b/src/SartAlgorithm.cpp @@ -151,6 +151,9 @@ bool CSartAlgorithm::initialize(const Config& _cfg) CC.markOptionParsed("ProjectionOrderList"); } + m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f); + CC.markOptionParsed("Relaxation"); + // create data objects m_pTotalRayLength = new CFloat32ProjectionData2D(m_pProjector->getProjectionGeometry()); m_pTotalPixelWeight = new CFloat32VolumeData2D(m_pProjector->getVolumeGeometry()); @@ -246,6 +249,7 @@ map<string,boost::any> CSartAlgorithm::getInformation() { map<string, boost::any> res; res["ProjectionOrder"] = getInformation("ProjectionOrder"); + res["Relaxation"] = getInformation("Relaxation"); return mergeMap<string,boost::any>(CReconstructionAlgorithm2D::getInformation(), res); }; @@ -253,6 +257,8 @@ map<string,boost::any> CSartAlgorithm::getInformation() // Information - Specific boost::any CSartAlgorithm::getInformation(std::string _sIdentifier) { + if (_sIdentifier == "Relaxation") + return m_fLambda; if (_sIdentifier == "ProjectionOrder") { vector<float32> res; for (int i = 0; i < m_iProjectionCount; i++) { @@ -272,9 +278,8 @@ void CSartAlgorithm::run(int _iNrIterations) m_bShouldAbort = false; - int iIteration = 0; - // data projectors + CDataProjectorInterface* pFirstForwardProjector; CDataProjectorInterface* pForwardProjector; CDataProjectorInterface* pBackProjector; @@ -286,13 +291,13 @@ void CSartAlgorithm::run(int _iNrIterations) m_pProjector, SinogramMaskPolicy(m_pSinogramMask), // sinogram mask ReconstructionMaskPolicy(m_pReconstructionMask), // reconstruction mask - SIRTBPPolicy(m_pReconstruction, m_pDiffSinogram, m_pTotalPixelWeight, m_pTotalRayLength), // SIRT backprojection + SIRTBPPolicy(m_pReconstruction, m_pDiffSinogram, m_pTotalPixelWeight, m_pTotalRayLength, m_fLambda), // SIRT backprojection m_bUseSinogramMask, m_bUseReconstructionMask, true // options on/off ); // first time forward projection data projector, // also computes total pixel weight and total ray length - pForwardProjector = dispatchDataProjector( + pFirstForwardProjector = dispatchDataProjector( m_pProjector, SinogramMaskPolicy(m_pSinogramMask), // sinogram mask ReconstructionMaskPolicy(m_pReconstructionMask), // reconstruction mask @@ -303,16 +308,30 @@ void CSartAlgorithm::run(int _iNrIterations) m_bUseSinogramMask, m_bUseReconstructionMask, true // options on/off ); + // forward projection data projector + pForwardProjector = dispatchDataProjector( + m_pProjector, + SinogramMaskPolicy(m_pSinogramMask), // sinogram mask + ReconstructionMaskPolicy(m_pReconstructionMask), // reconstruction mask + CombinePolicy<DiffFPPolicy, TotalPixelWeightPolicy>( // 2 basic operations + DiffFPPolicy(m_pReconstruction, m_pDiffSinogram, m_pSinogram), // forward projection with difference calculation + TotalPixelWeightPolicy(m_pTotalPixelWeight)), // calculate the total pixel weights + m_bUseSinogramMask, m_bUseReconstructionMask, true // options on/off + ); + // iteration loop - for (; iIteration < _iNrIterations && !m_bShouldAbort; ++iIteration) { + for (int iIteration = 0; iIteration < _iNrIterations && !m_bShouldAbort; ++iIteration) { int iProjection = m_piProjectionOrder[m_iIterationCount % m_iProjectionCount]; // forward projection and difference calculation m_pTotalPixelWeight->setData(0.0f); - pForwardProjector->projectSingleProjection(iProjection); + if (iIteration < m_iProjectionCount) + pFirstForwardProjector->projectSingleProjection(iProjection); + else + pForwardProjector->projectSingleProjection(iProjection); // backprojection pBackProjector->projectSingleProjection(iProjection); // update iteration count @@ -325,6 +344,7 @@ void CSartAlgorithm::run(int _iNrIterations) } + ASTRA_DELETE(pFirstForwardProjector); ASTRA_DELETE(pForwardProjector); ASTRA_DELETE(pBackProjector); diff --git a/src/SirtAlgorithm.cpp b/src/SirtAlgorithm.cpp index d9f3a65..ff25648 100644 --- a/src/SirtAlgorithm.cpp +++ b/src/SirtAlgorithm.cpp @@ -76,6 +76,7 @@ void CSirtAlgorithm::_clear() m_pDiffSinogram = NULL; m_pTmpVolume = NULL; + m_fLambda = 1.0f; m_iIterationCount = 0; } @@ -91,6 +92,7 @@ void CSirtAlgorithm::clear() ASTRA_DELETE(m_pDiffSinogram); ASTRA_DELETE(m_pTmpVolume); + m_fLambda = 1.0f; m_iIterationCount = 0; } @@ -128,6 +130,9 @@ bool CSirtAlgorithm::initialize(const Config& _cfg) return false; } + m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f); + CC.markOptionParsed("Relaxation"); + // init data objects and data projectors _init(); @@ -152,6 +157,8 @@ bool CSirtAlgorithm::initialize(CProjector2D* _pProjector, m_pSinogram = _pSinogram; m_pReconstruction = _pReconstruction; + m_fLambda = 1.0f; + // init data objects and data projectors _init(); @@ -248,7 +255,7 @@ void CSirtAlgorithm::run(int _iNrIterations) x = 1.0f / x; else x = 0.0f; - pfT[i] = x; + pfT[i] = m_fLambda * x; } pfT = m_pTotalRayLength->getData(); for (int i = 0; i < m_pTotalRayLength->getSize(); ++i) { @@ -296,7 +303,7 @@ void CSirtAlgorithm::run(int _iNrIterations) m_pTmpVolume->setData(0.0f); pBackProjector->project(); - // divide by pixel weights + // multiply with relaxation factor divided by pixel weights (*m_pTmpVolume) *= (*m_pTotalPixelWeight); (*m_pReconstruction) += (*m_pTmpVolume); diff --git a/src/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; diff --git a/src/XMLNode.cpp b/src/XMLNode.cpp index 40a9b22..cf268c2 100644 --- a/src/XMLNode.cpp +++ b/src/XMLNode.cpp @@ -158,7 +158,7 @@ vector<string> XMLNode::getContentArray() const vector<string> res(iSize); // loop all list item nodes list<XMLNode> nodes = getNodes("ListItem"); - for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); it++) { + for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); ++it) { int iIndex = it->getAttributeNumerical("index"); string sValue = it->getAttribute("value"); ASTRA_ASSERT(iIndex < iSize); @@ -290,7 +290,7 @@ vector<float32> XMLNode::getOptionNumericalArray(string _sKey) const if (!hasOption(_sKey)) return vector<float32>(); list<XMLNode> nodes = getNodes("Option"); - for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); it++) { + for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); ++it) { if (it->getAttribute("key") == _sKey) { vector<float32> vals = it->getContentNumericalArray(); return vals; |