summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2016-04-18 11:51:32 +0200
committerWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2016-04-18 11:51:32 +0200
commitc366f2b07ce16c4ccdafc7cc4199fdac2d3ffef2 (patch)
treee9dc99f16607021d60ddb24c63ba7438e41a235d
parentb8ee38bdada2067f4351b27d841e68580bcbff8e (diff)
parent547def0ea6e3eab07b7e4c48cee6d6a81f6155e1 (diff)
downloadastra-c366f2b07ce16c4ccdafc7cc4199fdac2d3ffef2.tar.gz
astra-c366f2b07ce16c4ccdafc7cc4199fdac2d3ffef2.tar.bz2
astra-c366f2b07ce16c4ccdafc7cc4199fdac2d3ffef2.tar.xz
astra-c366f2b07ce16c4ccdafc7cc4199fdac2d3ffef2.zip
Merge branch 'master' into aniso
-rw-r--r--.travis.yml45
-rw-r--r--build/linux/Makefile.in22
-rw-r--r--build/msvc/gen.py7
-rw-r--r--cuda/2d/sart.cu7
-rw-r--r--cuda/2d/sart.h4
-rw-r--r--cuda/2d/sirt.cu8
-rw-r--r--cuda/2d/sirt.h4
-rw-r--r--cuda/3d/astra3d.cu13
-rw-r--r--cuda/3d/astra3d.h2
-rw-r--r--cuda/3d/mem3d.cu19
-rw-r--r--cuda/3d/mem3d.h1
-rw-r--r--cuda/3d/sirt3d.cu8
-rw-r--r--cuda/3d/sirt3d.h5
-rw-r--r--include/astra/Algorithm.h4
-rw-r--r--include/astra/ArtAlgorithm.h4
-rw-r--r--include/astra/AstraObjectFactory.h7
-rw-r--r--include/astra/AstraObjectManager.h121
-rw-r--r--include/astra/CompositeGeometryManager.h12
-rw-r--r--include/astra/CudaReconstructionAlgorithm2D.h3
-rw-r--r--include/astra/CudaSartAlgorithm.h10
-rw-r--r--include/astra/CudaSirtAlgorithm.h10
-rw-r--r--include/astra/CudaSirtAlgorithm3D.h3
-rw-r--r--include/astra/DataProjectorPolicies.h4
-rw-r--r--include/astra/DataProjectorPolicies.inl6
-rw-r--r--include/astra/PluginAlgorithm.h57
-rw-r--r--include/astra/Projector2D.h4
-rw-r--r--include/astra/Projector3D.h4
-rw-r--r--include/astra/SartAlgorithm.h8
-rw-r--r--include/astra/SirtAlgorithm.h9
-rw-r--r--include/astra/Utilities.h20
-rw-r--r--matlab/mex/astra_mex_c.cpp51
-rw-r--r--matlab/mex/astra_mex_plugin_c.cpp86
-rw-r--r--matlab/mex/mexInitFunctions.cpp8
-rw-r--r--matlab/tools/opTomo.m81
-rw-r--r--python/astra/PyIndexManager.pxd40
-rw-r--r--python/astra/algorithm_c.pyx37
-rw-r--r--python/astra/astra.py20
-rw-r--r--python/astra/astra_c.pyx70
-rw-r--r--python/astra/data2d_c.pyx36
-rw-r--r--python/astra/data3d.py2
-rw-r--r--python/astra/data3d_c.pyx36
-rw-r--r--python/astra/experimental.pyx10
-rw-r--r--python/astra/extrautils.pyx35
-rw-r--r--python/astra/functions.py2
-rw-r--r--python/astra/log_c.pyx36
-rw-r--r--python/astra/matrix_c.pyx36
-rw-r--r--python/astra/optomo.py4
-rw-r--r--python/astra/plugin_c.pyx53
-rw-r--r--python/astra/projector3d_c.pyx36
-rw-r--r--python/astra/projector_c.pyx36
-rw-r--r--python/astra/src/PythonPluginAlgorithm.cpp372
-rw-r--r--python/astra/src/PythonPluginAlgorithm.h88
-rw-r--r--python/astra/utils.pyx48
-rw-r--r--python/builder.py40
-rw-r--r--python/conda/build.sh10
-rw-r--r--python/conda/libastra/build.sh15
-rw-r--r--python/conda/libastra/meta.yaml22
-rw-r--r--python/conda/meta.yaml7
-rw-r--r--src/ArtAlgorithm.cpp10
-rw-r--r--src/AstraObjectManager.cpp16
-rw-r--r--src/CompositeGeometryManager.cpp272
-rw-r--r--src/CudaDataOperationAlgorithm.cpp2
-rw-r--r--src/CudaReconstructionAlgorithm2D.cpp22
-rw-r--r--src/CudaSartAlgorithm.cpp17
-rw-r--r--src/CudaSirtAlgorithm.cpp41
-rw-r--r--src/CudaSirtAlgorithm3D.cpp8
-rw-r--r--src/FilteredBackProjectionAlgorithm.cpp14
-rw-r--r--src/Float32VolumeData3DMemory.cpp1
-rw-r--r--src/PluginAlgorithm.cpp367
-rw-r--r--src/SartAlgorithm.cpp32
-rw-r--r--src/SirtAlgorithm.cpp11
-rw-r--r--src/Utilities.cpp2
-rw-r--r--src/XMLNode.cpp4
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 {
* &lt;Option key="UseMinConstraint" value="yes"/&gt;\n
* &lt;Option key="UseMaxConstraint" value="yes"/&gt;\n
* &lt;Option key="MaxConstraintValue" value="1024"/&gt;\n
+ * &lt;Option key="Relaxation" value="1"/&gt;\n
* &lt;/Algorithm&gt;
* }
*
@@ -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;