Merge branch 'master' into FDK

20 files changed, 450 insertions, 408 deletions

diff --git a/README.md b/README.md
index b34cf38..fb1f641 100644
--- a/README.md
+++ b/README.md
@@ -30,7 +30,7 @@ cd build/linux
 ./autogen.sh   # when building a git version
 ./configure --with-cuda=/usr/local/cuda \
             --with-matlab=/usr/local/MATLAB/R2012a \
-            --with-python
+            --with-python \
             --prefix=/usr/local/astra
 make
 make install
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 aab755a..3dd1c22 100644
--- a/cuda/2d/fft.cu
+++ b/cuda/2d/fft.cu
@@ -304,44 +304,42 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
 	float * pfFilt = new float[_iFFTFourierDetectorCount];
 	float * pfW = new float[_iFFTFourierDetectorCount];
 
-#if 1
+	// 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;
+			}
+		}
 
-	for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
-	{
-		float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount;
+		pfData[0] = 0.25f;
 
-		// filt = 2*( 0:(order/2) )./order;
-		pfFilt[iDetectorIndex] = 2.0f * fRelIndex;
-		//pfFilt[iDetectorIndex] = 1.0f;
+		cdft(2*_iFFTRealDetectorCount, -1, pfData, ip, w);
+		delete[] ip;
+		delete[] w;
 
-		// w = 2*pi*(0:size(filt,2)-1)/order
-		pfW[iDetectorIndex] = 3.1415f * 2.0f * fRelIndex;
-	}
-#else
-
-	float *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;
-		}
+		iFilterCacheSize = _iFFTRealDetectorCount;
 	}
 
-	pfData[0] = 0.25f;
-
-	cdft(2*_iFFTRealDetectorCount, -1, pfData, ip, w);
-
 	for(int iDetectorIndex = 0; iDetectorIndex < _iFFTFourierDetectorCount; iDetectorIndex++)
 	{
 		float fRelIndex = (float)iDetectorIndex / (float)_iFFTRealDetectorCount;
@@ -350,11 +348,6 @@ void genFilter(E_FBPFILTER _eFilter, float _fD, int _iProjectionCount,
 		pfW[iDetectorIndex] = M_PI * 2.0f * fRelIndex;
 	}
 
-	delete[] pfData;
-	delete[] ip;
-	delete[] w;
-#endif
-
 	switch(_eFilter)
 	{
 		case FILTER_RAMLAK:
@@ -906,5 +899,4 @@ void downloadDebugFilterReal(float * _pfHostSinogram, int _iProjectionCount,
 	free(pfHostFilter);
 }
 
-
 #endif
diff --git a/include/astra/AstraObjectManager.h b/include/astra/AstraObjectManager.h
index ad89c2a..9faecbe 100644
--- a/include/astra/AstraObjectManager.h
+++ b/include/astra/AstraObjectManager.h
@@ -60,7 +60,7 @@ public:
 };
 
 
-class CAstraIndexManager : public Singleton<CAstraIndexManager> {
+class _AstraExport CAstraIndexManager : public Singleton<CAstraIndexManager> {
 public:
 	CAstraIndexManager() : m_iLastIndex(0) { }
 
diff --git a/include/astra/Fourier.h b/include/astra/Fourier.h
index ff26f82..68f9f38 100644
--- a/include/astra/Fourier.h
+++ b/include/astra/Fourier.h
@@ -76,7 +76,7 @@ namespace astra {
             }
         .
 */
-void cdft(int n, int isgn, float32 *a, int *ip, float32 *w);
+_AstraExport void cdft(int n, int isgn, float32 *a, int *ip, float32 *w);
 
 }
 
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 c5b4d3b..7c4f8e6 100644
--- a/src/CompositeGeometryManager.cpp
+++ b/src/CompositeGeometryManager.cpp
@@ -197,6 +197,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;
@@ -237,119 +300,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;
 }
@@ -584,7 +663,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);
 
@@ -631,7 +712,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);
 
@@ -678,7 +761,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);
 
@@ -743,62 +828,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())
@@ -838,7 +877,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;
@@ -880,7 +923,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/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp
index 90efd52..70462f7 100644
--- a/src/FilteredBackProjectionAlgorithm.cpp
+++ b/src/FilteredBackProjectionAlgorithm.cpp
@@ -275,7 +275,7 @@ void CFilteredBackProjectionAlgorithm::performFiltering(CFloat32ProjectionData2D
 
 
 	float32* pf = new float32[2 * iAngleCount * zpDetector];
-	int *ip = new int[int(2+sqrt(zpDetector)+1)];
+	int *ip = new int[int(2+sqrt((float)zpDetector)+1)];
 	ip[0]=0;
 	float32 *w = new float32[zpDetector/2];
 
diff --git a/src/Fourier.cpp b/src/Fourier.cpp
index 5ca22e6..c33f7bd 100644
--- a/src/Fourier.cpp
+++ b/src/Fourier.cpp
@@ -27,6 +27,7 @@ $Id$
 */
 
 #include "astra/Fourier.h"
+#include <cmath>
 
 namespace astra {
 
@@ -320,11 +321,45 @@ Appendix :
 */
 
 
-void cdft(int n, int isgn, float32 *a, int *ip, float32 *w)
+static int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w);
+static void bitrv208(float32 *a);
+static void bitrv208neg(float32 *a);
+static void bitrv216(float32 *a);
+static void bitrv216neg(float32 *a);
+static void bitrv2conj(int n, int *ip, float32 *a);
+static void bitrv2(int n, int *ip, float32 *a);
+static void cftb040(float32 *a);
+static void cftb1st(int n, float32 *a, float32 *w);
+static void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w);
+static void cftf040(float32 *a);
+static void cftf081(float32 *a, float32 *w);
+static void cftf082(float32 *a, float32 *w);
+static void cftf161(float32 *a, float32 *w);
+static void cftf162(float32 *a, float32 *w);
+static void cftf1st(int n, float32 *a, float32 *w);
+static void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
+static void cftfx41(int n, float32 *a, int nw, float32 *w);
+static void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w);
+static void cftmdl1(int n, float32 *a, float32 *w);
+static void cftmdl2(int n, float32 *a, float32 *w);
+static void *cftrec1_th(void *p);
+static void *cftrec2_th(void *p);
+static void cftrec4(int n, float32 *a, int nw, float32 *w);
+static void cftx020(float32 *a);
+static void dctsub(int n, float32 *a, int nc, float32 *c);
+static void dstsub(int n, float32 *a, int nc, float32 *c);
+static void makect(int nc, int *ip, float32 *c);
+static void makeipt(int nw, int *ip);
+static void makewt(int nw, int *ip, float32 *w);
+static void rftbsub(int n, float32 *a, int nc, float32 *c);
+static void rftfsub(int n, float32 *a, int nc, float32 *c);
+#ifdef USE_CDFT_THREADS
+static void cftrec4_th(int n, float32 *a, int nw, float32 *w);
+#endif /* USE_CDFT_THREADS */
+    
+  
+_AstraExport void cdft(int n, int isgn, float32 *a, int *ip, float32 *w)
 {
-    void makewt(int nw, int *ip, float32 *w);
-    void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
-    void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w);
     int nw;
     
     nw = ip[0];
@@ -340,14 +375,8 @@ void cdft(int n, int isgn, float32 *a, int *ip, float32 *w)
 }
 
 
-void rdft(int n, int isgn, float32 *a, int *ip, float32 *w)
+_AstraExport void rdft(int n, int isgn, float32 *a, int *ip, float32 *w)
 {
-    void makewt(int nw, int *ip, float32 *w);
-    void makect(int nc, int *ip, float32 *c);
-    void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
-    void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w);
-    void rftfsub(int n, float32 *a, int nc, float32 *c);
-    void rftbsub(int n, float32 *a, int nc, float32 *c);
     int nw, nc;
     float32 xi;
     
@@ -384,15 +413,8 @@ void rdft(int n, int isgn, float32 *a, int *ip, float32 *w)
 }
 
 
-void ddct(int n, int isgn, float32 *a, int *ip, float32 *w)
+_AstraExport void ddct(int n, int isgn, float32 *a, int *ip, float32 *w)
 {
-    void makewt(int nw, int *ip, float32 *w);
-    void makect(int nc, int *ip, float32 *c);
-    void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
-    void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w);
-    void rftfsub(int n, float32 *a, int nc, float32 *c);
-    void rftbsub(int n, float32 *a, int nc, float32 *c);
-    void dctsub(int n, float32 *a, int nc, float32 *c);
     int j, nw, nc;
     float32 xr;
     
@@ -440,15 +462,8 @@ void ddct(int n, int isgn, float32 *a, int *ip, float32 *w)
 }
 
 
-void ddst(int n, int isgn, float32 *a, int *ip, float32 *w)
+_AstraExport void ddst(int n, int isgn, float32 *a, int *ip, float32 *w)
 {
-    void makewt(int nw, int *ip, float32 *w);
-    void makect(int nc, int *ip, float32 *c);
-    void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
-    void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w);
-    void rftfsub(int n, float32 *a, int nc, float32 *c);
-    void rftbsub(int n, float32 *a, int nc, float32 *c);
-    void dstsub(int n, float32 *a, int nc, float32 *c);
     int j, nw, nc;
     float32 xr;
     
@@ -496,13 +511,8 @@ void ddst(int n, int isgn, float32 *a, int *ip, float32 *w)
 }
 
 
-void dfct(int n, float32 *a, float32 *t, int *ip, float32 *w)
+_AstraExport void dfct(int n, float32 *a, float32 *t, int *ip, float32 *w)
 {
-    void makewt(int nw, int *ip, float32 *w);
-    void makect(int nc, int *ip, float32 *c);
-    void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
-    void rftfsub(int n, float32 *a, int nc, float32 *c);
-    void dctsub(int n, float32 *a, int nc, float32 *c);
     int j, k, l, m, mh, nw, nc;
     float32 xr, xi, yr, yi;
     
@@ -589,13 +599,8 @@ void dfct(int n, float32 *a, float32 *t, int *ip, float32 *w)
 }
 
 
-void dfst(int n, float32 *a, float32 *t, int *ip, float32 *w)
+_AstraExport void dfst(int n, float32 *a, float32 *t, int *ip, float32 *w)
 {
-    void makewt(int nw, int *ip, float32 *w);
-    void makect(int nc, int *ip, float32 *c);
-    void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w);
-    void rftfsub(int n, float32 *a, int nc, float32 *c);
-    void dstsub(int n, float32 *a, int nc, float32 *c);
     int j, k, l, m, mh, nw, nc;
     float32 xr, xi, yr, yi;
     
@@ -675,12 +680,8 @@ void dfst(int n, float32 *a, float32 *t, int *ip, float32 *w)
 
 /* -------- initializing routines -------- */
 
-
-#include <math.h>
-
-void makewt(int nw, int *ip, float32 *w)
+static void makewt(int nw, int *ip, float32 *w)
 {
-    void makeipt(int nw, int *ip);
     int j, nwh, nw0, nw1;
     float32 delta, wn4r, wk1r, wk1i, wk3r, wk3i;
     
@@ -739,7 +740,7 @@ void makewt(int nw, int *ip, float32 *w)
 }
 
 
-void makeipt(int nw, int *ip)
+static void makeipt(int nw, int *ip)
 {
     int j, l, m, m2, p, q;
     
@@ -759,7 +760,7 @@ void makeipt(int nw, int *ip)
 }
 
 
-void makect(int nc, int *ip, float32 *c)
+static void makect(int nc, int *ip, float32 *c)
 {
     int j, nch;
     float32 delta;
@@ -835,23 +836,8 @@ void makect(int nc, int *ip, float32 *c)
 #endif /* USE_CDFT_WINTHREADS */
 
 
-void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w)
+static void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w)
 {
-    void bitrv2(int n, int *ip, float32 *a);
-    void bitrv216(float32 *a);
-    void bitrv208(float32 *a);
-    void cftf1st(int n, float32 *a, float32 *w);
-    void cftrec4(int n, float32 *a, int nw, float32 *w);
-    void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w);
-    void cftfx41(int n, float32 *a, int nw, float32 *w);
-    void cftf161(float32 *a, float32 *w);
-    void cftf081(float32 *a, float32 *w);
-    void cftf040(float32 *a);
-    void cftx020(float32 *a);
-#ifdef USE_CDFT_THREADS
-    void cftrec4_th(int n, float32 *a, int nw, float32 *w);
-#endif /* USE_CDFT_THREADS */
-    
     if (n > 8) {
         if (n > 32) {
             cftf1st(n, a, &w[nw - (n >> 2)]);
@@ -883,23 +869,8 @@ void cftfsub(int n, float32 *a, int *ip, int nw, float32 *w)
 }
 
 
-void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w)
+static void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w)
 {
-    void bitrv2conj(int n, int *ip, float32 *a);
-    void bitrv216neg(float32 *a);
-    void bitrv208neg(float32 *a);
-    void cftb1st(int n, float32 *a, float32 *w);
-    void cftrec4(int n, float32 *a, int nw, float32 *w);
-    void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w);
-    void cftfx41(int n, float32 *a, int nw, float32 *w);
-    void cftf161(float32 *a, float32 *w);
-    void cftf081(float32 *a, float32 *w);
-    void cftb040(float32 *a);
-    void cftx020(float32 *a);
-#ifdef USE_CDFT_THREADS
-    void cftrec4_th(int n, float32 *a, int nw, float32 *w);
-#endif /* USE_CDFT_THREADS */
-    
     if (n > 8) {
         if (n > 32) {
             cftb1st(n, a, &w[nw - (n >> 2)]);
@@ -931,7 +902,7 @@ void cftbsub(int n, float32 *a, int *ip, int nw, float32 *w)
 }
 
 
-void bitrv2(int n, int *ip, float32 *a)
+static void bitrv2(int n, int *ip, float32 *a)
 {
     int j, j1, k, k1, l, m, nh, nm;
     float32 xr, xi, yr, yi;
@@ -1278,7 +1249,7 @@ void bitrv2(int n, int *ip, float32 *a)
 }
 
 
-void bitrv2conj(int n, int *ip, float32 *a)
+static void bitrv2conj(int n, int *ip, float32 *a)
 {
     int j, j1, k, k1, l, m, nh, nm;
     float32 xr, xi, yr, yi;
@@ -1633,7 +1604,7 @@ void bitrv2conj(int n, int *ip, float32 *a)
 }
 
 
-void bitrv216(float32 *a)
+static void bitrv216(float32 *a)
 {
     float32 x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
         x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i, 
@@ -1690,7 +1661,7 @@ void bitrv216(float32 *a)
 }
 
 
-void bitrv216neg(float32 *a)
+static void bitrv216neg(float32 *a)
 {
     float32 x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
         x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i, 
@@ -1760,7 +1731,7 @@ void bitrv216neg(float32 *a)
 }
 
 
-void bitrv208(float32 *a)
+static void bitrv208(float32 *a)
 {
     float32 x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i;
     
@@ -1783,7 +1754,7 @@ void bitrv208(float32 *a)
 }
 
 
-void bitrv208neg(float32 *a)
+static void bitrv208neg(float32 *a)
 {
     float32 x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i, 
         x5r, x5i, x6r, x6i, x7r, x7i;
@@ -1819,7 +1790,7 @@ void bitrv208neg(float32 *a)
 }
 
 
-void cftf1st(int n, float32 *a, float32 *w)
+static void cftf1st(int n, float32 *a, float32 *w)
 {
     int j, j0, j1, j2, j3, k, m, mh;
     float32 wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, 
@@ -2025,7 +1996,7 @@ void cftf1st(int n, float32 *a, float32 *w)
 }
 
 
-void cftb1st(int n, float32 *a, float32 *w)
+static void cftb1st(int n, float32 *a, float32 *w)
 {
     int j, j0, j1, j2, j3, k, m, mh;
     float32 wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i, 
@@ -2242,10 +2213,8 @@ struct cdft_arg_st {
 typedef struct cdft_arg_st cdft_arg_t;
 
 
-void cftrec4_th(int n, float32 *a, int nw, float32 *w)
+static void cftrec4_th(int n, float32 *a, int nw, float32 *w)
 {
-    void *cftrec1_th(void *p);
-    void *cftrec2_th(void *p);
     int i, idiv4, m, nthread;
     cdft_thread_t th[4];
     cdft_arg_t ag[4];
@@ -2276,11 +2245,8 @@ void cftrec4_th(int n, float32 *a, int nw, float32 *w)
 }
 
 
-void *cftrec1_th(void *p)
+static void *cftrec1_th(void *p)
 {
-    int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w);
-    void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w);
-    void cftmdl1(int n, float32 *a, float32 *w);
     int isplt, j, k, m, n, n0, nw;
     float32 *a, *w;
     
@@ -2305,11 +2271,8 @@ void *cftrec1_th(void *p)
 }
 
 
-void *cftrec2_th(void *p)
+static void *cftrec2_th(void *p)
 {
-    int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w);
-    void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w);
-    void cftmdl2(int n, float32 *a, float32 *w);
     int isplt, j, k, m, n, n0, nw;
     float32 *a, *w;
     
@@ -2337,11 +2300,8 @@ void *cftrec2_th(void *p)
 #endif /* USE_CDFT_THREADS */
 
 
-void cftrec4(int n, float32 *a, int nw, float32 *w)
+static void cftrec4(int n, float32 *a, int nw, float32 *w)
 {
-    int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w);
-    void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w);
-    void cftmdl1(int n, float32 *a, float32 *w);
     int isplt, j, k, m;
     
     m = n;
@@ -2361,8 +2321,6 @@ void cftrec4(int n, float32 *a, int nw, float32 *w)
 
 int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w)
 {
-    void cftmdl1(int n, float32 *a, float32 *w);
-    void cftmdl2(int n, float32 *a, float32 *w);
     int i, isplt, m;
     
     if ((k & 3) != 0) {
@@ -2394,15 +2352,8 @@ int cfttree(int n, int j, int k, float32 *a, int nw, float32 *w)
 }
 
 
-void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w)
+static void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w)
 {
-    void cftmdl1(int n, float32 *a, float32 *w);
-    void cftmdl2(int n, float32 *a, float32 *w);
-    void cftf161(float32 *a, float32 *w);
-    void cftf162(float32 *a, float32 *w);
-    void cftf081(float32 *a, float32 *w);
-    void cftf082(float32 *a, float32 *w);
-    
     if (n == 512) {
         cftmdl1(128, a, &w[nw - 64]);
         cftf161(a, &w[nw - 8]);
@@ -2459,7 +2410,7 @@ void cftleaf(int n, int isplt, float32 *a, int nw, float32 *w)
 }
 
 
-void cftmdl1(int n, float32 *a, float32 *w)
+static void cftmdl1(int n, float32 *a, float32 *w)
 {
     int j, j0, j1, j2, j3, k, m, mh;
     float32 wn4r, wk1r, wk1i, wk3r, wk3i;
@@ -2569,7 +2520,7 @@ void cftmdl1(int n, float32 *a, float32 *w)
 }
 
 
-void cftmdl2(int n, float32 *a, float32 *w)
+static void cftmdl2(int n, float32 *a, float32 *w)
 {
     int j, j0, j1, j2, j3, k, kr, m, mh;
     float32 wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i;
@@ -2703,13 +2654,8 @@ void cftmdl2(int n, float32 *a, float32 *w)
 }
 
 
-void cftfx41(int n, float32 *a, int nw, float32 *w)
+static void cftfx41(int n, float32 *a, int nw, float32 *w)
 {
-    void cftf161(float32 *a, float32 *w);
-    void cftf162(float32 *a, float32 *w);
-    void cftf081(float32 *a, float32 *w);
-    void cftf082(float32 *a, float32 *w);
-    
     if (n == 128) {
         cftf161(a, &w[nw - 8]);
         cftf162(&a[32], &w[nw - 32]);
@@ -2724,7 +2670,7 @@ void cftfx41(int n, float32 *a, int nw, float32 *w)
 }
 
 
-void cftf161(float32 *a, float32 *w)
+static void cftf161(float32 *a, float32 *w)
 {
     float32 wn4r, wk1r, wk1i, 
         x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
@@ -2883,7 +2829,7 @@ void cftf161(float32 *a, float32 *w)
 }
 
 
-void cftf162(float32 *a, float32 *w)
+static void cftf162(float32 *a, float32 *w)
 {
     float32 wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i, 
         x0r, x0i, x1r, x1i, x2r, x2i, 
@@ -3066,7 +3012,7 @@ void cftf162(float32 *a, float32 *w)
 }
 
 
-void cftf081(float32 *a, float32 *w)
+static void cftf081(float32 *a, float32 *w)
 {
     float32 wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i, 
         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
@@ -3128,7 +3074,7 @@ void cftf081(float32 *a, float32 *w)
 }
 
 
-void cftf082(float32 *a, float32 *w)
+static void cftf082(float32 *a, float32 *w)
 {
     float32 wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i, 
         y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i, 
@@ -3200,7 +3146,7 @@ void cftf082(float32 *a, float32 *w)
 }
 
 
-void cftf040(float32 *a)
+static void cftf040(float32 *a)
 {
     float32 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
     
@@ -3223,7 +3169,7 @@ void cftf040(float32 *a)
 }
 
 
-void cftb040(float32 *a)
+static void cftb040(float32 *a)
 {
     float32 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
     
@@ -3246,7 +3192,7 @@ void cftb040(float32 *a)
 }
 
 
-void cftx020(float32 *a)
+static void cftx020(float32 *a)
 {
     float32 x0r, x0i;
     
@@ -3259,7 +3205,7 @@ void cftx020(float32 *a)
 }
 
 
-void rftfsub(int n, float32 *a, int nc, float32 *c)
+static void rftfsub(int n, float32 *a, int nc, float32 *c)
 {
     int j, k, kk, ks, m;
     float32 wkr, wki, xr, xi, yr, yi;
@@ -3284,7 +3230,7 @@ void rftfsub(int n, float32 *a, int nc, float32 *c)
 }
 
 
-void rftbsub(int n, float32 *a, int nc, float32 *c)
+static void rftbsub(int n, float32 *a, int nc, float32 *c)
 {
     int j, k, kk, ks, m;
     float32 wkr, wki, xr, xi, yr, yi;
@@ -3309,7 +3255,7 @@ void rftbsub(int n, float32 *a, int nc, float32 *c)
 }
 
 
-void dctsub(int n, float32 *a, int nc, float32 *c)
+static void dctsub(int n, float32 *a, int nc, float32 *c)
 {
     int j, k, kk, ks, m;
     float32 wkr, wki, xr;
@@ -3330,7 +3276,7 @@ void dctsub(int n, float32 *a, int nc, float32 *c)
 }
 
 
-void dstsub(int n, float32 *a, int nc, float32 *c)
+static void dstsub(int n, float32 *a, int nc, float32 *c)
 {
     int j, k, kk, ks, m;
     float32 wkr, wki, xr;
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);
+}
 
 }