summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--build/linux/Makefile.in6
-rw-r--r--cuda/2d/fft.cu46
-rw-r--r--include/astra/Singleton.h28
-rw-r--r--include/astra/Utilities.h3
-rw-r--r--matlab/mex/mexHelpFunctions.cpp28
-rw-r--r--python/astra/data3d_c.pyx2
-rw-r--r--python/astra/extrautils.pyx2
-rw-r--r--python/astra/optomo.py96
-rw-r--r--python/astra/src/PythonPluginAlgorithm.cpp10
-rw-r--r--samples/python/s009_projection_matrix.py2
-rw-r--r--samples/python/s015_fp_bp.py6
-rw-r--r--samples/python/s017_OpTomo.py2
-rw-r--r--samples/python/s018_plugin.py34
-rw-r--r--src/CompositeGeometryManager.cpp321
-rw-r--r--src/Utilities.cpp34
15 files changed, 378 insertions, 242 deletions
diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in
index a199bf6..f10f482 100644
--- a/build/linux/Makefile.in
+++ b/build/linux/Makefile.in
@@ -281,13 +281,15 @@ ifeq ($(python),yes)
py: libastra.la
$(MKDIR) python/build
$(MKDIR) python/finalbuild
- cd $(srcdir)/../../python; CPPFLAGS="${PYCPPFLAGS}" LDFLAGS="${PYLDFLAGS}" $(PYTHON) builder.py build --build-base=$(abs_top_builddir)/python/build install \
+ # Note: setting CC to CXX is intentional. Python uses CC for compilation even if input is C++.
+ cd $(srcdir)/../../python; CXX="${CXX}" CC="${CXX}" CPPFLAGS="${PYCPPFLAGS}" LDFLAGS="${PYLDFLAGS}" $(PYTHON) builder.py build --build-base=$(abs_top_builddir)/python/build install \
--install-base=$(abs_top_builddir)/python/finalbuild --install-headers=$(abs_top_builddir)/python/finalbuild --install-purelib=$(abs_top_builddir)/python/finalbuild \
--install-platlib=$(abs_top_builddir)/python/finalbuild --install-scripts=$(abs_top_builddir)/python/finalbuild --install-data=$(abs_top_builddir)/python/finalbuild
python-root-install: libastra.la
$(MKDIR) python/build
- cd $(srcdir)/../../python; CPPFLAGS="${PYCPPFLAGS}" LDFLAGS="${PYLDFLAGS}" $(PYTHON) builder.py build --build-base=$(abs_top_builddir)/python/build install
+ # Note: setting CC to CXX is intentional. Python uses CC for compilation even if input is C++.
+ cd $(srcdir)/../../python; CXX="${CXX}" CC="${CXX}" CPPFLAGS="${PYCPPFLAGS}" LDFLAGS="${PYLDFLAGS}" $(PYTHON) builder.py build --build-base=$(abs_top_builddir)/python/build install
endif
diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu
index 2bfd493..2d259a9 100644
--- a/cuda/2d/fft.cu
+++ b/cuda/2d/fft.cu
@@ -35,6 +35,7 @@ $Id$
#include <fstream>
#include "../../include/astra/Logging.h"
+#include "astra/Fourier.h"
using namespace astra;
@@ -303,16 +304,48 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
float * pfFilt = new float[_iFFTFourierDetectorCount];
float * pfW = new float[_iFFTFourierDetectorCount];
+ // We cache one Fourier transform for repeated FBP's of the same size
+ static float *pfData = 0;
+ static int iFilterCacheSize = 0;
+
+ if (!pfData || iFilterCacheSize != _iFFTRealDetectorCount) {
+ // Compute filter in spatial domain
+
+ delete[] pfData;
+ pfData = new float[2*_iFFTRealDetectorCount];
+ int *ip = new int[int(2+sqrt(_iFFTRealDetectorCount)+1)];
+ ip[0] = 0;
+ float32 *w = new float32[_iFFTRealDetectorCount/2];
+
+ for (int i = 0; i < _iFFTRealDetectorCount; ++i) {
+ pfData[2*i+1] = 0.0f;
+
+ if (i & 1) {
+ int j = i;
+ if (2*j > _iFFTRealDetectorCount)
+ j = _iFFTRealDetectorCount - j;
+ float f = M_PI * j;
+ pfData[2*i] = -1 / (f*f);
+ } else {
+ pfData[2*i] = 0.0f;
+ }
+ }
+
+ pfData[0] = 0.25f;
+
+ cdft(2*_iFFTRealDetectorCount, -1, pfData, ip, w);
+ delete[] ip;
+ delete[] w;
+
+ iFilterCacheSize = _iFFTRealDetectorCount;
+ }
+
for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
{
float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount;
- // filt = 2*( 0:(order/2) )./order;
- pfFilt[iDetectorIndex] = 2.0f * fRelIndex;
- //pfFilt[iDetectorIndex] = 1.0f;
-
- // w = 2*pi*(0:size(filt,2)-1)/order
- pfW[iDetectorIndex] = 3.1415f * 2.0f * fRelIndex;
+ pfFilt[iDetectorIndex] = 2.0f * pfData[2*iDetectorIndex];
+ pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex;
}
switch(_eFilter)
@@ -866,5 +899,4 @@ void downloadDebugFilterReal(float * _pfHostSinogram, int _iProjectionCount,
free(pfHostFilter);
}
-
#endif
diff --git a/include/astra/Singleton.h b/include/astra/Singleton.h
index a256187..9d3c088 100644
--- a/include/astra/Singleton.h
+++ b/include/astra/Singleton.h
@@ -45,11 +45,7 @@ class Singleton {
public:
// constructor
- Singleton() {
- assert(!m_singleton);
- int offset = (uintptr_t)(T*)1 - (uintptr_t)(Singleton<T>*)(T*)1;
- m_singleton = (T*)((uintptr_t)this + offset);
- };
+ Singleton() { }
// destructor
virtual ~Singleton() {
@@ -57,15 +53,17 @@ class Singleton {
m_singleton = 0;
}
+ static void construct();
+
// get singleton
static T& getSingleton() {
if (!m_singleton)
- m_singleton = new T();
+ construct();
return *m_singleton;
}
static T* getSingletonPtr() {
if (!m_singleton)
- m_singleton = new T();
+ construct();
return m_singleton;
}
@@ -76,11 +74,23 @@ class Singleton {
};
-#define DEFINE_SINGLETON(T) template<> T* Singleton<T >::m_singleton = 0
+// We specifically avoid defining construct() in the header.
+// That way, the call to new is always executed by code inside libastra.
+// This avoids the situation where a singleton gets created by a copy
+// of the constructor linked into an object file outside of libastra, such
+// as a .mex file, which would then also cause the vtable to be outside of
+// libastra. This situation would cause issues when .mex files are unloaded.
+
+#define DEFINE_SINGLETON(T) \
+template<> void Singleton<T >::construct() { assert(!m_singleton); m_singleton = new T(); } \
+template<> T* Singleton<T >::m_singleton = 0
+
// This is a hack to support statements like
// DEFINE_SINGLETON2(CTemplatedClass<C1, C2>);
-#define DEFINE_SINGLETON2(A,B) template<> A,B* Singleton<A,B >::m_singleton = 0
+#define DEFINE_SINGLETON2(A,B) \
+template<> void Singleton<A,B >::construct() { assert(!m_singleton); m_singleton = new A,B(); } \
+template<> A,B* Singleton<A,B >::m_singleton = 0
} // end namespace
diff --git a/include/astra/Utilities.h b/include/astra/Utilities.h
index 8d7c44d..22adfe2 100644
--- a/include/astra/Utilities.h
+++ b/include/astra/Utilities.h
@@ -85,6 +85,9 @@ _AstraExport std::string doubleToString(double f);
template<typename T>
_AstraExport std::string toString(T f);
+
+_AstraExport void splitString(std::vector<std::string> &items, const std::string& s, const char *delim);
+
}
diff --git a/matlab/mex/mexHelpFunctions.cpp b/matlab/mex/mexHelpFunctions.cpp
index 13c4ade..d957aea 100644
--- a/matlab/mex/mexHelpFunctions.cpp
+++ b/matlab/mex/mexHelpFunctions.cpp
@@ -33,11 +33,6 @@ $Id$
#include "mexHelpFunctions.h"
#include "astra/Utilities.h"
-#include <algorithm>
-#include <boost/algorithm/string.hpp>
-#include <boost/algorithm/string/split.hpp>
-#include <boost/algorithm/string/classification.hpp>
-
using namespace std;
using namespace astra;
@@ -362,8 +357,8 @@ mxArray* stringToMxArray(std::string input)
// split rows
std::vector<std::string> row_strings;
std::vector<std::string> col_strings;
- boost::split(row_strings, input, boost::is_any_of(";"));
- boost::split(col_strings, row_strings[0], boost::is_any_of(","));
+ StringUtil::splitString(row_strings, input, ";");
+ StringUtil::splitString(col_strings, row_strings[0], ",");
// get dimensions
int rows = row_strings.size();
@@ -375,7 +370,7 @@ mxArray* stringToMxArray(std::string input)
// loop elements
for (unsigned int row = 0; row < rows; row++) {
- boost::split(col_strings, row_strings[row], boost::is_any_of(","));
+ StringUtil::splitString(col_strings, row_strings[row], ",");
// check size
for (unsigned int col = 0; col < col_strings.size(); col++) {
out[col*rows + row] = StringUtil::stringToFloat(col_strings[col]);
@@ -389,7 +384,7 @@ mxArray* stringToMxArray(std::string input)
// split
std::vector<std::string> items;
- boost::split(items, input, boost::is_any_of(","));
+ StringUtil::splitString(items, input, ",");
// init matrix
mxArray* pVector = mxCreateDoubleMatrix(1, items.size(), mxREAL);
@@ -402,16 +397,13 @@ mxArray* stringToMxArray(std::string input)
return pVector;
}
- // number
- char* end;
- double content = ::strtod(input.c_str(), &end);
- bool isnumber = !*end;
- if (isnumber) {
- return mxCreateDoubleScalar(content);
+ try {
+ // number
+ return mxCreateDoubleScalar(StringUtil::stringToDouble(input));
+ } catch (const StringUtil::bad_cast &) {
+ // string
+ return mxCreateString(input.c_str());
}
-
- // string
- return mxCreateString(input.c_str());
}
//-----------------------------------------------------------------------------------------
// turn a c++ map into a matlab struct
diff --git a/python/astra/data3d_c.pyx b/python/astra/data3d_c.pyx
index 207d9a5..811d1e4 100644
--- a/python/astra/data3d_c.pyx
+++ b/python/astra/data3d_c.pyx
@@ -264,7 +264,7 @@ def store(i,data):
def dimensions(i):
cdef CFloat32Data3D * pDataObject = getObject(i)
- return (pDataObject.getWidth(),pDataObject.getHeight(),pDataObject.getDepth())
+ return (pDataObject.getDepth(),pDataObject.getHeight(),pDataObject.getWidth())
def delete(ids):
try:
diff --git a/python/astra/extrautils.pyx b/python/astra/extrautils.pyx
index 5bc315f..2c7771e 100644
--- a/python/astra/extrautils.pyx
+++ b/python/astra/extrautils.pyx
@@ -22,6 +22,8 @@
# along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#
# -----------------------------------------------------------------------
+# distutils: language = c++
+
def clipCircle(img):
cdef int i,j
diff --git a/python/astra/optomo.py b/python/astra/optomo.py
index dd10713..dde719e 100644
--- a/python/astra/optomo.py
+++ b/python/astra/optomo.py
@@ -111,21 +111,7 @@ class OpTomo(scipy.sparse.linalg.LinearOperator):
:param v: Volume to forward project.
:type v: :class:`numpy.ndarray`
"""
- v = self.__checkArray(v, self.vshape)
- vid = self.data_mod.link('-vol',self.vg,v)
- s = np.zeros(self.sshape,dtype=np.float32)
- sid = self.data_mod.link('-sino',self.pg,s)
-
- cfg = creators.astra_dict('FP'+self.appendString)
- cfg['ProjectionDataId'] = sid
- cfg['VolumeDataId'] = vid
- cfg['ProjectorId'] = self.proj_id
- fp_id = algorithm.create(cfg)
- algorithm.run(fp_id)
-
- algorithm.delete(fp_id)
- self.data_mod.delete([vid,sid])
- return s.flatten()
+ return self.FP(v, out=None).ravel()
def rmatvec(self,s):
"""Implements the transpose operator.
@@ -133,21 +119,7 @@ class OpTomo(scipy.sparse.linalg.LinearOperator):
:param s: The projection data.
:type s: :class:`numpy.ndarray`
"""
- s = self.__checkArray(s, self.sshape)
- sid = self.data_mod.link('-sino',self.pg,s)
- v = np.zeros(self.vshape,dtype=np.float32)
- vid = self.data_mod.link('-vol',self.vg,v)
-
- cfg = creators.astra_dict('BP'+self.appendString)
- cfg['ProjectionDataId'] = sid
- cfg['ReconstructionDataId'] = vid
- cfg['ProjectorId'] = self.proj_id
- bp_id = algorithm.create(cfg)
- algorithm.run(bp_id)
-
- algorithm.delete(bp_id)
- self.data_mod.delete([vid,sid])
- return v.flatten()
+ return self.BP(s, out=None).ravel()
def __mul__(self,v):
"""Provides easy forward operator by *.
@@ -189,6 +161,70 @@ class OpTomo(scipy.sparse.linalg.LinearOperator):
self.data_mod.delete([vid,sid])
return v
+ def FP(self,v,out=None):
+ """Perform forward projection.
+
+ Output must have the right 2D/3D shape. Input may also be flattened.
+
+ Output must also be contiguous and float32. This isn't required for the
+ input, but it is more efficient if it is.
+
+ :param v: Volume to forward project.
+ :type v: :class:`numpy.ndarray`
+ :param out: Array to store result in.
+ :type out: :class:`numpy.ndarray`
+ """
+
+ v = self.__checkArray(v, self.vshape)
+ vid = self.data_mod.link('-vol',self.vg,v)
+ if out is None:
+ out = np.zeros(self.sshape,dtype=np.float32)
+ sid = self.data_mod.link('-sino',self.pg,out)
+
+ cfg = creators.astra_dict('FP'+self.appendString)
+ cfg['ProjectionDataId'] = sid
+ cfg['VolumeDataId'] = vid
+ cfg['ProjectorId'] = self.proj_id
+ fp_id = algorithm.create(cfg)
+ algorithm.run(fp_id)
+
+ algorithm.delete(fp_id)
+ self.data_mod.delete([vid,sid])
+ return out
+
+ def BP(self,s,out=None):
+ """Perform backprojection.
+
+ Output must have the right 2D/3D shape. Input may also be flattened.
+
+ Output must also be contiguous and float32. This isn't required for the
+ input, but it is more efficient if it is.
+
+ :param : The projection data.
+ :type s: :class:`numpy.ndarray`
+ :param out: Array to store result in.
+ :type out: :class:`numpy.ndarray`
+ """
+ s = self.__checkArray(s, self.sshape)
+ sid = self.data_mod.link('-sino',self.pg,s)
+ if out is None:
+ out = np.zeros(self.vshape,dtype=np.float32)
+ vid = self.data_mod.link('-vol',self.vg,out)
+
+ cfg = creators.astra_dict('BP'+self.appendString)
+ cfg['ProjectionDataId'] = sid
+ cfg['ReconstructionDataId'] = vid
+ cfg['ProjectorId'] = self.proj_id
+ bp_id = algorithm.create(cfg)
+ algorithm.run(bp_id)
+
+ algorithm.delete(bp_id)
+ self.data_mod.delete([vid,sid])
+ return out
+
+
+
+
class OpTomoTranspose(scipy.sparse.linalg.LinearOperator):
"""This object provides the transpose operation (``.T``) of the OpTomo object.
diff --git a/python/astra/src/PythonPluginAlgorithm.cpp b/python/astra/src/PythonPluginAlgorithm.cpp
index 617c0f4..893db94 100644
--- a/python/astra/src/PythonPluginAlgorithm.cpp
+++ b/python/astra/src/PythonPluginAlgorithm.cpp
@@ -31,8 +31,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
#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>
@@ -146,7 +144,7 @@ CPythonPluginAlgorithmFactory::~CPythonPluginAlgorithmFactory(){
PyObject * getClassFromString(std::string str){
std::vector<std::string> items;
- boost::split(items, str, boost::is_any_of("."));
+ StringUtil::splitString(items, str, ".");
PyObject *pyclass = PyImport_ImportModule(items[0].c_str());
if(pyclass==NULL){
logPythonError();
@@ -303,10 +301,10 @@ PyObject * pyStringFromString(std::string str){
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(";"));
+ StringUtil::splitString(rows, str, ";");
PyObject *mat = PyList_New(rows.size());
for(unsigned int i=0; i<rows.size(); i++){
- boost::split(row, rows[i], boost::is_any_of(","));
+ StringUtil::splitString(row, rows[i], ",");
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])));
@@ -317,7 +315,7 @@ PyObject* stringToPythonValue(std::string str){
}
if(str.find(",")!=std::string::npos){
std::vector<std::string> vec;
- boost::split(vec, str, boost::is_any_of(","));
+ StringUtil::splitString(vec, str, ",");
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])));
diff --git a/samples/python/s009_projection_matrix.py b/samples/python/s009_projection_matrix.py
index c4c4557..e20d58c 100644
--- a/samples/python/s009_projection_matrix.py
+++ b/samples/python/s009_projection_matrix.py
@@ -46,7 +46,7 @@ W = astra.matrix.get(matrix_id)
# Manually use this projection matrix to do a projection:
import scipy.io
P = scipy.io.loadmat('phantom.mat')['phantom256']
-s = W.dot(P.flatten())
+s = W.dot(P.ravel())
s = np.reshape(s, (len(proj_geom['ProjectionAngles']),proj_geom['DetectorCount']))
import pylab
diff --git a/samples/python/s015_fp_bp.py b/samples/python/s015_fp_bp.py
index fa0bf86..ff0b30a 100644
--- a/samples/python/s015_fp_bp.py
+++ b/samples/python/s015_fp_bp.py
@@ -46,12 +46,12 @@ class astra_wrap(object):
def matvec(self,v):
sid, s = astra.create_sino(np.reshape(v,(vol_geom['GridRowCount'],vol_geom['GridColCount'])),self.proj_id)
astra.data2d.delete(sid)
- return s.flatten()
+ return s.ravel()
def rmatvec(self,v):
bid, b = astra.create_backprojection(np.reshape(v,(len(proj_geom['ProjectionAngles']),proj_geom['DetectorCount'],)),self.proj_id)
astra.data2d.delete(bid)
- return b.flatten()
+ return b.ravel()
vol_geom = astra.create_vol_geom(256, 256)
proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False))
@@ -65,7 +65,7 @@ proj_id = astra.create_projector('cuda',proj_geom,vol_geom)
sinogram_id, sinogram = astra.create_sino(P, proj_id)
# Reshape the sinogram into a vector
-b = sinogram.flatten()
+b = sinogram.ravel()
# Call lsqr with ASTRA FP and BP
import scipy.sparse.linalg
diff --git a/samples/python/s017_OpTomo.py b/samples/python/s017_OpTomo.py
index 967fa64..214e9a7 100644
--- a/samples/python/s017_OpTomo.py
+++ b/samples/python/s017_OpTomo.py
@@ -50,7 +50,7 @@ pylab.figure(2)
pylab.imshow(sinogram)
# Run the lsqr linear solver
-output = scipy.sparse.linalg.lsqr(W, sinogram.flatten(), iter_lim=150)
+output = scipy.sparse.linalg.lsqr(W, sinogram.ravel(), iter_lim=150)
rec = output[0].reshape([256, 256])
pylab.figure(3)
diff --git a/samples/python/s018_plugin.py b/samples/python/s018_plugin.py
index 31cca95..85b5486 100644
--- a/samples/python/s018_plugin.py
+++ b/samples/python/s018_plugin.py
@@ -30,30 +30,38 @@ import six
# Define the plugin class (has to subclass astra.plugin.base)
# Note that usually, these will be defined in a separate package/module
-class SIRTPlugin(astra.plugin.base):
- """Example of an ASTRA plugin class, implementing a simple 2D SIRT algorithm.
+class LandweberPlugin(astra.plugin.base):
+ """Example of an ASTRA plugin class, implementing a simple 2D Landweber algorithm.
Options:
- 'rel_factor': relaxation factor (optional)
+ 'Relaxation': relaxation factor (optional)
"""
# The astra_name variable defines the name to use to
# call the plugin from ASTRA
- astra_name = "SIRT-PLUGIN"
+ astra_name = "LANDWEBER-PLUGIN"
- def initialize(self,cfg, rel_factor = 1):
+ def initialize(self,cfg, Relaxation = 1):
self.W = astra.OpTomo(cfg['ProjectorId'])
self.vid = cfg['ReconstructionDataId']
self.sid = cfg['ProjectionDataId']
- self.rel = rel_factor
+ self.rel = Relaxation
def run(self, its):
v = astra.data2d.get_shared(self.vid)
s = astra.data2d.get_shared(self.sid)
+ tv = np.zeros(v.shape, dtype=np.float32)
+ ts = np.zeros(s.shape, dtype=np.float32)
W = self.W
for i in range(its):
- v[:] += self.rel*(W.T*(s - (W*v).reshape(s.shape))).reshape(v.shape)/s.size
+ W.FP(v,out=ts)
+ ts -= s # ts = W*v - s
+
+ W.BP(ts,out=tv)
+ tv *= self.rel / s.size
+
+ v -= tv # v = v - rel * W'*(W*v-s) / s.size
if __name__=='__main__':
@@ -75,20 +83,20 @@ if __name__=='__main__':
# First we import the package that contains the plugin
import s018_plugin
# Then, we register the plugin class with ASTRA
- astra.plugin.register(s018_plugin.SIRTPlugin)
+ astra.plugin.register(s018_plugin.LandweberPlugin)
# Get a list of registered plugins
six.print_(astra.plugin.get_registered())
# To get help on a registered plugin, use get_help
- six.print_(astra.plugin.get_help('SIRT-PLUGIN'))
+ six.print_(astra.plugin.get_help('LANDWEBER-PLUGIN'))
# Create data structures
sid = astra.data2d.create('-sino', proj_geom, sinogram)
vid = astra.data2d.create('-vol', vol_geom)
# Create config using plugin name
- cfg = astra.astra_dict('SIRT-PLUGIN')
+ cfg = astra.astra_dict('LANDWEBER-PLUGIN')
cfg['ProjectorId'] = proj_id
cfg['ProjectionDataId'] = sid
cfg['ReconstructionDataId'] = vid
@@ -103,18 +111,18 @@ if __name__=='__main__':
rec = astra.data2d.get(vid)
# Options for the plugin go in cfg['option']
- cfg = astra.astra_dict('SIRT-PLUGIN')
+ cfg = astra.astra_dict('LANDWEBER-PLUGIN')
cfg['ProjectorId'] = proj_id
cfg['ProjectionDataId'] = sid
cfg['ReconstructionDataId'] = vid
cfg['option'] = {}
- cfg['option']['rel_factor'] = 1.5
+ cfg['option']['Relaxation'] = 1.5
alg_id_rel = astra.algorithm.create(cfg)
astra.algorithm.run(alg_id_rel, 100)
rec_rel = astra.data2d.get(vid)
# We can also use OpTomo to call the plugin
- rec_op = W.reconstruct('SIRT-PLUGIN', sinogram, 100, extraOptions={'rel_factor':1.5})
+ rec_op = W.reconstruct('LANDWEBER-PLUGIN', sinogram, 100, extraOptions={'Relaxation':1.5})
import pylab as pl
pl.gray()
diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp
index 084ba8c..c63faaa 100644
--- a/src/CompositeGeometryManager.cpp
+++ b/src/CompositeGeometryManager.cpp
@@ -196,6 +196,69 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div
return true;
}
+
+static std::pair<double, double> reduceProjectionVertical(const CVolumeGeometry3D* pVolGeom, const CProjectionGeometry3D* pProjGeom)
+{
+ double vmin_g, vmax_g;
+
+ // reduce self to only cover intersection with projection of VolumePart
+ // (Project corners of volume, take bounding box)
+
+ for (int i = 0; i < pProjGeom->getProjectionCount(); ++i) {
+
+ double vol_u[8];
+ double vol_v[8];
+
+ double pixx = pVolGeom->getPixelLengthX();
+ double pixy = pVolGeom->getPixelLengthY();
+ double pixz = pVolGeom->getPixelLengthZ();
+
+ // TODO: Is 0.5 sufficient?
+ double xmin = pVolGeom->getWindowMinX() - 0.5 * pixx;
+ double xmax = pVolGeom->getWindowMaxX() + 0.5 * pixx;
+ double ymin = pVolGeom->getWindowMinY() - 0.5 * pixy;
+ double ymax = pVolGeom->getWindowMaxY() + 0.5 * pixy;
+ double zmin = pVolGeom->getWindowMinZ() - 0.5 * pixz;
+ double zmax = pVolGeom->getWindowMaxZ() + 0.5 * pixz;
+
+ pProjGeom->projectPoint(xmin, ymin, zmin, i, vol_u[0], vol_v[0]);
+ pProjGeom->projectPoint(xmin, ymin, zmax, i, vol_u[1], vol_v[1]);
+ pProjGeom->projectPoint(xmin, ymax, zmin, i, vol_u[2], vol_v[2]);
+ pProjGeom->projectPoint(xmin, ymax, zmax, i, vol_u[3], vol_v[3]);
+ pProjGeom->projectPoint(xmax, ymin, zmin, i, vol_u[4], vol_v[4]);
+ pProjGeom->projectPoint(xmax, ymin, zmax, i, vol_u[5], vol_v[5]);
+ pProjGeom->projectPoint(xmax, ymax, zmin, i, vol_u[6], vol_v[6]);
+ pProjGeom->projectPoint(xmax, ymax, zmax, i, vol_u[7], vol_v[7]);
+
+ double vmin = vol_v[0];
+ double vmax = vol_v[0];
+
+ for (int j = 1; j < 8; ++j) {
+ if (vol_v[j] < vmin)
+ vmin = vol_v[j];
+ if (vol_v[j] > vmax)
+ vmax = vol_v[j];
+ }
+
+ if (i == 0 || vmin < vmin_g)
+ vmin_g = vmin;
+ if (i == 0 || vmax > vmax_g)
+ vmax_g = vmax;
+ }
+
+ if (vmin_g < -1.0)
+ vmin_g = -1.0;
+ if (vmax_g > pProjGeom->getDetectorRowCount())
+ vmax_g = pProjGeom->getDetectorRowCount();
+
+ if (vmin_g >= vmax_g)
+ vmin_g = vmax_g = 0.0;
+
+ return std::pair<double, double>(vmin_g, vmax_g);
+
+}
+
+
CCompositeGeometryManager::CPart::CPart(const CPart& other)
{
eType = other.eType;
@@ -236,119 +299,135 @@ size_t CCompositeGeometryManager::CPart::getSize()
}
+static bool testVolumeRange(const std::pair<double, double>& fullRange,
+ const CVolumeGeometry3D *pVolGeom,
+ const CProjectionGeometry3D *pProjGeom,
+ int zmin, int zmax)
+{
+ double pixz = pVolGeom->getPixelLengthZ();
+
+ CVolumeGeometry3D test(pVolGeom->getGridColCount(),
+ pVolGeom->getGridRowCount(),
+ zmax - zmin,
+ pVolGeom->getWindowMinX(),
+ pVolGeom->getWindowMinY(),
+ pVolGeom->getWindowMinZ() + zmin * pixz,
+ pVolGeom->getWindowMaxX(),
+ pVolGeom->getWindowMaxY(),
+ pVolGeom->getWindowMinZ() + zmax * pixz);
+
+
+ std::pair<double, double> subRange = reduceProjectionVertical(&test, pProjGeom);
+
+ // empty
+ if (subRange.first == subRange.second)
+ return true;
+
+ // fully outside of fullRange
+ if (subRange.first >= fullRange.second || subRange.second <= fullRange.first)
+ return true;
+
+ return false;
+}
+
CCompositeGeometryManager::CPart* CCompositeGeometryManager::CVolumePart::reduce(const CPart *_other)
{
const CProjectionPart *other = dynamic_cast<const CProjectionPart *>(_other);
assert(other);
- // TODO: Is 0.5 sufficient?
- double umin = -0.5;
- double umax = other->pGeom->getDetectorColCount() + 0.5;
- double vmin = -0.5;
- double vmax = other->pGeom->getDetectorRowCount() + 0.5;
-
- double uu[4];
- double vv[4];
- uu[0] = umin; vv[0] = vmin;
- uu[1] = umin; vv[1] = vmax;
- uu[2] = umax; vv[2] = vmin;
- uu[3] = umax; vv[3] = vmax;
-
- double pixx = pGeom->getPixelLengthX();
- double pixy = pGeom->getPixelLengthY();
- double pixz = pGeom->getPixelLengthZ();
- double xmin = pGeom->getWindowMinX() - 0.5 * pixx;
- double xmax = pGeom->getWindowMaxX() + 0.5 * pixx;
- double ymin = pGeom->getWindowMinY() - 0.5 * pixy;
- double ymax = pGeom->getWindowMaxY() + 0.5 * pixy;
-
- // NB: Flipped
- double zmax = pGeom->getWindowMinZ() - 2.5 * pixz;
- double zmin = pGeom->getWindowMaxZ() + 2.5 * pixz;
-
- // TODO: This isn't as tight as it could be.
- // In particular it won't detect the detector being
- // missed entirely on the u side.
-
- for (int i = 0; i < other->pGeom->getProjectionCount(); ++i) {
- for (int j = 0; j < 4; ++j) {
- double px, py, pz;
-
- other->pGeom->backprojectPointX(i, uu[j], vv[j], xmin, py, pz);
- //ASTRA_DEBUG("%f %f (%f - %f)", py, pz, ymin, ymax);
- if (pz < zmin) zmin = pz;
- if (pz > zmax) zmax = pz;
- other->pGeom->backprojectPointX(i, uu[j], vv[j], xmax, py, pz);
- //ASTRA_DEBUG("%f %f (%f - %f)", py, pz, ymin, ymax);
- if (pz < zmin) zmin = pz;
- if (pz > zmax) zmax = pz;
-
- other->pGeom->backprojectPointY(i, uu[j], vv[j], ymin, px, pz);
- //ASTRA_DEBUG("%f %f (%f - %f)", px, pz, xmin, xmax);
- if (pz < zmin) zmin = pz;
- if (pz > zmax) zmax = pz;
- other->pGeom->backprojectPointY(i, uu[j], vv[j], ymax, px, pz);
- //ASTRA_DEBUG("%f %f (%f - %f)", px, pz, xmin, xmax);
- if (pz < zmin) zmin = pz;
- if (pz > zmax) zmax = pz;
+ std::pair<double, double> fullRange = reduceProjectionVertical(pGeom, other->pGeom);
+
+ int top_slice = 0, bottom_slice = 0;
+
+ if (fullRange.first < fullRange.second) {
+
+
+ // TOP SLICE
+
+ int zmin = 0;
+ int zmax = pGeom->getGridSliceCount()-1; // (Don't try empty region)
+
+ // Setting top slice to zmin is always valid.
+
+ while (zmin < zmax) {
+ int zmid = (zmin + zmax + 1) / 2;
+
+ bool ok = testVolumeRange(fullRange, pGeom, other->pGeom,
+ 0, zmid);
+
+ ASTRA_DEBUG("binsearch min: [%d,%d], %d, %s", zmin, zmax, zmid, ok ? "ok" : "removed too much");
+
+ if (ok)
+ zmin = zmid;
+ else
+ zmax = zmid - 1;
}
- }
- //ASTRA_DEBUG("coord extent: %f - %f", zmin, zmax);
+ top_slice = zmin;
- // Clip both zmin and zmax to get rid of extreme (or infinite) values
- // NB: When individual pz values are +/-Inf, the sign is determined
- // by ray direction and on which side of the face the ray passes.
- if (zmin < pGeom->getWindowMinZ() - 2*pixz)
- zmin = pGeom->getWindowMinZ() - 2*pixz;
- if (zmin > pGeom->getWindowMaxZ() + 2*pixz)
- zmin = pGeom->getWindowMaxZ() + 2*pixz;
- if (zmax < pGeom->getWindowMinZ() - 2*pixz)
- zmax = pGeom->getWindowMinZ() - 2*pixz;
- if (zmax > pGeom->getWindowMaxZ() + 2*pixz)
- zmax = pGeom->getWindowMaxZ() + 2*pixz;
- zmin = (zmin - pixz - pGeom->getWindowMinZ()) / pixz;
- zmax = (zmax + pixz - pGeom->getWindowMinZ()) / pixz;
+ // BOTTOM SLICE
- int _zmin = (int)floor(zmin);
- int _zmax = (int)ceil(zmax);
+ zmin = top_slice + 1; // (Don't try empty region)
+ zmax = pGeom->getGridSliceCount();
- //ASTRA_DEBUG("index extent: %d - %d", _zmin, _zmax);
+ // Setting bottom slice to zmax is always valid
+
+ while (zmin < zmax) {
+ int zmid = (zmin + zmax) / 2;
+
+ bool ok = testVolumeRange(fullRange, pGeom, other->pGeom,
+ zmid, pGeom->getGridSliceCount());
+
+ ASTRA_DEBUG("binsearch max: [%d,%d], %d, %s", zmin, zmax, zmid, ok ? "ok" : "removed too much");
+
+ if (ok)
+ zmax = zmid;
+ else
+ zmin = zmid + 1;
+
+ }
- if (_zmin < 0)
- _zmin = 0;
- if (_zmax > pGeom->getGridSliceCount())
- _zmax = pGeom->getGridSliceCount();
+ bottom_slice = zmax;
- if (_zmax <= _zmin) {
- _zmin = _zmax = 0;
}
- //ASTRA_DEBUG("adjusted extent: %d - %d", _zmin, _zmax);
+
+ ASTRA_DEBUG("found extent: %d - %d", top_slice, bottom_slice);
+
+ top_slice -= 1;
+ if (top_slice < 0)
+ top_slice = 0;
+ bottom_slice += 1;
+ if (bottom_slice >= pGeom->getGridSliceCount())
+ bottom_slice = pGeom->getGridSliceCount();
+
+ ASTRA_DEBUG("adjusted extent: %d - %d", top_slice, bottom_slice);
+
+ double pixz = pGeom->getPixelLengthZ();
CVolumePart *sub = new CVolumePart();
sub->subX = this->subX;
sub->subY = this->subY;
- sub->subZ = this->subZ + _zmin;
+ sub->subZ = this->subZ + top_slice;
sub->pData = pData;
- if (_zmin == _zmax) {
+ if (top_slice == bottom_slice) {
sub->pGeom = 0;
} else {
sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(),
pGeom->getGridRowCount(),
- _zmax - _zmin,
+ bottom_slice - top_slice,
pGeom->getWindowMinX(),
pGeom->getWindowMinY(),
- pGeom->getWindowMinZ() + _zmin * pixz,
+ pGeom->getWindowMinZ() + top_slice * pixz,
pGeom->getWindowMaxX(),
pGeom->getWindowMaxY(),
- pGeom->getWindowMinZ() + _zmax * pixz);
+ pGeom->getWindowMinZ() + bottom_slice * pixz);
}
- ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d", this->subZ, this->subZ + pGeom->getGridSliceCount(), this->subZ + _zmin, this->subZ + _zmax);
+ ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d ( %f - %f )", this->subZ, this->subZ + pGeom->getGridSliceCount(), this->subZ + top_slice, this->subZ + bottom_slice, pGeom->getWindowMinZ() + top_slice * pixz, pGeom->getWindowMinZ() + bottom_slice * pixz);
return sub;
}
@@ -583,7 +662,9 @@ void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::T
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -630,7 +711,9 @@ void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::T
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -677,7 +760,9 @@ void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::T
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
@@ -742,62 +827,16 @@ void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, s
}
+
CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::reduce(const CPart *_other)
{
const CVolumePart *other = dynamic_cast<const CVolumePart *>(_other);
assert(other);
- double vmin_g, vmax_g;
-
- // reduce self to only cover intersection with projection of VolumePart
- // (Project corners of volume, take bounding box)
-
- for (int i = 0; i < pGeom->getProjectionCount(); ++i) {
-
- double vol_u[8];
- double vol_v[8];
-
- double pixx = other->pGeom->getPixelLengthX();
- double pixy = other->pGeom->getPixelLengthY();
- double pixz = other->pGeom->getPixelLengthZ();
-
- // TODO: Is 0.5 sufficient?
- double xmin = other->pGeom->getWindowMinX() - 0.5 * pixx;
- double xmax = other->pGeom->getWindowMaxX() + 0.5 * pixx;
- double ymin = other->pGeom->getWindowMinY() - 0.5 * pixy;
- double ymax = other->pGeom->getWindowMaxY() + 0.5 * pixy;
- double zmin = other->pGeom->getWindowMinZ() - 0.5 * pixz;
- double zmax = other->pGeom->getWindowMaxZ() + 0.5 * pixz;
-
- pGeom->projectPoint(xmin, ymin, zmin, i, vol_u[0], vol_v[0]);
- pGeom->projectPoint(xmin, ymin, zmax, i, vol_u[1], vol_v[1]);
- pGeom->projectPoint(xmin, ymax, zmin, i, vol_u[2], vol_v[2]);
- pGeom->projectPoint(xmin, ymax, zmax, i, vol_u[3], vol_v[3]);
- pGeom->projectPoint(xmax, ymin, zmin, i, vol_u[4], vol_v[4]);
- pGeom->projectPoint(xmax, ymin, zmax, i, vol_u[5], vol_v[5]);
- pGeom->projectPoint(xmax, ymax, zmin, i, vol_u[6], vol_v[6]);
- pGeom->projectPoint(xmax, ymax, zmax, i, vol_u[7], vol_v[7]);
-
- double vmin = vol_v[0];
- double vmax = vol_v[0];
-
- for (int j = 1; j < 8; ++j) {
- if (vol_v[j] < vmin)
- vmin = vol_v[j];
- if (vol_v[j] > vmax)
- vmax = vol_v[j];
- }
-
- if (i == 0 || vmin < vmin_g)
- vmin_g = vmin;
- if (i == 0 || vmax > vmax_g)
- vmax_g = vmax;
- }
-
- // fprintf(stderr, "v extent: %f %f\n", vmin_g, vmax_g);
-
- int _vmin = (int)floor(vmin_g - 1.0f);
- int _vmax = (int)ceil(vmax_g + 1.0f);
+ std::pair<double, double> r = reduceProjectionVertical(other->pGeom, pGeom);
+ // fprintf(stderr, "v extent: %f %f\n", r.first, r.second);
+ int _vmin = (int)floor(r.first - 1.0);
+ int _vmax = (int)ceil(r.second + 1.0);
if (_vmin < 0)
_vmin = 0;
if (_vmax > pGeom->getDetectorRowCount())
@@ -837,7 +876,11 @@ void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManage
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
+
+ ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
for (int x = -(rem / 2); x < sliceCount; x += blockSize) {
int newsubX = x;
@@ -879,7 +922,11 @@ void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManage
size_t m = std::min(maxSize / sliceSize, maxDim);
size_t blockSize = computeLinearSplit(m, div, sliceCount);
- int rem = sliceCount % blockSize;
+ int rem = blockSize - (sliceCount % blockSize);
+ if (rem == blockSize)
+ rem = 0;
+
+ ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
for (int z = -(rem / 2); z < sliceCount; z += blockSize) {
int newsubZ = z;
diff --git a/src/Utilities.cpp b/src/Utilities.cpp
index c9740bf..8b0ca94 100644
--- a/src/Utilities.cpp
+++ b/src/Utilities.cpp
@@ -28,10 +28,6 @@ $Id$
#include "astra/Utilities.h"
-#include <boost/algorithm/string.hpp>
-#include <boost/algorithm/string/split.hpp>
-#include <boost/algorithm/string/classification.hpp>
-
#include <sstream>
#include <locale>
#include <iomanip>
@@ -84,18 +80,16 @@ std::vector<double> stringToDoubleVector(const std::string &s)
template<typename T>
std::vector<T> stringToVector(const std::string& s)
{
- // split
- std::vector<std::string> items;
- boost::split(items, s, boost::is_any_of(",;"));
-
- // init list
std::vector<T> out;
- out.resize(items.size());
+ size_t current = 0;
+ size_t next;
+ do {
+ next = s.find_first_of(",;", current);
+ std::string t = s.substr(current, next - current);
+ out.push_back(stringTo<T>(t));
+ current = next + 1;
+ } while (next != std::string::npos);
- // loop elements
- for (unsigned int i = 0; i < items.size(); i++) {
- out[i] = stringTo<T>(items[i]);
- }
return out;
}
@@ -120,6 +114,18 @@ std::string doubleToString(double f)
template<> std::string toString(float f) { return floatToString(f); }
template<> std::string toString(double f) { return doubleToString(f); }
+void splitString(std::vector<std::string> &items, const std::string& s,
+ const char *delim)
+{
+ items.clear();
+ size_t current = 0;
+ size_t next;
+ do {
+ next = s.find_first_of(",;", current);
+ items.push_back(s.substr(current, next - current));
+ current = next + 1;
+ } while (next != std::string::npos);
+}
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-08 18:03:19 +0000