Merge branch 'master' into FDK

4 files changed, 281 insertions, 282 deletions

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);
+}
 
 }