Merge branch 'master'

27 files changed, 2259 insertions, 567 deletions

diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp
new file mode 100644
index 0000000..9be4797
--- /dev/null
+++ b/src/CompositeGeometryManager.cpp
@@ -0,0 +1,993 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+           2014-2015, CWI, Amsterdam
+
+Contact: astra@uantwerpen.be
+Website: http://sf.net/projects/astra-toolbox
+
+This file is part of the ASTRA Toolbox.
+
+
+The ASTRA Toolbox is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+The ASTRA Toolbox is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+*/
+
+#include "astra/CompositeGeometryManager.h"
+
+#ifdef ASTRA_CUDA
+
+#include "astra/GeometryUtil3D.h"
+#include "astra/VolumeGeometry3D.h"
+#include "astra/ConeProjectionGeometry3D.h"
+#include "astra/ConeVecProjectionGeometry3D.h"
+#include "astra/ParallelProjectionGeometry3D.h"
+#include "astra/ParallelVecProjectionGeometry3D.h"
+#include "astra/Projector3D.h"
+#include "astra/CudaProjector3D.h"
+#include "astra/Float32ProjectionData3DMemory.h"
+#include "astra/Float32VolumeData3DMemory.h"
+#include "astra/Logging.h"
+
+#include "../cuda/3d/mem3d.h"
+
+#include <cstring>
+
+namespace astra {
+
+// JOB:
+//  
+// VolumePart
+// ProjectionPart
+// FP-or-BP
+// SET-or-ADD
+
+
+// Running a set of jobs:
+//
+// [ Assume OUTPUT Parts in a single JobSet don't alias?? ]
+// Group jobs by output Part
+// One thread per group?
+
+// Automatically split parts if too large
+// Performance model for odd-sized tasks?
+// Automatically split parts if not enough tasks to fill available GPUs
+
+
+// Splitting:
+// Constraints:
+//   number of sub-parts divisible by N
+//   max size of sub-parts
+
+// For splitting on both input and output side:
+//   How to divide up memory? (Optimization problem; compute/benchmark)
+//   (First approach: 0.5/0.5)
+
+
+
+bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div, TJobSet &split)
+{
+	split.clear();
+
+	for (TJobSet::const_iterator i = jobs.begin(); i != jobs.end(); ++i)
+	{
+		CPart* pOutput = i->first;
+		const TJobList &L = i->second;
+
+		// 1. Split output part
+		// 2. Per sub-part:
+		//    a. reduce input part
+		//    b. split input part
+		//    c. create jobs for new (input,output) subparts
+
+		TPartList splitOutput = pOutput->split(maxSize/3, div);
+
+		for (TJobList::const_iterator j = L.begin(); j != L.end(); ++j)
+		{
+			const SJob &job = *j;
+
+			for (TPartList::iterator i_out = splitOutput.begin();
+			     i_out != splitOutput.end(); ++i_out)
+			{
+				boost::shared_ptr<CPart> outputPart = *i_out;
+
+				SJob newjob;
+				newjob.pOutput = outputPart;
+				newjob.eType = j->eType;
+				newjob.eMode = j->eMode;
+				newjob.pProjector = j->pProjector;
+
+				CPart* input = job.pInput->reduce(outputPart.get());
+
+				if (input->getSize() == 0) {
+					ASTRA_DEBUG("Empty input");
+					newjob.eType = SJob::JOB_NOP;
+					split[outputPart.get()].push_back(newjob);
+					continue;
+				}
+
+				size_t remainingSize = ( maxSize - outputPart->getSize() ) / 2;
+
+				TPartList splitInput = input->split(remainingSize, 1);
+				delete input;
+				ASTRA_DEBUG("Input split into %d parts", splitInput.size());
+
+				for (TPartList::iterator i_in = splitInput.begin();
+				     i_in != splitInput.end(); ++i_in)
+				{
+					newjob.pInput = *i_in;
+
+					split[outputPart.get()].push_back(newjob);
+
+					// Second and later (input) parts should always be added to
+					// output of first (input) part.
+					newjob.eMode = SJob::MODE_ADD;
+				}
+
+			
+			}
+
+		}
+	}
+
+	return true;
+}
+
+CCompositeGeometryManager::CPart::CPart(const CPart& other)
+{
+	eType = other.eType;
+	pData = other.pData;
+	subX = other.subX;
+	subY = other.subY;
+	subZ = other.subZ;
+}
+
+CCompositeGeometryManager::CVolumePart::CVolumePart(const CVolumePart& other)
+ : CPart(other)
+{
+	pGeom = other.pGeom->clone();
+}
+
+CCompositeGeometryManager::CVolumePart::~CVolumePart()
+{
+	delete pGeom;
+}
+
+void CCompositeGeometryManager::CVolumePart::getDims(size_t &x, size_t &y, size_t &z)
+{
+	if (!pGeom) {
+		x = y = z = 0;
+		return;
+	}
+
+	x = pGeom->getGridColCount();
+	y = pGeom->getGridRowCount();
+	z = pGeom->getGridSliceCount();
+}
+
+size_t CCompositeGeometryManager::CPart::getSize()
+{
+	size_t x, y, z;
+	getDims(x, y, z);
+	return x * y * z;
+}
+
+
+
+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;
+		}
+	}
+
+	//ASTRA_DEBUG("coord extent: %f - %f", zmin, zmax);
+
+	zmin = (zmin - pixz - pGeom->getWindowMinZ()) / pixz;
+	zmax = (zmax + pixz - pGeom->getWindowMinZ()) / pixz;
+
+	int _zmin = (int)floor(zmin);
+	int _zmax = (int)ceil(zmax);
+
+	//ASTRA_DEBUG("index extent: %d - %d", _zmin, _zmax);
+
+	if (_zmin < 0)
+		_zmin = 0;
+	if (_zmax > pGeom->getGridSliceCount())
+		_zmax = pGeom->getGridSliceCount();
+
+	if (_zmax <= _zmin) {
+		_zmin = _zmax = 0;
+	}
+	//ASTRA_DEBUG("adjusted extent: %d - %d", _zmin, _zmax);
+
+	CVolumePart *sub = new CVolumePart();
+	sub->subX = this->subX;
+	sub->subY = this->subY;
+	sub->subZ = this->subZ + _zmin;
+	sub->pData = pData;
+
+	if (_zmin == _zmax) {
+		sub->pGeom = 0;
+	} else {
+		sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(),
+		                                   pGeom->getGridRowCount(),
+		                                   _zmax - _zmin,
+		                                   pGeom->getWindowMinX(),
+		                                   pGeom->getWindowMinY(),
+		                                   pGeom->getWindowMinZ() + _zmin * pixz,
+		                                   pGeom->getWindowMaxX(),
+		                                   pGeom->getWindowMaxY(),
+		                                   pGeom->getWindowMinZ() + _zmax * pixz);
+	}
+
+	ASTRA_DEBUG("Reduce volume from %d - %d to %d - %d", this->subZ, this->subZ +  pGeom->getGridSliceCount(), this->subZ + _zmin, this->subZ + _zmax);
+
+	return sub;
+}
+
+
+
+static size_t ceildiv(size_t a, size_t b) {
+    return (a + b - 1) / b;
+}
+
+static size_t computeVerticalSplit(size_t maxBlock, int div, size_t sliceCount)
+{
+    size_t blockSize = maxBlock;
+    size_t blockCount = ceildiv(sliceCount, blockSize);
+
+    // Increase number of blocks to be divisible by div
+    size_t divCount = div * ceildiv(blockCount, div);
+
+    // If divCount is above sqrt(number of slices), then
+    // we can't guarantee divisibility by div, but let's try anyway
+    if (ceildiv(sliceCount, ceildiv(sliceCount, divCount)) % div == 0) {
+        blockCount = divCount;
+    } else {
+        // If divisibility isn't achievable, we may want to optimize
+        // differently.
+        // TODO: Figure out how to model and optimize this.
+    }
+
+    // Final adjustment to make blocks more evenly sized
+    // (This can't make the blocks larger)
+    blockSize = ceildiv(sliceCount, blockCount); 
+
+    ASTRA_DEBUG("%ld %ld -> %ld * %ld", sliceCount, maxBlock, blockCount, blockSize);
+
+    assert(blockSize <= maxBlock);
+    assert((divCount * divCount > sliceCount) || (blockCount % div) == 0);
+
+    return blockSize;
+}
+
+template<class V, class P>
+static V* getProjectionVectors(const P* geom);
+
+template<>
+SConeProjection* getProjectionVectors(const CConeProjectionGeometry3D* pProjGeom)
+{
+	return genConeProjections(pProjGeom->getProjectionCount(),
+	                          pProjGeom->getDetectorColCount(),
+	                          pProjGeom->getDetectorRowCount(),
+	                          pProjGeom->getOriginSourceDistance(),
+	                          pProjGeom->getOriginDetectorDistance(),
+	                          pProjGeom->getDetectorSpacingX(),
+	                          pProjGeom->getDetectorSpacingY(),
+	                          pProjGeom->getProjectionAngles());
+}
+
+template<>
+SConeProjection* getProjectionVectors(const CConeVecProjectionGeometry3D* pProjGeom)
+{
+	int nth = pProjGeom->getProjectionCount();
+
+	SConeProjection* pProjs = new SConeProjection[nth];
+	for (int i = 0; i < nth; ++i)
+		pProjs[i] = pProjGeom->getProjectionVectors()[i];
+
+	return pProjs;
+}
+
+template<>
+SPar3DProjection* getProjectionVectors(const CParallelProjectionGeometry3D* pProjGeom)
+{
+	return genPar3DProjections(pProjGeom->getProjectionCount(),
+	                           pProjGeom->getDetectorColCount(),
+	                           pProjGeom->getDetectorRowCount(),
+	                           pProjGeom->getDetectorSpacingX(),
+	                           pProjGeom->getDetectorSpacingY(),
+	                           pProjGeom->getProjectionAngles());
+}
+
+template<>
+SPar3DProjection* getProjectionVectors(const CParallelVecProjectionGeometry3D* pProjGeom)
+{
+	int nth = pProjGeom->getProjectionCount();
+
+	SPar3DProjection* pProjs = new SPar3DProjection[nth];
+	for (int i = 0; i < nth; ++i)
+		pProjs[i] = pProjGeom->getProjectionVectors()[i];
+
+	return pProjs;
+}
+
+
+template<class V>
+static void translateProjectionVectors(V* pProjs, int count, double dv)
+{
+	for (int i = 0; i < count; ++i) {
+		pProjs[i].fDetSX += dv * pProjs[i].fDetVX;
+		pProjs[i].fDetSY += dv * pProjs[i].fDetVY;
+		pProjs[i].fDetSZ += dv * pProjs[i].fDetVZ;
+	}
+}
+
+
+
+static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry3D* pProjGeom, int v, int size)
+{
+	// First convert to vectors, then translate, then convert into new object
+
+	const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(pProjGeom);
+	const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(pProjGeom);
+	const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(pProjGeom);
+	const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(pProjGeom);
+
+	if (conegeom || conevec3dgeom) {
+		SConeProjection* pConeProjs;
+		if (conegeom) {
+			pConeProjs = getProjectionVectors<SConeProjection>(conegeom);
+		} else {
+			pConeProjs = getProjectionVectors<SConeProjection>(conevec3dgeom);
+		}
+
+		translateProjectionVectors(pConeProjs, pProjGeom->getProjectionCount(), v);
+
+		CProjectionGeometry3D* ret = new CConeVecProjectionGeometry3D(pProjGeom->getProjectionCount(),
+		                                                              size,
+		                                                              pProjGeom->getDetectorColCount(),
+		                                                              pConeProjs);
+
+
+		delete[] pConeProjs;
+		return ret;
+	} else {
+		assert(par3dgeom || parvec3dgeom);
+		SPar3DProjection* pParProjs;
+		if (par3dgeom) {
+			pParProjs = getProjectionVectors<SPar3DProjection>(par3dgeom);
+		} else {
+			pParProjs = getProjectionVectors<SPar3DProjection>(parvec3dgeom);
+		}
+
+		translateProjectionVectors(pParProjs, pProjGeom->getProjectionCount(), v);
+
+		CProjectionGeometry3D* ret = new CParallelVecProjectionGeometry3D(pProjGeom->getProjectionCount(),
+		                                                                  size,
+		                                                                  pProjGeom->getDetectorColCount(),
+		                                                                  pParProjs);
+
+		delete[] pParProjs;
+		return ret;
+	}
+
+}
+
+
+
+// split self into sub-parts:
+// - each no bigger than maxSize
+// - number of sub-parts is divisible by div
+// - maybe all approximately the same size?
+CCompositeGeometryManager::TPartList CCompositeGeometryManager::CVolumePart::split(size_t maxSize, int div)
+{
+	TPartList ret;
+
+	if (true) {
+		// Split in vertical direction only at first, until we figure out
+		// a model for splitting in other directions
+
+		size_t sliceSize = ((size_t) pGeom->getGridColCount()) * pGeom->getGridRowCount();
+		int sliceCount = pGeom->getGridSliceCount();
+		size_t blockSize = computeVerticalSplit(maxSize / sliceSize, div, sliceCount);
+
+		int rem = sliceCount % blockSize;
+
+		ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
+
+		for (int z = -(rem / 2); z < sliceCount; z += blockSize) {
+			int newsubZ = z;
+			if (newsubZ < 0) newsubZ = 0;
+			int endZ = z + blockSize;
+			if (endZ > sliceCount) endZ = sliceCount;
+			int size = endZ - newsubZ;
+
+			CVolumePart *sub = new CVolumePart();
+			sub->subX = this->subX;
+			sub->subY = this->subY;
+			sub->subZ = this->subZ + newsubZ;
+
+			ASTRA_DEBUG("VolumePart split %d %d %d -> %p", sub->subX, sub->subY, sub->subZ, (void*)sub);
+
+			double shift = pGeom->getPixelLengthZ() * newsubZ;
+
+			sub->pData = pData;
+			sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(),
+			                                   pGeom->getGridRowCount(),
+			                                   size,
+			                                   pGeom->getWindowMinX(),
+			                                   pGeom->getWindowMinY(),
+			                                   pGeom->getWindowMinZ() + shift,
+			                                   pGeom->getWindowMaxX(),
+			                                   pGeom->getWindowMaxY(),
+			                                   pGeom->getWindowMinZ() + shift + size * pGeom->getPixelLengthZ());
+
+			ret.push_back(boost::shared_ptr<CPart>(sub));
+		}
+	}
+
+	return ret;
+}
+
+CCompositeGeometryManager::CVolumePart* CCompositeGeometryManager::CVolumePart::clone() const
+{
+	return new CVolumePart(*this);
+}
+
+CCompositeGeometryManager::CProjectionPart::CProjectionPart(const CProjectionPart& other)
+ : CPart(other)
+{
+	pGeom = other.pGeom->clone();
+}
+
+CCompositeGeometryManager::CProjectionPart::~CProjectionPart()
+{
+	delete pGeom;
+}
+
+void CCompositeGeometryManager::CProjectionPart::getDims(size_t &x, size_t &y, size_t &z)
+{
+	if (!pGeom) {
+		x = y = z = 0;
+		return;
+	}
+
+	x = pGeom->getDetectorColCount();
+	y = pGeom->getProjectionCount();
+	z = pGeom->getDetectorRowCount();
+}
+
+
+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);
+	if (_vmin < 0)
+		_vmin = 0;
+	if (_vmax > pGeom->getDetectorRowCount())
+		_vmax = pGeom->getDetectorRowCount();
+
+	if (_vmin >= _vmax) {
+		_vmin = _vmax = 0;
+	}
+
+	CProjectionPart *sub = new CProjectionPart();
+	sub->subX = this->subX;
+	sub->subY = this->subY;
+	sub->subZ = this->subZ + _vmin;
+
+	sub->pData = pData;
+
+	if (_vmin == _vmax) {
+		sub->pGeom = 0;
+	} else {
+		sub->pGeom = getSubProjectionGeometry(pGeom, _vmin, _vmax - _vmin);
+	}
+
+	ASTRA_DEBUG("Reduce projection from %d - %d to %d - %d", this->subZ, this->subZ + pGeom->getDetectorRowCount(), this->subZ + _vmin, this->subZ + _vmax);
+
+	return sub;
+}
+
+
+CCompositeGeometryManager::TPartList CCompositeGeometryManager::CProjectionPart::split(size_t maxSize, int div)
+{
+	TPartList ret;
+
+	if (true) {
+		// Split in vertical direction only at first, until we figure out
+		// a model for splitting in other directions
+
+		size_t sliceSize = ((size_t) pGeom->getDetectorColCount()) * pGeom->getProjectionCount();
+		int sliceCount = pGeom->getDetectorRowCount();
+		size_t blockSize = computeVerticalSplit(maxSize / sliceSize, div, sliceCount);
+
+		int rem = sliceCount % blockSize;
+
+		for (int z = -(rem / 2); z < sliceCount; z += blockSize) {
+			int newsubZ = z;
+			if (newsubZ < 0) newsubZ = 0;
+			int endZ = z + blockSize;
+			if (endZ > sliceCount) endZ = sliceCount;
+			int size = endZ - newsubZ;
+
+			CProjectionPart *sub = new CProjectionPart();
+			sub->subX = this->subX;
+			sub->subY = this->subY;
+			sub->subZ = this->subZ + newsubZ;
+
+			ASTRA_DEBUG("ProjectionPart split %d %d %d -> %p", sub->subX, sub->subY, sub->subZ, (void*)sub);
+
+			sub->pData = pData;
+
+			sub->pGeom = getSubProjectionGeometry(pGeom, newsubZ, size);
+
+			ret.push_back(boost::shared_ptr<CPart>(sub));
+		}
+	}
+
+	return ret;
+
+}
+
+CCompositeGeometryManager::CProjectionPart* CCompositeGeometryManager::CProjectionPart::clone() const
+{
+	return new CProjectionPart(*this);
+}
+
+
+bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData,
+                                     CFloat32ProjectionData3DMemory *pProjData)
+{
+	ASTRA_DEBUG("CCompositeGeometryManager::doFP");
+	// Create single job for FP
+	// Run result
+
+	CVolumePart *input = new CVolumePart();
+	input->pData = pVolData;
+	input->subX = 0;
+	input->subY = 0;
+	input->subZ = 0;
+	input->pGeom = pVolData->getGeometry()->clone();
+	ASTRA_DEBUG("Main FP VolumePart -> %p", (void*)input);
+
+	CProjectionPart *output = new CProjectionPart();
+	output->pData = pProjData;
+	output->subX = 0;
+	output->subY = 0;
+	output->subZ = 0;
+	output->pGeom = pProjData->getGeometry()->clone();
+	ASTRA_DEBUG("Main FP ProjectionPart -> %p", (void*)output);
+
+	SJob FP;
+	FP.pInput = boost::shared_ptr<CPart>(input);
+	FP.pOutput = boost::shared_ptr<CPart>(output);
+	FP.pProjector = pProjector;
+	FP.eType = SJob::JOB_FP;
+	FP.eMode = SJob::MODE_SET;
+
+	TJobList L;
+	L.push_back(FP);
+
+	return doJobs(L);
+}
+
+bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData,
+                                     CFloat32ProjectionData3DMemory *pProjData)
+{
+	ASTRA_DEBUG("CCompositeGeometryManager::doBP");
+	// Create single job for BP
+	// Run result
+
+	CProjectionPart *input = new CProjectionPart();
+	input->pData = pProjData;
+	input->subX = 0;
+	input->subY = 0;
+	input->subZ = 0;
+	input->pGeom = pProjData->getGeometry()->clone();
+
+	CVolumePart *output = new CVolumePart();
+	output->pData = pVolData;
+	output->subX = 0;
+	output->subY = 0;
+	output->subZ = 0;
+	output->pGeom = pVolData->getGeometry()->clone();
+
+	SJob BP;
+	BP.pInput = boost::shared_ptr<CPart>(input);
+	BP.pOutput = boost::shared_ptr<CPart>(output);
+	BP.pProjector = pProjector;
+	BP.eType = SJob::JOB_BP;
+	BP.eMode = SJob::MODE_SET;
+
+	TJobList L;
+	L.push_back(BP);
+
+	return doJobs(L);
+}
+
+bool CCompositeGeometryManager::doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData)
+{
+	ASTRA_DEBUG("CCompositeGeometryManager::doFP, multi-volume");
+
+	std::vector<CFloat32VolumeData3DMemory *>::const_iterator i;
+	std::vector<boost::shared_ptr<CPart> > inputs;
+
+	for (i = volData.begin(); i != volData.end(); ++i) {
+		CVolumePart *input = new CVolumePart();
+		input->pData = *i;
+		input->subX = 0;
+		input->subY = 0;
+		input->subZ = 0;
+		input->pGeom = (*i)->getGeometry()->clone();
+
+		inputs.push_back(boost::shared_ptr<CPart>(input));
+	}
+
+	std::vector<CFloat32ProjectionData3DMemory *>::const_iterator j;
+	std::vector<boost::shared_ptr<CPart> > outputs;
+
+	for (j = projData.begin(); j != projData.end(); ++j) {
+		CProjectionPart *output = new CProjectionPart();
+		output->pData = *j;
+		output->subX = 0;
+		output->subY = 0;
+		output->subZ = 0;
+		output->pGeom = (*j)->getGeometry()->clone();
+
+		outputs.push_back(boost::shared_ptr<CPart>(output));
+	}
+
+	std::vector<boost::shared_ptr<CPart> >::iterator i2;
+	std::vector<boost::shared_ptr<CPart> >::iterator j2;
+	TJobList L;
+
+	for (i2 = outputs.begin(); i2 != outputs.end(); ++i2) {
+		SJob FP;
+		FP.eMode = SJob::MODE_SET;
+		for (j2 = inputs.begin(); j2 != inputs.end(); ++j2) {
+			FP.pInput = *j2;
+			FP.pOutput = *i2;
+			FP.pProjector = pProjector;
+			FP.eType = SJob::JOB_FP;
+			L.push_back(FP);
+
+			// Set first, add rest
+			FP.eMode = SJob::MODE_ADD;
+		}
+	}
+
+	return doJobs(L);
+}
+
+bool CCompositeGeometryManager::doBP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData)
+{
+	ASTRA_DEBUG("CCompositeGeometryManager::doBP, multi-volume");
+
+
+	std::vector<CFloat32VolumeData3DMemory *>::const_iterator i;
+	std::vector<boost::shared_ptr<CPart> > outputs;
+
+	for (i = volData.begin(); i != volData.end(); ++i) {
+		CVolumePart *output = new CVolumePart();
+		output->pData = *i;
+		output->subX = 0;
+		output->subY = 0;
+		output->subZ = 0;
+		output->pGeom = (*i)->getGeometry()->clone();
+
+		outputs.push_back(boost::shared_ptr<CPart>(output));
+	}
+
+	std::vector<CFloat32ProjectionData3DMemory *>::const_iterator j;
+	std::vector<boost::shared_ptr<CPart> > inputs;
+
+	for (j = projData.begin(); j != projData.end(); ++j) {
+		CProjectionPart *input = new CProjectionPart();
+		input->pData = *j;
+		input->subX = 0;
+		input->subY = 0;
+		input->subZ = 0;
+		input->pGeom = (*j)->getGeometry()->clone();
+
+		inputs.push_back(boost::shared_ptr<CPart>(input));
+	}
+
+	std::vector<boost::shared_ptr<CPart> >::iterator i2;
+	std::vector<boost::shared_ptr<CPart> >::iterator j2;
+	TJobList L;
+
+	for (i2 = outputs.begin(); i2 != outputs.end(); ++i2) {
+		SJob BP;
+		BP.eMode = SJob::MODE_SET;
+		for (j2 = inputs.begin(); j2 != inputs.end(); ++j2) {
+			BP.pInput = *j2;
+			BP.pOutput = *i2;
+			BP.pProjector = pProjector;
+			BP.eType = SJob::JOB_BP;
+			L.push_back(BP);
+
+			// Set first, add rest
+			BP.eMode = SJob::MODE_ADD;
+		}
+	}
+
+	return doJobs(L);
+}
+
+
+
+
+bool CCompositeGeometryManager::doJobs(TJobList &jobs)
+{
+	ASTRA_DEBUG("CCompositeGeometryManager::doJobs");
+
+	// Sort job list into job set by output part
+	TJobSet jobset;
+
+	for (TJobList::iterator i = jobs.begin(); i != jobs.end(); ++i) {
+		jobset[i->pOutput.get()].push_back(*i);
+	}
+
+	size_t maxSize = astraCUDA3d::availableGPUMemory();
+	if (maxSize == 0) {
+		ASTRA_WARN("Unable to get available GPU memory. Defaulting to 1GB.");
+		maxSize = 1024 * 1024 * 1024;
+	} else {
+		ASTRA_DEBUG("Detected %lu bytes of GPU memory", maxSize);
+	}
+	maxSize = (maxSize * 9) / 10;
+
+	maxSize /= sizeof(float);
+	int div = 1;
+
+	// TODO: Multi-GPU support
+
+	// Split jobs to fit
+	TJobSet split;
+	splitJobs(jobset, maxSize, div, split);
+	jobset.clear();
+
+	// Run jobs
+	
+	for (TJobSet::iterator iter = split.begin(); iter != split.end(); ++iter) {
+
+		CPart* output = iter->first;
+		TJobList& L = iter->second;
+
+		assert(!L.empty());
+
+		bool zero = L.begin()->eMode == SJob::MODE_SET;
+
+		size_t outx, outy, outz;
+		output->getDims(outx, outy, outz);
+
+		if (L.begin()->eType == SJob::JOB_NOP) {
+			// just zero output?
+			if (zero) {
+				for (size_t z = 0; z < outz; ++z) {
+					for (size_t y = 0; y < outy; ++y) {
+						float* ptr = output->pData->getData();
+						ptr += (z + output->subX) * (size_t)output->pData->getHeight() * (size_t)output->pData->getWidth();
+						ptr += (y + output->subY) * (size_t)output->pData->getWidth();
+						ptr += output->subX;
+						memset(ptr, 0, sizeof(float) * outx);
+					}
+				}
+			}
+			continue;
+		}
+
+
+		astraCUDA3d::SSubDimensions3D dstdims;
+		dstdims.nx = output->pData->getWidth();
+		dstdims.pitch = dstdims.nx;
+		dstdims.ny = output->pData->getHeight();
+		dstdims.nz = output->pData->getDepth();
+		dstdims.subnx = outx;
+		dstdims.subny = outy;
+		dstdims.subnz = outz;
+		ASTRA_DEBUG("dstdims: %d,%d,%d in %d,%d,%d", dstdims.subnx, dstdims.subny, dstdims.subnz, dstdims.nx, dstdims.ny, dstdims.nz);
+		dstdims.subx = output->subX;
+		dstdims.suby = output->subY;
+		dstdims.subz = output->subZ;
+		float *dst = output->pData->getData();
+
+		astraCUDA3d::MemHandle3D outputMem = astraCUDA3d::allocateGPUMemory(outx, outy, outz, zero ? astraCUDA3d::INIT_ZERO : astraCUDA3d::INIT_NO);
+		bool ok = outputMem;
+
+		for (TJobList::iterator i = L.begin(); i != L.end(); ++i) {
+			SJob &j = *i;
+
+			assert(j.pInput);
+
+			CCudaProjector3D *projector = dynamic_cast<CCudaProjector3D*>(j.pProjector);
+			Cuda3DProjectionKernel projKernel = ker3d_default;
+			int detectorSuperSampling = 1;
+			int voxelSuperSampling = 1;
+			if (projector) {
+				projKernel = projector->getProjectionKernel();
+				detectorSuperSampling = projector->getDetectorSuperSampling();
+				voxelSuperSampling = projector->getVoxelSuperSampling();
+			}
+
+			size_t inx, iny, inz;
+			j.pInput->getDims(inx, iny, inz);
+			astraCUDA3d::MemHandle3D inputMem = astraCUDA3d::allocateGPUMemory(inx, iny, inz, astraCUDA3d::INIT_NO);
+
+			astraCUDA3d::SSubDimensions3D srcdims;
+			srcdims.nx = j.pInput->pData->getWidth();
+			srcdims.pitch = srcdims.nx;
+			srcdims.ny = j.pInput->pData->getHeight();
+			srcdims.nz = j.pInput->pData->getDepth();
+			srcdims.subnx = inx;
+			srcdims.subny = iny;
+			srcdims.subnz = inz;
+			srcdims.subx = j.pInput->subX;
+			srcdims.suby = j.pInput->subY;
+			srcdims.subz = j.pInput->subZ;
+			const float *src = j.pInput->pData->getDataConst();
+
+			ok = astraCUDA3d::copyToGPUMemory(src, inputMem, srcdims);
+			if (!ok) ASTRA_ERROR("Error copying input data to GPU");
+
+			if (j.eType == SJob::JOB_FP) {
+				assert(dynamic_cast<CVolumePart*>(j.pInput.get()));
+				assert(dynamic_cast<CProjectionPart*>(j.pOutput.get()));
+
+				ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing FP");
+
+				ok = astraCUDA3d::FP(((CProjectionPart*)j.pOutput.get())->pGeom, outputMem, ((CVolumePart*)j.pInput.get())->pGeom, inputMem, detectorSuperSampling, projKernel);
+				if (!ok) ASTRA_ERROR("Error performing sub-FP");
+				ASTRA_DEBUG("CCompositeGeometryManager::doJobs: FP done");
+			} else if (j.eType == SJob::JOB_BP) {
+				assert(dynamic_cast<CVolumePart*>(j.pOutput.get()));
+				assert(dynamic_cast<CProjectionPart*>(j.pInput.get()));
+
+				ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing BP");
+
+				ok = astraCUDA3d::BP(((CProjectionPart*)j.pInput.get())->pGeom, inputMem, ((CVolumePart*)j.pOutput.get())->pGeom, outputMem, voxelSuperSampling);
+				if (!ok) ASTRA_ERROR("Error performing sub-BP");
+				ASTRA_DEBUG("CCompositeGeometryManager::doJobs: BP done");
+			} else {
+				assert(false);
+			}
+
+			ok = astraCUDA3d::freeGPUMemory(inputMem);
+			if (!ok) ASTRA_ERROR("Error freeing GPU memory");
+
+		}
+
+		ok = astraCUDA3d::copyFromGPUMemory(dst, outputMem, dstdims);
+		if (!ok) ASTRA_ERROR("Error copying output data from GPU");
+		
+		ok = astraCUDA3d::freeGPUMemory(outputMem);
+		if (!ok) ASTRA_ERROR("Error freeing GPU memory");
+	}
+
+	return true;
+}
+
+
+
+}
+
+#endif
diff --git a/src/ConeProjectionGeometry3D.cpp b/src/ConeProjectionGeometry3D.cpp
index dd22eba..18f0f8a 100644
--- a/src/ConeProjectionGeometry3D.cpp
+++ b/src/ConeProjectionGeometry3D.cpp
@@ -29,6 +29,7 @@ $Id$
 #include "astra/ConeProjectionGeometry3D.h"
 
 #include "astra/Logging.h"
+#include "astra/GeometryUtil3D.h"
 
 #include <boost/lexical_cast.hpp>
 #include <cstring>
@@ -230,14 +231,14 @@ CVector3D CConeProjectionGeometry3D::getProjectionDirection(int _iProjectionInde
 	return ret;
 }
 
-void CConeProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ,
-                                                 int iAngleIndex,
-                                                 float32 &fU, float32 &fV) const
+void CConeProjectionGeometry3D::projectPoint(double fX, double fY, double fZ,
+                                             int iAngleIndex,
+                                             double &fU, double &fV) const
 {
 	ASTRA_ASSERT(iAngleIndex >= 0);
 	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
 
-	float alpha = m_pfProjectionAngles[iAngleIndex];
+	double alpha = m_pfProjectionAngles[iAngleIndex];
 
 	// Project point onto optical axis
 
@@ -245,14 +246,14 @@ void CConeProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ,
 	// Vector source->origin is (-sin(alpha), cos(alpha))
 
 	// Distance from source, projected on optical axis
-	float fD = -sin(alpha) * fX + cos(alpha) * fY + m_fOriginSourceDistance;
+	double fD = -sin(alpha) * fX + cos(alpha) * fY + m_fOriginSourceDistance;
 
 	// Scale fZ to detector plane
 	fV = detectorOffsetYToRowIndexFloat( (fZ * (m_fOriginSourceDistance + m_fOriginDetectorDistance)) / fD );
 
 
 	// Orthogonal distance in XY-plane to optical axis
-	float fS = cos(alpha) * fX + sin(alpha) * fY;
+	double fS = cos(alpha) * fX + sin(alpha) * fY;
 
 	// Scale fS to detector plane
 	fU = detectorOffsetXToColIndexFloat( (fS * (m_fOriginSourceDistance + m_fOriginDetectorDistance)) / fD );
@@ -261,5 +262,84 @@ void CConeProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ,
 
 }
 
+void CConeProjectionGeometry3D::backprojectPointX(int iAngleIndex, double fU, double fV,
+	                               double fX, double &fY, double &fZ) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	SConeProjection *projs = genConeProjections(1, m_iDetectorColCount, m_iDetectorRowCount,
+	                                           m_fOriginSourceDistance,
+	                                           m_fOriginDetectorDistance,
+	                                           m_fDetectorSpacingX, m_fDetectorSpacingY,
+	                                           &m_pfProjectionAngles[iAngleIndex]);
+
+	SConeProjection &proj = projs[0];
+
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fX - proj.fSrcX) / (px - proj.fSrcX);
+
+	fY = proj.fSrcY + a * (py - proj.fSrcY);
+	fZ = proj.fSrcZ + a * (pz - proj.fSrcZ);
+
+	delete[] projs;
+}
+
+void CConeProjectionGeometry3D::backprojectPointY(int iAngleIndex, double fU, double fV,
+	                               double fY, double &fX, double &fZ) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	SConeProjection *projs = genConeProjections(1, m_iDetectorColCount, m_iDetectorRowCount,
+	                                           m_fOriginSourceDistance,
+	                                           m_fOriginDetectorDistance,
+	                                           m_fDetectorSpacingX, m_fDetectorSpacingY,
+	                                           &m_pfProjectionAngles[iAngleIndex]);
+
+	SConeProjection &proj = projs[0];
+
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fY - proj.fSrcY) / (py - proj.fSrcY);
+
+	fX = proj.fSrcX + a * (px - proj.fSrcX);
+	fZ = proj.fSrcZ + a * (pz - proj.fSrcZ);
+
+	delete[] projs;
+}
+
+void CConeProjectionGeometry3D::backprojectPointZ(int iAngleIndex, double fU, double fV,
+	                               double fZ, double &fX, double &fY) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	SConeProjection *projs = genConeProjections(1, m_iDetectorColCount, m_iDetectorRowCount,
+	                                           m_fOriginSourceDistance,
+	                                           m_fOriginDetectorDistance,
+	                                           m_fDetectorSpacingX, m_fDetectorSpacingY,
+	                                           &m_pfProjectionAngles[iAngleIndex]);
+
+	SConeProjection &proj = projs[0];
+
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fZ - proj.fSrcZ) / (pz - proj.fSrcZ);
+
+	fX = proj.fSrcX + a * (px - proj.fSrcX);
+	fY = proj.fSrcY + a * (py - proj.fSrcY);
+
+	delete[] projs;
+}
+
+
 
 } // end namespace astra
diff --git a/src/ConeVecProjectionGeometry3D.cpp b/src/ConeVecProjectionGeometry3D.cpp
index 47ed630..86e3bd6 100644
--- a/src/ConeVecProjectionGeometry3D.cpp
+++ b/src/ConeVecProjectionGeometry3D.cpp
@@ -241,9 +241,9 @@ CVector3D CConeVecProjectionGeometry3D::getProjectionDirection(int _iProjectionI
 	return CVector3D(p.fDetSX + (u+0.5)*p.fDetUX + (v+0.5)*p.fDetVX - p.fSrcX, p.fDetSY + (u+0.5)*p.fDetUY + (v+0.5)*p.fDetVY - p.fSrcY, p.fDetSZ + (u+0.5)*p.fDetUZ + (v+0.5)*p.fDetVZ - p.fSrcZ);
 }
 
-void CConeVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ,
+void CConeVecProjectionGeometry3D::projectPoint(double fX, double fY, double fZ,
                                                  int iAngleIndex,
-                                                 float32 &fU, float32 &fV) const
+                                                 double &fU, double &fV) const
 {
 	ASTRA_ASSERT(iAngleIndex >= 0);
 	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
@@ -262,6 +262,60 @@ void CConeVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32
 }
 
 
+void CConeVecProjectionGeometry3D::backprojectPointX(int iAngleIndex, double fU, double fV,
+	                               double fX, double &fY, double &fZ) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	SConeProjection &proj = m_pProjectionAngles[iAngleIndex];
+
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fX - proj.fSrcX) / (px - proj.fSrcX);
+
+	fY = proj.fSrcY + a * (py - proj.fSrcY);
+	fZ = proj.fSrcZ + a * (pz - proj.fSrcZ);
+}
+
+void CConeVecProjectionGeometry3D::backprojectPointY(int iAngleIndex, double fU, double fV,
+	                               double fY, double &fX, double &fZ) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	SConeProjection &proj = m_pProjectionAngles[iAngleIndex];
+
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fY - proj.fSrcY) / (py - proj.fSrcY);
+
+	fX = proj.fSrcX + a * (px - proj.fSrcX);
+	fZ = proj.fSrcZ + a * (pz - proj.fSrcZ);
+}
+
+void CConeVecProjectionGeometry3D::backprojectPointZ(int iAngleIndex, double fU, double fV,
+	                               double fZ, double &fX, double &fY) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	SConeProjection &proj = m_pProjectionAngles[iAngleIndex];
+
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fZ - proj.fSrcZ) / (pz - proj.fSrcZ);
+
+	fX = proj.fSrcX + a * (px - proj.fSrcX);
+	fY = proj.fSrcY + a * (py - proj.fSrcY);
+}
+
 //----------------------------------------------------------------------------------------
 
 bool CConeVecProjectionGeometry3D::_check()
diff --git a/src/CudaBackProjectionAlgorithm.cpp b/src/CudaBackProjectionAlgorithm.cpp
index 365e058..a73f895 100644
--- a/src/CudaBackProjectionAlgorithm.cpp
+++ b/src/CudaBackProjectionAlgorithm.cpp
@@ -76,10 +76,9 @@ bool CCudaBackProjectionAlgorithm::initialize(const Config& _cfg)
 // Initialize - C++
 bool CCudaBackProjectionAlgorithm::initialize(CProjector2D* _pProjector,
                                      CFloat32ProjectionData2D* _pSinogram, 
-                                     CFloat32VolumeData2D* _pReconstruction,
-                                     int _iGPUindex, int _iPixelSuperSampling)
+                                     CFloat32VolumeData2D* _pReconstruction)
 {
-	m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, 1, _iPixelSuperSampling);
+	m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction);
 
 	if (!m_bIsInitialized)
 		return false;
diff --git a/src/CudaBackProjectionAlgorithm3D.cpp b/src/CudaBackProjectionAlgorithm3D.cpp
index e8e0433..ce8e111 100644
--- a/src/CudaBackProjectionAlgorithm3D.cpp
+++ b/src/CudaBackProjectionAlgorithm3D.cpp
@@ -37,6 +37,9 @@ $Id$
 #include "astra/ParallelProjectionGeometry3D.h"
 #include "astra/ParallelVecProjectionGeometry3D.h"
 #include "astra/ConeVecProjectionGeometry3D.h"
+#include "astra/CompositeGeometryManager.h"
+
+#include "astra/Logging.h"
 
 #include "../cuda/3d/astra3d.h"
 
@@ -87,6 +90,24 @@ bool CCudaBackProjectionAlgorithm3D::_check()
 }
 
 //---------------------------------------------------------------------------------------
+void CCudaBackProjectionAlgorithm3D::initializeFromProjector()
+{
+	m_iVoxelSuperSampling = 1;
+	m_iGPUIndex = -1;
+
+	CCudaProjector3D* pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector);
+	if (!pCudaProjector) {
+		if (m_pProjector) {
+			ASTRA_WARN("non-CUDA Projector3D passed to BP3D_CUDA");
+		}
+	} else {
+		m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling();
+		m_iGPUIndex = pCudaProjector->getGPUIndex();
+	}
+
+}
+
+//---------------------------------------------------------------------------------------
 // Initialize - Config
 bool CCudaBackProjectionAlgorithm3D::initialize(const Config& _cfg)
 {
@@ -103,32 +124,21 @@ bool CCudaBackProjectionAlgorithm3D::initialize(const Config& _cfg)
 		return false;
 	}
 
-	CCudaProjector3D* pCudaProjector = 0;
-	pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector);
-	if (!pCudaProjector) {
-		// TODO: Report
-	}
+	initializeFromProjector();
 
-	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
-	CC.markOptionParsed("GPUindex");
-
-
-	m_iVoxelSuperSampling = 1;
-	if (pCudaProjector)
-		m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling();
+	// Deprecated options
 	m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
 	CC.markOptionParsed("VoxelSuperSampling");
+	CC.markOptionParsed("GPUIndex");
+	if (!_cfg.self.hasOption("GPUIndex"))
+		CC.markOptionParsed("GPUindex");
 
-	CFloat32ProjectionData3DMemory* pSinoMem = dynamic_cast<CFloat32ProjectionData3DMemory*>(m_pSinogram);
-	ASTRA_ASSERT(pSinoMem);
-	const CProjectionGeometry3D* projgeom = pSinoMem->getGeometry();
-const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(projgeom);
-	const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom);
-	if (parvec3dgeom || par3dgeom) {
-		// This option is only supported for Par3D currently
-		m_bSIRTWeighting = _cfg.self.getOptionBool("SIRTWeighting", false);
-		CC.markOptionParsed("SIRTWeighting");
-	}
+
+
+	m_bSIRTWeighting = _cfg.self.getOptionBool("SIRTWeighting", false);
+	CC.markOptionParsed("SIRTWeighting");
 
 	// success
 	m_bIsInitialized = _check();
@@ -151,6 +161,8 @@ bool CCudaBackProjectionAlgorithm3D::initialize(CProjector3D* _pProjector,
 	m_pSinogram = _pSinogram;
 	m_pReconstruction = _pReconstruction;
 
+	initializeFromProjector();
+
 	// success
 	m_bIsInitialized = _check();
 	return m_bIsInitialized;
@@ -184,88 +196,24 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations)
 	ASTRA_ASSERT(pReconMem);
 
 	const CProjectionGeometry3D* projgeom = pSinoMem->getGeometry();
-	const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(projgeom);
-	const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(projgeom);
-	const CConeVecProjectionGeometry3D* conevecgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(projgeom);
-	const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom);
 	const CVolumeGeometry3D& volgeom = *pReconMem->getGeometry();
 
-	if (conegeom) {
-		astraCudaConeBP(pReconMem->getData(), pSinoMem->getDataConst(),
-		                volgeom.getGridColCount(),
-		                volgeom.getGridRowCount(),
-		                volgeom.getGridSliceCount(),
-		                conegeom->getProjectionCount(),
-		                conegeom->getDetectorColCount(),
-		                conegeom->getDetectorRowCount(),
-		                conegeom->getOriginSourceDistance(),
-		                conegeom->getOriginDetectorDistance(),
-		                conegeom->getDetectorSpacingX(),
-		                conegeom->getDetectorSpacingY(),
-		                conegeom->getProjectionAngles(),
-		                m_iGPUIndex, m_iVoxelSuperSampling);
-	} else if (par3dgeom) {
-		if (!m_bSIRTWeighting) {
-			astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(),
-			                 volgeom.getGridColCount(),
-			                 volgeom.getGridRowCount(),
-			                 volgeom.getGridSliceCount(),
-			                 par3dgeom->getProjectionCount(),
-			                 par3dgeom->getDetectorColCount(),
-			                 par3dgeom->getDetectorRowCount(),
-			                 par3dgeom->getDetectorSpacingX(),
-			                 par3dgeom->getDetectorSpacingY(),
-			                 par3dgeom->getProjectionAngles(),
-			                 m_iGPUIndex, m_iVoxelSuperSampling);
-		} else {
-			astraCudaPar3DBP_SIRTWeighted(pReconMem->getData(),
-			                 pSinoMem->getDataConst(),
-			                 volgeom.getGridColCount(),
-			                 volgeom.getGridRowCount(),
-			                 volgeom.getGridSliceCount(),
-			                 par3dgeom->getProjectionCount(),
-			                 par3dgeom->getDetectorColCount(),
-			                 par3dgeom->getDetectorRowCount(),
-			                 par3dgeom->getDetectorSpacingX(),
-			                 par3dgeom->getDetectorSpacingY(),
-			                 par3dgeom->getProjectionAngles(),
-			                 m_iGPUIndex, m_iVoxelSuperSampling);
-		}
-	} else if (parvec3dgeom) {
-		if (!m_bSIRTWeighting) {
-			astraCudaPar3DBP(pReconMem->getData(), pSinoMem->getDataConst(),
-			                 volgeom.getGridColCount(),
-			                 volgeom.getGridRowCount(),
-			                 volgeom.getGridSliceCount(),
-			                 parvec3dgeom->getProjectionCount(),
-			                 parvec3dgeom->getDetectorColCount(),
-			                 parvec3dgeom->getDetectorRowCount(),
-			                 parvec3dgeom->getProjectionVectors(),
-			                 m_iGPUIndex, m_iVoxelSuperSampling);
-		} else {
-			astraCudaPar3DBP_SIRTWeighted(pReconMem->getData(),
-			                 pSinoMem->getDataConst(),
-			                 volgeom.getGridColCount(),
-			                 volgeom.getGridRowCount(),
-			                 volgeom.getGridSliceCount(),
-			                 parvec3dgeom->getProjectionCount(),
-			                 parvec3dgeom->getDetectorColCount(),
-			                 parvec3dgeom->getDetectorRowCount(),
-			                 parvec3dgeom->getProjectionVectors(),
-			                 m_iGPUIndex, m_iVoxelSuperSampling);
-		}
-	} else if (conevecgeom) {
-		astraCudaConeBP(pReconMem->getData(), pSinoMem->getDataConst(),
-		                volgeom.getGridColCount(),
-		                volgeom.getGridRowCount(),
-		                volgeom.getGridSliceCount(),
-		                conevecgeom->getProjectionCount(),
-		                conevecgeom->getDetectorColCount(),
-		                conevecgeom->getDetectorRowCount(),
-		                conevecgeom->getProjectionVectors(),
-		                m_iGPUIndex, m_iVoxelSuperSampling);
+	if (m_bSIRTWeighting) {
+		astraCudaBP_SIRTWeighted(pReconMem->getData(),
+		                         pSinoMem->getDataConst(),
+		                         &volgeom, projgeom,
+		                         m_iGPUIndex, m_iVoxelSuperSampling);
 	} else {
-		ASTRA_ASSERT(false);
+
+#if 1
+		CCompositeGeometryManager cgm;
+
+		cgm.doBP(m_pProjector, pReconMem, pSinoMem);
+#else
+		astraCudaBP(pReconMem->getData(), pSinoMem->getDataConst(),
+		            &volgeom, projgeom,
+		            m_iGPUIndex, m_iVoxelSuperSampling);
+#endif
 	}
 
 }
diff --git a/src/CudaCglsAlgorithm.cpp b/src/CudaCglsAlgorithm.cpp
index 0cedff6..9dd4f78 100644
--- a/src/CudaCglsAlgorithm.cpp
+++ b/src/CudaCglsAlgorithm.cpp
@@ -77,11 +77,9 @@ bool CCudaCglsAlgorithm::initialize(const Config& _cfg)
 // Initialize - C++
 bool CCudaCglsAlgorithm::initialize(CProjector2D* _pProjector,
                                     CFloat32ProjectionData2D* _pSinogram, 
-                                    CFloat32VolumeData2D* _pReconstruction,
-                                    int _iGPUindex, int _iDetectorSuperSampling,
-                                    int _iPixelSuperSampling)
+                                    CFloat32VolumeData2D* _pReconstruction)
 {
-	m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, _iDetectorSuperSampling, _iPixelSuperSampling);
+	m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction);
 
 	if (!m_bIsInitialized)
 		return false;
diff --git a/src/CudaCglsAlgorithm3D.cpp b/src/CudaCglsAlgorithm3D.cpp
index f527dc5..abc18d1 100644
--- a/src/CudaCglsAlgorithm3D.cpp
+++ b/src/CudaCglsAlgorithm3D.cpp
@@ -37,6 +37,8 @@ $Id$
 #include "astra/ParallelVecProjectionGeometry3D.h"
 #include "astra/ConeVecProjectionGeometry3D.h"
 
+#include "astra/Logging.h"
+
 #include "../cuda/3d/astra3d.h"
 
 using namespace std;
@@ -90,6 +92,26 @@ bool CCudaCglsAlgorithm3D::_check()
 }
 
 //---------------------------------------------------------------------------------------
+void CCudaCglsAlgorithm3D::initializeFromProjector()
+{
+	m_iVoxelSuperSampling = 1;
+	m_iDetectorSuperSampling = 1;
+	m_iGPUIndex = -1;
+
+	CCudaProjector3D* pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector);
+	if (!pCudaProjector) {
+		if (m_pProjector) {
+			ASTRA_WARN("non-CUDA Projector3D passed to CGLS3D_CUDA");
+		}
+	} else {
+		m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling();
+		m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling();
+		m_iGPUIndex = pCudaProjector->getGPUIndex();
+	}
+
+}
+
+//---------------------------------------------------------------------------------------
 // Initialize - Config
 bool CCudaCglsAlgorithm3D::initialize(const Config& _cfg)
 {
@@ -107,27 +129,20 @@ bool CCudaCglsAlgorithm3D::initialize(const Config& _cfg)
 		return false;
 	}
 
-	CCudaProjector3D* pCudaProjector = 0;
-	pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector);
-	if (!pCudaProjector) {
-		// TODO: Report
-	}
+	initializeFromProjector();
 
-	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
-	CC.markOptionParsed("GPUindex");
-
-	m_iVoxelSuperSampling = 1;
-	m_iDetectorSuperSampling = 1;
-	if (pCudaProjector) {
-		// New interface
-		m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling();
-		m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling();
-	}
 	// Deprecated options
 	m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling);
 	m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
 	CC.markOptionParsed("VoxelSuperSampling");
 	CC.markOptionParsed("DetectorSuperSampling");
+	CC.markOptionParsed("GPUIndex");
+	if (!_cfg.self.hasOption("GPUIndex"))
+		CC.markOptionParsed("GPUindex");
+
+
 
 	m_pCgls = new AstraCGLS3d();
 
@@ -155,6 +170,8 @@ bool CCudaCglsAlgorithm3D::initialize(CProjector3D* _pProjector,
 	m_pSinogram = _pSinogram;
 	m_pReconstruction = _pReconstruction;
 
+	initializeFromProjector();
+
 	m_pCgls = new AstraCGLS3d;
 
 	m_bAstraCGLSInit = false;
@@ -187,9 +204,6 @@ void CCudaCglsAlgorithm3D::run(int _iNrIterations)
 	ASTRA_ASSERT(m_bIsInitialized);
 
 	const CProjectionGeometry3D* projgeom = m_pSinogram->getGeometry();
-	const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(projgeom);
-	const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom);
-	const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(projgeom);
 	const CVolumeGeometry3D& volgeom = *m_pReconstruction->getGeometry();
 
 	bool ok = true;
@@ -198,41 +212,7 @@ void CCudaCglsAlgorithm3D::run(int _iNrIterations)
 
 		ok &= m_pCgls->setGPUIndex(m_iGPUIndex);
 
-		ok &= m_pCgls->setReconstructionGeometry(volgeom.getGridColCount(),
-		                                         volgeom.getGridRowCount(),
-		                                         volgeom.getGridSliceCount());
-/*
-                                  unsigned int iProjAngles,
-                                  unsigned int iProjU,
-                                  unsigned int iProjV,
-                                  float fOriginSourceDistance,
-                                  float fOriginDetectorDistance,
-                                  float fDetUSize,
-                                  float fDetVSize,
-                                  const float *pfAngles)
-*/
-		if (conegeom) {
-			ok &= m_pCgls->setConeGeometry(conegeom->getProjectionCount(),
-			                               conegeom->getDetectorColCount(),
-			                               conegeom->getDetectorRowCount(),
-			                               conegeom->getOriginSourceDistance(),
-			                               conegeom->getOriginDetectorDistance(),
-			                               conegeom->getDetectorSpacingX(),
-			                               conegeom->getDetectorSpacingY(),
-			                               conegeom->getProjectionAngles());
-		} else if (parvec3dgeom) {
-			ok &= m_pCgls->setPar3DGeometry(parvec3dgeom->getProjectionCount(),
-			                                parvec3dgeom->getDetectorColCount(),
-			                                parvec3dgeom->getDetectorRowCount(),
-			                                parvec3dgeom->getProjectionVectors());
-		} else if (conevec3dgeom) {
-			ok &= m_pCgls->setConeGeometry(conevec3dgeom->getProjectionCount(),
-			                               conevec3dgeom->getDetectorColCount(),
-			                               conevec3dgeom->getDetectorRowCount(),
-			                               conevec3dgeom->getProjectionVectors());
-		} else {
-			ASTRA_ASSERT(false);
-		}
+		ok &= m_pCgls->setGeometry(&volgeom, projgeom);
 
 		ok &= m_pCgls->enableSuperSampling(m_iVoxelSuperSampling, m_iDetectorSuperSampling);
 
diff --git a/src/CudaEMAlgorithm.cpp b/src/CudaEMAlgorithm.cpp
index 5c71f3d..d0afd80 100644
--- a/src/CudaEMAlgorithm.cpp
+++ b/src/CudaEMAlgorithm.cpp
@@ -76,11 +76,9 @@ bool CCudaEMAlgorithm::initialize(const Config& _cfg)
 // Initialize - C++
 bool CCudaEMAlgorithm::initialize(CProjector2D* _pProjector,
                                      CFloat32ProjectionData2D* _pSinogram, 
-                                     CFloat32VolumeData2D* _pReconstruction,
-                                     int _iGPUindex, int _iDetectorSuperSampling,
-                                     int _iPixelSuperSampling)
+                                     CFloat32VolumeData2D* _pReconstruction)
 {
-	m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, _iDetectorSuperSampling, _iPixelSuperSampling);
+	m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction);
 
 	if (!m_bIsInitialized)
 		return false;
diff --git a/src/CudaFDKAlgorithm3D.cpp b/src/CudaFDKAlgorithm3D.cpp
index 667d926..1316daa 100644
--- a/src/CudaFDKAlgorithm3D.cpp
+++ b/src/CudaFDKAlgorithm3D.cpp
@@ -35,6 +35,8 @@ $Id$
 #include "astra/CudaProjector3D.h"
 #include "astra/ConeProjectionGeometry3D.h"
 
+#include "astra/Logging.h"
+
 #include "../cuda/3d/astra3d.h"
 
 using namespace std;
@@ -85,6 +87,24 @@ bool CCudaFDKAlgorithm3D::_check()
 }
 
 //---------------------------------------------------------------------------------------
+void CCudaFDKAlgorithm3D::initializeFromProjector()
+{
+	m_iVoxelSuperSampling = 1;
+	m_iGPUIndex = -1;
+
+	CCudaProjector3D* pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector);
+	if (!pCudaProjector) {
+		if (m_pProjector) {
+			ASTRA_WARN("non-CUDA Projector3D passed to FDK_CUDA");
+		}
+	} else {
+		m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling();
+		m_iGPUIndex = pCudaProjector->getGPUIndex();
+	}
+
+}
+
+//---------------------------------------------------------------------------------------
 // Initialize - Config
 bool CCudaFDKAlgorithm3D::initialize(const Config& _cfg)
 {
@@ -101,20 +121,18 @@ bool CCudaFDKAlgorithm3D::initialize(const Config& _cfg)
 		return false;
 	}
 
-	CCudaProjector3D* pCudaProjector = 0;
-	pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector);
-	if (!pCudaProjector) {
-		// TODO: Report
-	}
-
-	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
-	CC.markOptionParsed("GPUindex");
+	initializeFromProjector();
 
-	m_iVoxelSuperSampling = 1;
-	if (pCudaProjector)
-		m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling();
+	// Deprecated options
 	m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
 	CC.markOptionParsed("VoxelSuperSampling");
+	CC.markOptionParsed("GPUIndex");
+	if (!_cfg.self.hasOption("GPUIndex"))
+		CC.markOptionParsed("GPUindex");
+
+
 
 	m_bShortScan = _cfg.self.getOptionBool("ShortScan", false);
 	CC.markOptionParsed("ShortScan");
@@ -182,17 +200,7 @@ void CCudaFDKAlgorithm3D::run(int _iNrIterations)
 	bool ok = true;
 
 	ok = astraCudaFDK(pReconMem->getData(), pSinoMem->getDataConst(),
-	                  volgeom.getGridColCount(),
-	                  volgeom.getGridRowCount(),
-	                  volgeom.getGridSliceCount(),
-	                  conegeom->getProjectionCount(),
-	                  conegeom->getDetectorColCount(),
-	                  conegeom->getDetectorRowCount(),
-	                  conegeom->getOriginSourceDistance(),
-	                  conegeom->getOriginDetectorDistance(),
-	                  conegeom->getDetectorSpacingX(),
-	                  conegeom->getDetectorSpacingY(),
-	                  conegeom->getProjectionAngles(),
+	                  &volgeom, conegeom,
 	                  m_bShortScan, m_iGPUIndex, m_iVoxelSuperSampling);
 
 	ASTRA_ASSERT(ok);
diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp
index aac96d6..bcd70c4 100644
--- a/src/CudaFilteredBackProjectionAlgorithm.cpp
+++ b/src/CudaFilteredBackProjectionAlgorithm.cpp
@@ -67,6 +67,24 @@ CCudaFilteredBackProjectionAlgorithm::~CCudaFilteredBackProjectionAlgorithm()
 	}
 }
 
+void CCudaFilteredBackProjectionAlgorithm::initializeFromProjector()
+{
+	m_iPixelSuperSampling = 1;
+	m_iGPUIndex = -1;
+
+	// Projector
+	CCudaProjector2D* pCudaProjector = dynamic_cast<CCudaProjector2D*>(m_pProjector);
+	if (!pCudaProjector) {
+		if (m_pProjector) {
+			ASTRA_WARN("non-CUDA Projector2D passed to FBP_CUDA");
+		}
+	} else {
+		m_iPixelSuperSampling = pCudaProjector->getVoxelSuperSampling();
+		m_iGPUIndex = pCudaProjector->getGPUIndex();
+	}
+
+}
+
 bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
 {
 	ASTRA_ASSERT(_cfg.self);
@@ -163,35 +181,30 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
 	}
 	CC.markNodeParsed("FilterD"); // TODO: Only for some types!
 
-	// GPU number
-	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
-	CC.markOptionParsed("GPUindex");
-
-	m_iPixelSuperSampling = 1;
-	if (pCudaProjector) {
-		// New interface
-		m_iPixelSuperSampling = pCudaProjector->getVoxelSuperSampling();
-	}
-	// Deprecated options
-	m_iPixelSuperSampling = (int)_cfg.self.getOptionNumerical("PixelSuperSampling", m_iPixelSuperSampling);
-	CC.markOptionParsed("PixelSuperSampling");
-
-
 	// Fan beam short scan mode
 	if (m_pSinogram && dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pSinogram->getGeometry())) {
 		m_bShortScan = (int)_cfg.self.getOptionBool("ShortScan", false);
 		CC.markOptionParsed("ShortScan");
 	}
 
+	initializeFromProjector();
+
+	// Deprecated options
+	m_iPixelSuperSampling = (int)_cfg.self.getOptionNumerical("PixelSuperSampling", m_iPixelSuperSampling);
+	CC.markOptionParsed("PixelSuperSampling");
 
+	// GPU number
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
+	CC.markOptionParsed("GPUIndex");
+	if (!_cfg.self.hasOption("GPUIndex"))
+		CC.markOptionParsed("GPUindex");
 
 
 	m_pFBP = new AstraFBP;
 	m_bAstraFBPInit = false;
 
-	// success
-	m_bIsInitialized = true;
-	return m_bIsInitialized;
+	return check();
 }
 
 bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D * _pSinogram, CFloat32VolumeData2D * _pReconstruction, E_FBPFILTER _eFilter, const float * _pfFilter /* = NULL */, int _iFilterWidth /* = 0 */, int _iGPUIndex /* = 0 */, float _fFilterParameter /* = -1.0f */)
@@ -241,7 +254,7 @@ bool CCudaFilteredBackProjectionAlgorithm::initialize(CFloat32ProjectionData2D *
 
 	m_fFilterParameter = _fFilterParameter;
 
-	return m_bIsInitialized;
+	return check();
 }
 
 void CCudaFilteredBackProjectionAlgorithm::run(int _iNrIterations /* = 0 */)
@@ -361,7 +374,7 @@ bool CCudaFilteredBackProjectionAlgorithm::check()
 	ASTRA_CONFIG_CHECK(m_pReconstruction->isInitialized(), "FBP_CUDA", "Reconstruction Data Object Not Initialized.");
 
 	// check gpu index
-	ASTRA_CONFIG_CHECK(m_iGPUIndex >= -1, "FBP_CUDA", "GPUIndex must be a non-negative integer.");
+	ASTRA_CONFIG_CHECK(m_iGPUIndex >= -1, "FBP_CUDA", "GPUIndex must be a non-negative integer or -1.");
 	// check pixel supersampling
 	ASTRA_CONFIG_CHECK(m_iPixelSuperSampling >= 0, "FBP_CUDA", "PixelSuperSampling must be a non-negative integer.");
 
diff --git a/src/CudaForwardProjectionAlgorithm.cpp b/src/CudaForwardProjectionAlgorithm.cpp
index b382f2e..9ca13ae 100644
--- a/src/CudaForwardProjectionAlgorithm.cpp
+++ b/src/CudaForwardProjectionAlgorithm.cpp
@@ -38,8 +38,11 @@ $Id$
 #include <boost/lexical_cast.hpp>
 
 #include "astra/AstraObjectManager.h"
+#include "astra/ParallelProjectionGeometry2D.h"
 #include "astra/FanFlatProjectionGeometry2D.h"
 #include "astra/FanFlatVecProjectionGeometry2D.h"
+#include "astra/Float32ProjectionData2D.h"
+#include "astra/Float32VolumeData2D.h"
 #include "astra/CudaProjector2D.h"
 
 #include "astra/Logging.h"
@@ -66,6 +69,24 @@ CCudaForwardProjectionAlgorithm::~CCudaForwardProjectionAlgorithm()
 }
 
 //---------------------------------------------------------------------------------------
+void CCudaForwardProjectionAlgorithm::initializeFromProjector()
+{
+	m_iDetectorSuperSampling = 1;
+	m_iGPUIndex = -1;
+
+	// Projector
+	CCudaProjector2D* pCudaProjector = dynamic_cast<CCudaProjector2D*>(m_pProjector);
+	if (!pCudaProjector) {
+		if (m_pProjector) {
+			ASTRA_WARN("non-CUDA Projector2D passed to FP_CUDA");
+		}
+	} else {
+		m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling();
+		m_iGPUIndex = pCudaProjector->getGPUIndex();
+	}
+}
+
+//---------------------------------------------------------------------------------------
 // Initialize - Config
 bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg)
 {
@@ -74,14 +95,9 @@ bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg)
 
 	// Projector
 	XMLNode node = _cfg.self.getSingleNode("ProjectorId");
-	CCudaProjector2D* pCudaProjector = 0;
 	if (node) {
 		int id = boost::lexical_cast<int>(node.getContent());
-		CProjector2D *projector = CProjector2DManager::getSingleton().get(id);
-		pCudaProjector = dynamic_cast<CCudaProjector2D*>(projector);
-		if (!pCudaProjector) {
-			ASTRA_WARN("non-CUDA Projector2D passed to FP_CUDA");
-		}
+		m_pProjector = CProjector2DManager::getSingleton().get(id);
 	}
 	CC.markNodeParsed("ProjectorId");
 
@@ -101,22 +117,18 @@ bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg)
 	m_pVolume = dynamic_cast<CFloat32VolumeData2D*>(CData2DManager::getSingleton().get(id));
 	CC.markNodeParsed("VolumeDataId");
 
+	initializeFromProjector();
+
+	// Deprecated options
+	m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling);
+	CC.markOptionParsed("DetectorSuperSampling");
 	// GPU number
 	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
 	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
-	CC.markOptionParsed("GPUindex");
-	if (!_cfg.self.hasOption("GPUindex"))
-		CC.markOptionParsed("GPUIndex");
+	CC.markOptionParsed("GPUIndex");
+	if (!_cfg.self.hasOption("GPUIndex"))
+		CC.markOptionParsed("GPUindex");
 
-	// Detector supersampling factor
-	m_iDetectorSuperSampling = 1;
-	if (pCudaProjector) {
-		// New interface
-		m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling();
-	}
-	// Deprecated option
-	m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling);
-	CC.markOptionParsed("DetectorSuperSampling");
 
 
 	// return success
@@ -125,20 +137,16 @@ bool CCudaForwardProjectionAlgorithm::initialize(const Config& _cfg)
 
 //----------------------------------------------------------------------------------------
 // Initialize - C++
-bool CCudaForwardProjectionAlgorithm::initialize(CProjectionGeometry2D* _pProjectionGeometry,
-												 CVolumeGeometry2D* _pReconstructionGeometry,
+bool CCudaForwardProjectionAlgorithm::initialize(CProjector2D* _pProjector,
 												 CFloat32VolumeData2D* _pVolume,
-												 CFloat32ProjectionData2D* _pSinogram,
-												 int _iGPUindex, int _iDetectorSuperSampling)
+												 CFloat32ProjectionData2D* _pSinogram)
 {
 	// store classes
-	//m_pProjectionGeometry = _pProjectionGeometry;
-	//m_pReconstructionGeometry = _pReconstructionGeometry;
+	m_pProjector = _pProjector;
 	m_pVolume = _pVolume;
 	m_pSinogram = _pSinogram;
 
-	m_iDetectorSuperSampling = _iDetectorSuperSampling;
-	m_iGPUIndex = _iGPUindex;
+	initializeFromProjector();
 
 	// return success
 	return check();
diff --git a/src/CudaForwardProjectionAlgorithm3D.cpp b/src/CudaForwardProjectionAlgorithm3D.cpp
index 46dab12..209f5a5 100644
--- a/src/CudaForwardProjectionAlgorithm3D.cpp
+++ b/src/CudaForwardProjectionAlgorithm3D.cpp
@@ -40,6 +40,8 @@ $Id$
 #include "astra/ParallelVecProjectionGeometry3D.h"
 #include "astra/ConeVecProjectionGeometry3D.h"
 
+#include "astra/CompositeGeometryManager.h"
+
 #include "astra/Logging.h"
 
 #include "../cuda/3d/astra3d.h"
@@ -72,6 +74,23 @@ CCudaForwardProjectionAlgorithm3D::~CCudaForwardProjectionAlgorithm3D()
 }
 
 //---------------------------------------------------------------------------------------
+void CCudaForwardProjectionAlgorithm3D::initializeFromProjector()
+{
+	m_iDetectorSuperSampling = 1;
+	m_iGPUIndex = -1;
+
+	CCudaProjector3D* pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector);
+	if (!pCudaProjector) {
+		if (m_pProjector) {
+			ASTRA_WARN("non-CUDA Projector3D passed to FP3D_CUDA");
+		}
+	} else {
+		m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling();
+		m_iGPUIndex = pCudaProjector->getGPUIndex();
+	}
+}
+
+//---------------------------------------------------------------------------------------
 // Initialize - Config
 bool CCudaForwardProjectionAlgorithm3D::initialize(const Config& _cfg)
 {
@@ -97,29 +116,21 @@ bool CCudaForwardProjectionAlgorithm3D::initialize(const Config& _cfg)
 
 	// optional: projector
 	node = _cfg.self.getSingleNode("ProjectorId");
-	CCudaProjector3D* pCudaProjector = 0;
 	m_pProjector = 0;
 	if (node) {
 		id = boost::lexical_cast<int>(node.getContent());
 		m_pProjector = CProjector3DManager::getSingleton().get(id);
-		pCudaProjector = dynamic_cast<CCudaProjector3D*>(CProjector3DManager::getSingleton().get(id));
-		m_pProjector = pCudaProjector;
-		if (!pCudaProjector) {
-			// TODO: Report
-		}
 	}
 	CC.markNodeParsed("ProjectorId");
 
-	// GPU number
-	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
-	CC.markOptionParsed("GPUindex");
-
+	initializeFromProjector();
 
-	m_iDetectorSuperSampling = 1;
-	if (pCudaProjector)
-		m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling();
+	// Deprecated options
 	m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex);
 	CC.markOptionParsed("DetectorSuperSampling");
+	CC.markOptionParsed("GPUindex");
+
 
 	// success
 	m_bIsInitialized = check();
@@ -142,8 +153,15 @@ bool CCudaForwardProjectionAlgorithm3D::initialize(CProjector3D* _pProjector,
 	m_pProjections = _pProjections;
 	m_pVolume = _pVolume;
 
-	m_iDetectorSuperSampling = _iDetectorSuperSampling;
-	m_iGPUIndex = _iGPUindex;
+	CCudaProjector3D* pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector);
+	if (!pCudaProjector) {
+		// TODO: Report
+		m_iDetectorSuperSampling = _iDetectorSuperSampling;
+		m_iGPUIndex = _iGPUindex;
+	} else {
+		m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling();
+		m_iGPUIndex = pCudaProjector->getGPUIndex();
+	}
 
 	// success
 	m_bIsInitialized = check();
@@ -247,11 +265,13 @@ void CCudaForwardProjectionAlgorithm3D::run(int)
 	// check initialized
 	assert(m_bIsInitialized);
 
+#if 1
+	CCompositeGeometryManager cgm;
+
+	cgm.doFP(m_pProjector, m_pVolume, m_pProjections);
+
+#else
 	const CProjectionGeometry3D* projgeom = m_pProjections->getGeometry();
-	const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(projgeom);
-	const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(projgeom);
-	const CConeVecProjectionGeometry3D* conevecgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(projgeom);
-	const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom);
 	const CVolumeGeometry3D& volgeom = *m_pVolume->getGeometry();
 
 	Cuda3DProjectionKernel projKernel = ker3d_default;
@@ -279,58 +299,10 @@ void CCudaForwardProjectionAlgorithm3D::run(int)
 	}
 #endif
 
-	if (conegeom) {
-		astraCudaConeFP(m_pVolume->getDataConst(), m_pProjections->getData(),
-		                volgeom.getGridColCount(),
-		                volgeom.getGridRowCount(),
-		                volgeom.getGridSliceCount(),
-		                conegeom->getProjectionCount(),
-		                conegeom->getDetectorColCount(),
-		                conegeom->getDetectorRowCount(),
-		                conegeom->getOriginSourceDistance(),
-		                conegeom->getOriginDetectorDistance(),
-		                conegeom->getDetectorSpacingX(),
-		                conegeom->getDetectorSpacingY(),
-		                conegeom->getProjectionAngles(),
-		                m_iGPUIndex, m_iDetectorSuperSampling);
-	} else if (par3dgeom) {
-		astraCudaPar3DFP(m_pVolume->getDataConst(), m_pProjections->getData(),
-		                 volgeom.getGridColCount(),
-		                 volgeom.getGridRowCount(),
-		                 volgeom.getGridSliceCount(),
-		                 par3dgeom->getProjectionCount(),
-		                 par3dgeom->getDetectorColCount(),
-		                 par3dgeom->getDetectorRowCount(),
-		                 par3dgeom->getDetectorSpacingX(),
-		                 par3dgeom->getDetectorSpacingY(),
-		                 par3dgeom->getProjectionAngles(),
-		                 m_iGPUIndex, m_iDetectorSuperSampling,
-		                 projKernel);
-	} else if (parvec3dgeom) {
-		astraCudaPar3DFP(m_pVolume->getDataConst(), m_pProjections->getData(),
-		                 volgeom.getGridColCount(),
-		                 volgeom.getGridRowCount(),
-		                 volgeom.getGridSliceCount(),
-		                 parvec3dgeom->getProjectionCount(),
-		                 parvec3dgeom->getDetectorColCount(),
-		                 parvec3dgeom->getDetectorRowCount(),
-		                 parvec3dgeom->getProjectionVectors(),
-		                 m_iGPUIndex, m_iDetectorSuperSampling,
-		                 projKernel);
-	} else if (conevecgeom) {
-		astraCudaConeFP(m_pVolume->getDataConst(), m_pProjections->getData(),
-		                volgeom.getGridColCount(),
-		                volgeom.getGridRowCount(),
-		                volgeom.getGridSliceCount(),
-		                conevecgeom->getProjectionCount(),
-		                conevecgeom->getDetectorColCount(),
-		                conevecgeom->getDetectorRowCount(),
-		                conevecgeom->getProjectionVectors(),
-		                m_iGPUIndex, m_iDetectorSuperSampling);
-	} else {
-		ASTRA_ASSERT(false);
-	}
-
+	astraCudaFP(m_pVolume->getDataConst(), m_pProjections->getData(),
+	            &volgeom, projgeom,
+	            m_iGPUIndex, m_iDetectorSuperSampling, projKernel);
+#endif
 }
 
 
diff --git a/src/CudaProjector2D.cpp b/src/CudaProjector2D.cpp
index a26e32d..acf6000 100644
--- a/src/CudaProjector2D.cpp
+++ b/src/CudaProjector2D.cpp
@@ -61,6 +61,7 @@ void CCudaProjector2D::_clear()
 	m_projectionKernel = ker2d_default;
 	m_iVoxelSuperSampling = 1;
 	m_iDetectorSuperSampling = 1;
+	m_iGPUIndex = -1;
 }
 
 //----------------------------------------------------------------------------------------
@@ -125,18 +126,18 @@ bool CCudaProjector2D::initialize(const Config& _cfg)
 	m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", 1);
 	CC.markOptionParsed("DetectorSuperSampling");
 
+	// GPU number
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
+	CC.markOptionParsed("GPUIndex");
+	if (!_cfg.self.hasOption("GPUIndex"))
+		CC.markOptionParsed("GPUindex");
+
+
 	m_bIsInitialized = _check();
 	return m_bIsInitialized;
 }
 
-/*
-bool CProjector2D::initialize(astra::CProjectionGeometry2D *, astra::CVolumeGeometry2D *)
-{
-	ASTRA_ASSERT(false);
-
-	return false;
-}
-*/
 
 std::string CCudaProjector2D::description() const
 {
diff --git a/src/CudaProjector3D.cpp b/src/CudaProjector3D.cpp
index d2fd74c..bbfbd34 100644
--- a/src/CudaProjector3D.cpp
+++ b/src/CudaProjector3D.cpp
@@ -64,6 +64,7 @@ void CCudaProjector3D::_clear()
 	m_projectionKernel = ker3d_default;
 	m_iVoxelSuperSampling = 1;
 	m_iDetectorSuperSampling = 1;
+	m_iGPUIndex = -1;
 }
 
 //----------------------------------------------------------------------------------------
@@ -128,6 +129,12 @@ bool CCudaProjector3D::initialize(const Config& _cfg)
 	m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", 1);
 	CC.markOptionParsed("DetectorSuperSampling");
 
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
+	CC.markOptionParsed("GPUIndex");
+	if (!_cfg.self.hasOption("GPUIndex"))
+		CC.markOptionParsed("GPUindex");
+
 	m_bIsInitialized = _check();
 	return m_bIsInitialized;
 }
diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp
index 71b6637..71dddf7 100644
--- a/src/CudaReconstructionAlgorithm2D.cpp
+++ b/src/CudaReconstructionAlgorithm2D.cpp
@@ -84,111 +84,51 @@ void CCudaReconstructionAlgorithm2D::_clear()
 }
 
 //---------------------------------------------------------------------------------------
-// Initialize - Config
-bool CCudaReconstructionAlgorithm2D::initialize(const Config& _cfg)
+void CCudaReconstructionAlgorithm2D::initializeFromProjector()
 {
-	ASTRA_ASSERT(_cfg.self);
-	ConfigStackCheck<CAlgorithm> CC("CudaReconstructionAlgorithm2D", this, _cfg);
-
-	// if already initialized, clear first
-	if (m_bIsInitialized) {
-		clear();
-	}
+	m_iPixelSuperSampling = 1;
+	m_iDetectorSuperSampling = 1;
+	m_iGPUIndex = -1;
 
 	// Projector
-	XMLNode node = _cfg.self.getSingleNode("ProjectorId");
-	CCudaProjector2D* pCudaProjector = 0;
-	if (node) {
-		int id = boost::lexical_cast<int>(node.getContent());
-		CProjector2D *projector = CProjector2DManager::getSingleton().get(id);
-		pCudaProjector = dynamic_cast<CCudaProjector2D*>(projector);
-		if (!pCudaProjector) {
+	CCudaProjector2D* pCudaProjector = dynamic_cast<CCudaProjector2D*>(m_pProjector);
+	if (!pCudaProjector) {
+		if (m_pProjector) {
 			ASTRA_WARN("non-CUDA Projector2D passed");
 		}
-	}
-	CC.markNodeParsed("ProjectorId");
-
-
-	// sinogram data
-	node = _cfg.self.getSingleNode("ProjectionDataId");
-	ASTRA_CONFIG_CHECK(node, "CudaSirt2", "No ProjectionDataId tag specified.");
-	int id = boost::lexical_cast<int>(node.getContent());
-	m_pSinogram = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id));
-	CC.markNodeParsed("ProjectionDataId");
-
-	// reconstruction data
-	node = _cfg.self.getSingleNode("ReconstructionDataId");
-	ASTRA_CONFIG_CHECK(node, "CudaSirt2", "No ReconstructionDataId tag specified.");
-	id = boost::lexical_cast<int>(node.getContent());
-	m_pReconstruction = dynamic_cast<CFloat32VolumeData2D*>(CData2DManager::getSingleton().get(id));
-	CC.markNodeParsed("ReconstructionDataId");
-
-	// fixed mask
-	if (_cfg.self.hasOption("ReconstructionMaskId")) {
-		m_bUseReconstructionMask = true;
-		id = boost::lexical_cast<int>(_cfg.self.getOption("ReconstructionMaskId"));
-		m_pReconstructionMask = dynamic_cast<CFloat32VolumeData2D*>(CData2DManager::getSingleton().get(id));
-		ASTRA_CONFIG_CHECK(m_pReconstructionMask, "CudaReconstruction2D", "Invalid ReconstructionMaskId.");
-	}
-	CC.markOptionParsed("ReconstructionMaskId");
-	// fixed mask
-	if (_cfg.self.hasOption("SinogramMaskId")) {
-		m_bUseSinogramMask = true;
-		id = boost::lexical_cast<int>(_cfg.self.getOption("SinogramMaskId"));
-		m_pSinogramMask = dynamic_cast<CFloat32ProjectionData2D*>(CData2DManager::getSingleton().get(id));
-		ASTRA_CONFIG_CHECK(m_pSinogramMask, "CudaReconstruction2D", "Invalid SinogramMaskId.");
-	}
-	CC.markOptionParsed("SinogramMaskId");
-
-	// Constraints - NEW
-	if (_cfg.self.hasOption("MinConstraint")) {
-		m_bUseMinConstraint = true;
-		m_fMinValue = _cfg.self.getOptionNumerical("MinConstraint", 0.0f);
-		CC.markOptionParsed("MinConstraint");
-	} else {
-		// Constraint - OLD
-		m_bUseMinConstraint = _cfg.self.getOptionBool("UseMinConstraint", false);
-		CC.markOptionParsed("UseMinConstraint");
-		if (m_bUseMinConstraint) {
-			m_fMinValue = _cfg.self.getOptionNumerical("MinConstraintValue", 0.0f);
-			CC.markOptionParsed("MinConstraintValue");
-		}
-	}
-	if (_cfg.self.hasOption("MaxConstraint")) {
-		m_bUseMaxConstraint = true;
-		m_fMaxValue = _cfg.self.getOptionNumerical("MaxConstraint", 255.0f);
-		CC.markOptionParsed("MaxConstraint");
 	} else {
-		// Constraint - OLD
-		m_bUseMaxConstraint = _cfg.self.getOptionBool("UseMaxConstraint", false);
-		CC.markOptionParsed("UseMaxConstraint");
-		if (m_bUseMaxConstraint) {
-			m_fMaxValue = _cfg.self.getOptionNumerical("MaxConstraintValue", 0.0f);
-			CC.markOptionParsed("MaxConstraintValue");
-		}
-	}
-
-	// GPU number
-	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
-	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
-	CC.markOptionParsed("GPUindex");
-	if (!_cfg.self.hasOption("GPUindex"))
-		CC.markOptionParsed("GPUIndex");
-
-	// Supersampling factors
-	m_iDetectorSuperSampling = 1;
-	m_iPixelSuperSampling = 1;
-	if (pCudaProjector) {
-		// New interface
 		m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling();
 		m_iPixelSuperSampling = pCudaProjector->getVoxelSuperSampling();
+		m_iGPUIndex = pCudaProjector->getGPUIndex();
 	}
+}
+
+//---------------------------------------------------------------------------------------
+// Initialize - Config
+bool CCudaReconstructionAlgorithm2D::initialize(const Config& _cfg)
+{
+	ASTRA_ASSERT(_cfg.self);
+	ConfigStackCheck<CAlgorithm> CC("CudaReconstructionAlgorithm2D", this, _cfg);
+
+	m_bIsInitialized = CReconstructionAlgorithm2D::initialize(_cfg);
+
+	if (!m_bIsInitialized)
+		return false;
+
+	initializeFromProjector();
+
 	// Deprecated options
 	m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling);
 	m_iPixelSuperSampling = (int)_cfg.self.getOptionNumerical("PixelSuperSampling", m_iPixelSuperSampling);
 	CC.markOptionParsed("DetectorSuperSampling");
 	CC.markOptionParsed("PixelSuperSampling");
 
+	// GPU number
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
+	CC.markOptionParsed("GPUIndex");
+	if (!_cfg.self.hasOption("GPUIndex"))
+		CC.markOptionParsed("GPUindex");
 
 	return _check();
 }
@@ -199,32 +139,18 @@ bool CCudaReconstructionAlgorithm2D::initialize(CProjector2D* _pProjector,
                                      CFloat32ProjectionData2D* _pSinogram, 
                                      CFloat32VolumeData2D* _pReconstruction)
 {
-	return initialize(_pProjector, _pSinogram, _pReconstruction, 0, 1);
-}
-
-//---------------------------------------------------------------------------------------
-// Initialize - C++
-bool CCudaReconstructionAlgorithm2D::initialize(CProjector2D* _pProjector,
-                                     CFloat32ProjectionData2D* _pSinogram, 
-                                     CFloat32VolumeData2D* _pReconstruction,
-                                     int _iGPUindex,
-                                     int _iDetectorSuperSampling,
-                                     int _iPixelSuperSampling)
-{
 	// if already initialized, clear first
 	if (m_bIsInitialized) {
 		clear();
 	}
 	
-	m_pProjector = 0;
+	m_pProjector = _pProjector;
 	
 	// required classes
 	m_pSinogram = _pSinogram;
 	m_pReconstruction = _pReconstruction;
 
-	m_iDetectorSuperSampling = _iDetectorSuperSampling;
-	m_iPixelSuperSampling = _iPixelSuperSampling;
-	m_iGPUIndex = _iGPUindex;
+	initializeFromProjector();
 
 	return _check();
 }
@@ -234,40 +160,13 @@ bool CCudaReconstructionAlgorithm2D::initialize(CProjector2D* _pProjector,
 // Check
 bool CCudaReconstructionAlgorithm2D::_check() 
 {
-	// TODO: CLEAN UP
-
-
-	// check pointers
-	//ASTRA_CONFIG_CHECK(m_pProjector, "Reconstruction2D", "Invalid Projector Object.");
-	ASTRA_CONFIG_CHECK(m_pSinogram, "SIRT_CUDA", "Invalid Projection Data Object.");
-	ASTRA_CONFIG_CHECK(m_pReconstruction, "SIRT_CUDA", "Invalid Reconstruction Data Object.");
-
-	// check initializations
-	//ASTRA_CONFIG_CHECK(m_pProjector->isInitialized(), "Reconstruction2D", "Projector Object Not Initialized.");
-	ASTRA_CONFIG_CHECK(m_pSinogram->isInitialized(), "SIRT_CUDA", "Projection Data Object Not Initialized.");
-	ASTRA_CONFIG_CHECK(m_pReconstruction->isInitialized(), "SIRT_CUDA", "Reconstruction Data Object Not Initialized.");
-
-	ASTRA_CONFIG_CHECK(m_iDetectorSuperSampling >= 1, "SIRT_CUDA", "DetectorSuperSampling must be a positive integer.");
-	ASTRA_CONFIG_CHECK(m_iPixelSuperSampling >= 1, "SIRT_CUDA", "PixelSuperSampling must be a positive integer.");
-	ASTRA_CONFIG_CHECK(m_iGPUIndex >= -1, "SIRT_CUDA", "GPUIndex must be a non-negative integer.");
-
-	// check compatibility between projector and data classes
-//	ASTRA_CONFIG_CHECK(m_pSinogram->getGeometry()->isEqual(m_pProjector->getProjectionGeometry()), "SIRT_CUDA", "Projection Data not compatible with the specified Projector.");
-//	ASTRA_CONFIG_CHECK(m_pReconstruction->getGeometry()->isEqual(m_pProjector->getVolumeGeometry()), "SIRT_CUDA", "Reconstruction Data not compatible with the specified Projector.");
-
-	// todo: turn some of these back on
+	if (!CReconstructionAlgorithm2D::_check())
+		return false;
 
-// 	ASTRA_CONFIG_CHECK(m_pProjectionGeometry, "SIRT_CUDA", "ProjectionGeometry not specified.");
-// 	ASTRA_CONFIG_CHECK(m_pProjectionGeometry->isInitialized(), "SIRT_CUDA", "ProjectionGeometry not initialized.");
-// 	ASTRA_CONFIG_CHECK(m_pReconstructionGeometry, "SIRT_CUDA", "ReconstructionGeometry not specified.");
-// 	ASTRA_CONFIG_CHECK(m_pReconstructionGeometry->isInitialized(), "SIRT_CUDA", "ReconstructionGeometry not initialized.");
+	ASTRA_CONFIG_CHECK(m_iDetectorSuperSampling >= 1, "CudaReconstructionAlgorithm2D", "DetectorSuperSampling must be a positive integer.");
+	ASTRA_CONFIG_CHECK(m_iPixelSuperSampling >= 1, "CudaReconstructionAlgorithm2D", "PixelSuperSampling must be a positive integer.");
+	ASTRA_CONFIG_CHECK(m_iGPUIndex >= -1, "CudaReconstructionAlgorithm2D", "GPUIndex must be a non-negative integer or -1.");
 
-	// check dimensions
-	//ASTRA_CONFIG_CHECK(m_pSinogram->getAngleCount() == m_pProjectionGeometry->getProjectionAngleCount(), "SIRT_CUDA", "Sinogram data object size mismatch.");
-	//ASTRA_CONFIG_CHECK(m_pSinogram->getDetectorCount() == m_pProjectionGeometry->getDetectorCount(), "SIRT_CUDA", "Sinogram data object size mismatch.");
-	//ASTRA_CONFIG_CHECK(m_pReconstruction->getWidth() == m_pReconstructionGeometry->getGridColCount(), "SIRT_CUDA", "Reconstruction data object size mismatch.");
-	//ASTRA_CONFIG_CHECK(m_pReconstruction->getHeight() == m_pReconstructionGeometry->getGridRowCount(), "SIRT_CUDA", "Reconstruction data object size mismatch.");
-	
 	// check restrictions
 	// TODO: check restrictions built into cuda code
 
@@ -462,10 +361,18 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations)
 
 	ASTRA_ASSERT(ok);
 
-	if (m_bUseMinConstraint)
-		ok &= m_pAlgo->setMinConstraint(m_fMinValue);
-	if (m_bUseMaxConstraint)
-		ok &= m_pAlgo->setMaxConstraint(m_fMaxValue);
+	if (m_bUseMinConstraint) {
+		bool ret = m_pAlgo->setMinConstraint(m_fMinValue);
+		if (!ret) {
+			ASTRA_WARN("This algorithm ignores MinConstraint");
+		}
+	}
+	if (m_bUseMaxConstraint) {
+		bool ret= m_pAlgo->setMaxConstraint(m_fMaxValue);
+		if (!ret) {
+			ASTRA_WARN("This algorithm ignores MaxConstraint");
+		}
+	}
 
 	ok &= m_pAlgo->iterate(_iNrIterations);
 	ASTRA_ASSERT(ok);
diff --git a/src/CudaSartAlgorithm.cpp b/src/CudaSartAlgorithm.cpp
index 8c0c6d7..d202847 100644
--- a/src/CudaSartAlgorithm.cpp
+++ b/src/CudaSartAlgorithm.cpp
@@ -116,10 +116,9 @@ bool CCudaSartAlgorithm::initialize(const Config& _cfg)
 // Initialize - C++
 bool CCudaSartAlgorithm::initialize(CProjector2D* _pProjector,
                                      CFloat32ProjectionData2D* _pSinogram, 
-                                     CFloat32VolumeData2D* _pReconstruction,
-                                     int _iGPUindex, int _iDetectorSuperSampling)
+                                     CFloat32VolumeData2D* _pReconstruction)
 {
-	m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, _iDetectorSuperSampling, 1);
+	m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction);
 
 	if (!m_bIsInitialized)
 		return false;
diff --git a/src/CudaSirtAlgorithm.cpp b/src/CudaSirtAlgorithm.cpp
index d424915..ab0a418 100644
--- a/src/CudaSirtAlgorithm.cpp
+++ b/src/CudaSirtAlgorithm.cpp
@@ -98,11 +98,9 @@ bool CCudaSirtAlgorithm::initialize(const Config& _cfg)
 // Initialize - C++
 bool CCudaSirtAlgorithm::initialize(CProjector2D* _pProjector,
                                      CFloat32ProjectionData2D* _pSinogram, 
-                                     CFloat32VolumeData2D* _pReconstruction,
-                                     int _iGPUindex, int _iDetectorSuperSampling,
-                                     int _iPixelSuperSampling)
+                                     CFloat32VolumeData2D* _pReconstruction)
 {
-	m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction, _iGPUindex, _iDetectorSuperSampling, _iPixelSuperSampling);
+	m_bIsInitialized = CCudaReconstructionAlgorithm2D::initialize(_pProjector, _pSinogram, _pReconstruction);
 
 	if (!m_bIsInitialized)
 		return false;
diff --git a/src/CudaSirtAlgorithm3D.cpp b/src/CudaSirtAlgorithm3D.cpp
index abbb9fd..1fa0da2 100644
--- a/src/CudaSirtAlgorithm3D.cpp
+++ b/src/CudaSirtAlgorithm3D.cpp
@@ -38,6 +38,8 @@ $Id$
 #include "astra/ConeVecProjectionGeometry3D.h"
 #include "astra/CudaProjector3D.h"
 
+#include "astra/Logging.h"
+
 #include "../cuda/3d/astra3d.h"
 
 using namespace std;
@@ -90,7 +92,27 @@ bool CCudaSirtAlgorithm3D::_check()
 	return true;
 }
 
-//---------------------------------------------------------------------------------------
+//----------------------------------------------------------------------------------------
+void CCudaSirtAlgorithm3D::initializeFromProjector()
+{
+	m_iVoxelSuperSampling = 1;
+	m_iDetectorSuperSampling = 1;
+	m_iGPUIndex = -1;
+
+	CCudaProjector3D* pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector);
+	if (!pCudaProjector) {
+		if (m_pProjector) {
+			ASTRA_WARN("non-CUDA Projector3D passed to SIRT3D_CUDA");
+		}
+	} else {
+		m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling();
+		m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling();
+		m_iGPUIndex = pCudaProjector->getGPUIndex();
+	}
+
+}
+
+//--------------------------------------------------------------------------------------
 // Initialize - Config
 bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg)
 {
@@ -108,28 +130,20 @@ bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg)
 		return false;
 	}
 
-	CCudaProjector3D* pCudaProjector = 0;
-	pCudaProjector = dynamic_cast<CCudaProjector3D*>(m_pProjector);
-	if (!pCudaProjector) {
-		// TODO: Report
-	}
+	initializeFromProjector();
 
-	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
-	CC.markOptionParsed("GPUindex");
-
-
-	m_iVoxelSuperSampling = 1;
-	m_iDetectorSuperSampling = 1;
-	if (pCudaProjector) {
-		// New interface
-		m_iVoxelSuperSampling = pCudaProjector->getVoxelSuperSampling();
-		m_iDetectorSuperSampling = pCudaProjector->getDetectorSuperSampling();
-	}
 	// Deprecated options
 	m_iVoxelSuperSampling = (int)_cfg.self.getOptionNumerical("VoxelSuperSampling", m_iVoxelSuperSampling);
 	m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex);
+	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
 	CC.markOptionParsed("VoxelSuperSampling");
 	CC.markOptionParsed("DetectorSuperSampling");
+	CC.markOptionParsed("GPUIndex");
+	if (!_cfg.self.hasOption("GPUIndex"))
+		CC.markOptionParsed("GPUindex");
+
+
 
 	m_pSirt = new AstraSIRT3d();
 
@@ -189,10 +203,6 @@ void CCudaSirtAlgorithm3D::run(int _iNrIterations)
 	ASTRA_ASSERT(m_bIsInitialized);
 
 	const CProjectionGeometry3D* projgeom = m_pSinogram->getGeometry();
-	const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(projgeom);
-	const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(projgeom);
-	const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(projgeom);
-	const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(projgeom);
 	const CVolumeGeometry3D& volgeom = *m_pReconstruction->getGeometry();
 
 	bool ok = true;
@@ -201,39 +211,7 @@ void CCudaSirtAlgorithm3D::run(int _iNrIterations)
 
 		ok &= m_pSirt->setGPUIndex(m_iGPUIndex);
 
-		ok &= m_pSirt->setReconstructionGeometry(volgeom.getGridColCount(),
-		                                         volgeom.getGridRowCount(),
-		                                         volgeom.getGridSliceCount());
-
-		if (conegeom) {
-			ok &= m_pSirt->setConeGeometry(conegeom->getProjectionCount(),
-			                               conegeom->getDetectorColCount(),
-			                               conegeom->getDetectorRowCount(),
-			                               conegeom->getOriginSourceDistance(),
-			                               conegeom->getOriginDetectorDistance(),
-			                               conegeom->getDetectorSpacingX(),
-			                               conegeom->getDetectorSpacingY(),
-			                               conegeom->getProjectionAngles());
-		} else if (par3dgeom) {
-			ok &= m_pSirt->setPar3DGeometry(par3dgeom->getProjectionCount(),
-			                                par3dgeom->getDetectorColCount(),
-			                                par3dgeom->getDetectorRowCount(),
-			                                par3dgeom->getDetectorSpacingX(),
-			                                par3dgeom->getDetectorSpacingY(),
-			                                par3dgeom->getProjectionAngles());
-		} else if (parvec3dgeom) {
-			ok &= m_pSirt->setPar3DGeometry(parvec3dgeom->getProjectionCount(),
-			                                parvec3dgeom->getDetectorColCount(),
-			                                parvec3dgeom->getDetectorRowCount(),
-			                                parvec3dgeom->getProjectionVectors());
-		} else if (conevec3dgeom) {
-			ok &= m_pSirt->setConeGeometry(conevec3dgeom->getProjectionCount(),
-			                               conevec3dgeom->getDetectorColCount(),
-			                               conevec3dgeom->getDetectorRowCount(),
-			                               conevec3dgeom->getProjectionVectors());
-		} else {
-			ASTRA_ASSERT(false);
-		}
+		ok &= m_pSirt->setGeometry(&volgeom, projgeom);
 
 		ok &= m_pSirt->enableSuperSampling(m_iVoxelSuperSampling, m_iDetectorSuperSampling);
 
diff --git a/src/Float32ProjectionData3D.cpp b/src/Float32ProjectionData3D.cpp
index 2bd0447..680ad55 100644
--- a/src/Float32ProjectionData3D.cpp
+++ b/src/Float32ProjectionData3D.cpp
@@ -53,13 +53,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const CFloat32Pro
 	CProjectionGeometry3D * pThisGeometry = getGeometry();
 
 	int iProjectionCount = pThisGeometry->getProjectionCount();
+	int iDetectorCount = pThisGeometry->getDetectorTotCount();
 #ifdef _DEBUG
 	CProjectionGeometry3D * pDataGeometry = _data.getGeometry();
-	int iThisProjectionDetectorCount = pThisGeometry->getDetectorRowCount() * pThisGeometry->getDetectorColCount();
-	int iDataProjectionDetectorCount = pDataGeometry->getDetectorRowCount() * pDataGeometry->getDetectorColCount();
+	int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount();
 
 	ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount());
-	ASTRA_ASSERT(iThisProjectionDetectorCount == iDataProjectionDetectorCount);
+	ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount);
 #endif
 
 	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
@@ -67,7 +67,7 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const CFloat32Pro
 		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
 		CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex);
 
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
+		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
 		{
 			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
 			float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex];
@@ -91,13 +91,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const CFloat32Pro
 	CProjectionGeometry3D * pThisGeometry = getGeometry();
 
 	int iProjectionCount = pThisGeometry->getProjectionCount();
+	int iDetectorCount = pThisGeometry->getDetectorTotCount();
 #ifdef _DEBUG
 	CProjectionGeometry3D * pDataGeometry = _data.getGeometry();
-	int iThisProjectionDetectorCount = pThisGeometry->getDetectorRowCount() * pThisGeometry->getDetectorColCount();
-	int iDataProjectionDetectorCount = pDataGeometry->getDetectorRowCount() * pDataGeometry->getDetectorColCount();
+	int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount();
 
 	ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount());
-	ASTRA_ASSERT(iThisProjectionDetectorCount == iDataProjectionDetectorCount);
+	ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount);
 #endif
 
 	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
@@ -105,7 +105,7 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const CFloat32Pro
 		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
 		CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex);
 
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
+		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
 		{
 			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
 			float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex];
@@ -129,13 +129,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const CFloat32Pro
 	CProjectionGeometry3D * pThisGeometry = getGeometry();
 
 	int iProjectionCount = pThisGeometry->getProjectionCount();
+	int iDetectorCount = pThisGeometry->getDetectorTotCount();
 #ifdef _DEBUG
 	CProjectionGeometry3D * pDataGeometry = _data.getGeometry();
-	int iThisProjectionDetectorCount = pThisGeometry->getDetectorRowCount() * pThisGeometry->getDetectorColCount();
-	int iDataProjectionDetectorCount = pDataGeometry->getDetectorRowCount() * pDataGeometry->getDetectorColCount();
+	int iDataProjectionDetectorCount = pDataGeometry->getDetectorTotCount();
 
 	ASTRA_ASSERT(iProjectionCount == pDataGeometry->getProjectionCount());
-	ASTRA_ASSERT(iThisProjectionDetectorCount == iDataProjectionDetectorCount);
+	ASTRA_ASSERT(iDetectorCount == iDataProjectionDetectorCount);
 #endif
 
 	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
@@ -143,7 +143,7 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const CFloat32Pro
 		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
 		CFloat32VolumeData2D * pDataProjection = _data.fetchProjection(iProjectionIndex);
 
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
+		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
 		{
 			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
 			float32 fDataValue = pDataProjection->getDataConst()[iDetectorIndex];
@@ -167,12 +167,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator*=(const float32& _f
 	CProjectionGeometry3D * pThisGeometry = getGeometry();
 
 	int iProjectionCount = pThisGeometry->getProjectionCount();
+	int iDetectorCount = pThisGeometry->getDetectorTotCount();
 
 	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
 	{
 		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
 
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
+		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
 		{
 			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
 
@@ -194,12 +195,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator/=(const float32& _f
 	CProjectionGeometry3D * pThisGeometry = getGeometry();
 
 	int iProjectionCount = pThisGeometry->getProjectionCount();
+	int iDetectorCount = pThisGeometry->getDetectorTotCount();
 
 	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
 	{
 		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
 
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
+		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
 		{
 			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
 
@@ -221,12 +223,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator+=(const float32& _f
 	CProjectionGeometry3D * pThisGeometry = getGeometry();
 
 	int iProjectionCount = pThisGeometry->getProjectionCount();
+	int iDetectorCount = pThisGeometry->getDetectorTotCount();
 
 	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
 	{
 		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
 
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
+		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
 		{
 			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
 
@@ -248,12 +251,13 @@ CFloat32ProjectionData3D& CFloat32ProjectionData3D::operator-=(const float32& _f
 	CProjectionGeometry3D * pThisGeometry = getGeometry();
 
 	int iProjectionCount = pThisGeometry->getProjectionCount();
+	int iDetectorCount = pThisGeometry->getDetectorTotCount();
 
 	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
 	{
 		CFloat32VolumeData2D * pThisProjection = fetchProjection(iProjectionIndex);
 
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorIndex; iDetectorIndex++)
+		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
 		{
 			float32 fThisValue = pThisProjection->getData()[iDetectorIndex];
 
diff --git a/src/GeometryUtil3D.cpp b/src/GeometryUtil3D.cpp
index 52dd5a9..c6bfd8b 100644
--- a/src/GeometryUtil3D.cpp
+++ b/src/GeometryUtil3D.cpp
@@ -28,8 +28,96 @@ $Id$
 
 #include "astra/GeometryUtil3D.h"
 
+#include <cmath>
+
 namespace astra {
 
+
+SConeProjection* genConeProjections(unsigned int iProjAngles,
+                                    unsigned int iProjU,
+                                    unsigned int iProjV,
+                                    double fOriginSourceDistance,
+                                    double fOriginDetectorDistance,
+                                    double fDetUSize,
+                                    double fDetVSize,
+                                    const float *pfAngles)
+{
+	SConeProjection base;
+	base.fSrcX = 0.0f;
+	base.fSrcY = -fOriginSourceDistance;
+	base.fSrcZ = 0.0f;
+
+	base.fDetSX = iProjU * fDetUSize * -0.5f;
+	base.fDetSY = fOriginDetectorDistance;
+	base.fDetSZ = iProjV * fDetVSize * -0.5f;
+
+	base.fDetUX = fDetUSize;
+	base.fDetUY = 0.0f;
+	base.fDetUZ = 0.0f;
+
+	base.fDetVX = 0.0f;
+	base.fDetVY = 0.0f;
+	base.fDetVZ = fDetVSize;
+
+	SConeProjection* p = new SConeProjection[iProjAngles];
+
+#define ROTATE0(name,i,alpha) do { p[i].f##name##X = base.f##name##X * cos(alpha) - base.f##name##Y * sin(alpha); p[i].f##name##Y = base.f##name##X * sin(alpha) + base.f##name##Y * cos(alpha); p[i].f##name##Z = base.f##name##Z; } while(0)
+
+	for (unsigned int i = 0; i < iProjAngles; ++i) {
+		ROTATE0(Src, i, pfAngles[i]);
+		ROTATE0(DetS, i, pfAngles[i]);
+		ROTATE0(DetU, i, pfAngles[i]);
+		ROTATE0(DetV, i, pfAngles[i]);
+	}
+
+#undef ROTATE0
+
+	return p;
+}
+
+SPar3DProjection* genPar3DProjections(unsigned int iProjAngles,
+                                      unsigned int iProjU,
+                                      unsigned int iProjV,
+                                      double fDetUSize,
+                                      double fDetVSize,
+                                      const float *pfAngles)
+{
+	SPar3DProjection base;
+	base.fRayX = 0.0f;
+	base.fRayY = 1.0f;
+	base.fRayZ = 0.0f;
+
+	base.fDetSX = iProjU * fDetUSize * -0.5f;
+	base.fDetSY = 0.0f;
+	base.fDetSZ = iProjV * fDetVSize * -0.5f;
+
+	base.fDetUX = fDetUSize;
+	base.fDetUY = 0.0f;
+	base.fDetUZ = 0.0f;
+
+	base.fDetVX = 0.0f;
+	base.fDetVY = 0.0f;
+	base.fDetVZ = fDetVSize;
+
+	SPar3DProjection* p = new SPar3DProjection[iProjAngles];
+
+#define ROTATE0(name,i,alpha) do { p[i].f##name##X = base.f##name##X * cos(alpha) - base.f##name##Y * sin(alpha); p[i].f##name##Y = base.f##name##X * sin(alpha) + base.f##name##Y * cos(alpha); p[i].f##name##Z = base.f##name##Z; } while(0)
+
+	for (unsigned int i = 0; i < iProjAngles; ++i) {
+		ROTATE0(Ray, i, pfAngles[i]);
+		ROTATE0(DetS, i, pfAngles[i]);
+		ROTATE0(DetU, i, pfAngles[i]);
+		ROTATE0(DetV, i, pfAngles[i]);
+	}
+
+#undef ROTATE0
+
+	return p;
+}
+
+
+
+
 // (See declaration in header for (mathematical) description of these functions)
 
 
@@ -72,4 +160,88 @@ void computeBP_UV_Coeffs(const SConeProjection& proj, double &fUX, double &fUY,
 }
 
 
+// TODO: Handle cases of rays parallel to coordinate planes
+
+void backprojectPointX(const SPar3DProjection& proj, double fU, double fV,
+                       double fX, double &fY, double &fZ)
+{
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fX - px) / proj.fRayX;
+
+	fY = py + a * proj.fRayY;
+	fZ = pz + a * proj.fRayZ;
+}
+
+void backprojectPointY(const SPar3DProjection& proj, double fU, double fV,
+                       double fY, double &fX, double &fZ)
+{
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fY - py) / proj.fRayY;
+
+	fX = px + a * proj.fRayX;
+	fZ = pz + a * proj.fRayZ;
+
+}
+
+void backprojectPointZ(const SPar3DProjection& proj, double fU, double fV,
+                       double fZ, double &fX, double &fY)
+{
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fZ - pz) / proj.fRayZ;
+
+	fX = px + a * proj.fRayX;
+	fY = py + a * proj.fRayY;
+}
+
+
+
+void backprojectPointX(const SConeProjection& proj, double fU, double fV,
+                       double fX, double &fY, double &fZ)
+{
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fX - proj.fSrcX) / (px - proj.fSrcX);
+
+	fY = proj.fSrcY + a * (py - proj.fSrcY);
+	fZ = proj.fSrcZ + a * (pz - proj.fSrcZ);
+}
+
+void backprojectPointY(const SConeProjection& proj, double fU, double fV,
+                       double fY, double &fX, double &fZ)
+{
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fY - proj.fSrcY) / (py - proj.fSrcY);
+
+	fX = proj.fSrcX + a * (px - proj.fSrcX);
+	fZ = proj.fSrcZ + a * (pz - proj.fSrcZ);
+}
+
+void backprojectPointZ(const SConeProjection& proj, double fU, double fV,
+                       double fZ, double &fX, double &fY)
+{
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fZ - proj.fSrcZ) / (pz - proj.fSrcZ);
+
+	fX = proj.fSrcX + a * (px - proj.fSrcX);
+	fY = proj.fSrcY + a * (py - proj.fSrcY);
+}
+
+
 }
diff --git a/src/Globals.cpp b/src/Globals.cpp
index 813f9c9..904a459 100644
--- a/src/Globals.cpp
+++ b/src/Globals.cpp
@@ -28,5 +28,8 @@ $Id$
 
 #include "astra/Globals.h"
 
+namespace astra{
+    bool running_in_matlab=false;
+}
 // nothing to see here :)
 
diff --git a/src/Logging.cpp b/src/Logging.cpp
index 8290ca0..cd7e3f0 100644
--- a/src/Logging.cpp
+++ b/src/Logging.cpp
@@ -70,37 +70,65 @@ void CLogger::disable()
 void CLogger::debug(const char *sfile, int sline, const char *fmt, ...)
 {
 	_assureIsInitialized();
-	va_list ap;
-	va_start(ap, fmt);
-	if(m_bEnabledScreen) clog_debug(sfile,sline,0,fmt,ap);
-	if(m_bEnabledFile && m_bFileProvided) clog_debug(sfile,sline,1,fmt,ap);
+	va_list ap, apf;
+	if(m_bEnabledScreen){
+        va_start(ap, fmt);
+        clog_debug(sfile,sline,0,fmt,ap);
+        va_end(ap);
+    }
+	if(m_bEnabledFile && m_bFileProvided){
+        va_start(apf, fmt);
+        clog_debug(sfile,sline,1,fmt,apf);
+        va_end(apf);
+    }
 }
 
 void CLogger::info(const char *sfile, int sline, const char *fmt, ...)
 {
 	_assureIsInitialized();
-	va_list ap;
-	va_start(ap, fmt);
-	if(m_bEnabledScreen) clog_info(sfile,sline,0,fmt,ap);
-	if(m_bEnabledFile && m_bFileProvided) clog_info(sfile,sline,1,fmt,ap);
+	va_list ap, apf;
+	if(m_bEnabledScreen){
+        va_start(ap, fmt);
+        clog_info(sfile,sline,0,fmt,ap);
+        va_end(ap);
+    }
+	if(m_bEnabledFile && m_bFileProvided){
+        va_start(apf, fmt);
+        clog_info(sfile,sline,1,fmt,apf);
+        va_end(apf);
+    }
 }
 
 void CLogger::warn(const char *sfile, int sline, const char *fmt, ...)
 {
 	_assureIsInitialized();
-	va_list ap;
-	va_start(ap, fmt);
-	if(m_bEnabledScreen) clog_warn(sfile,sline,0,fmt,ap);
-	if(m_bEnabledFile && m_bFileProvided) clog_warn(sfile,sline,1,fmt,ap);
+	va_list ap, apf;
+	if(m_bEnabledScreen){
+        va_start(ap, fmt);
+        clog_warn(sfile,sline,0,fmt,ap);
+        va_end(ap);
+    }
+	if(m_bEnabledFile && m_bFileProvided){
+        va_start(apf, fmt);
+        clog_warn(sfile,sline,1,fmt,apf);
+        va_end(apf);
+    }
 }
 
 void CLogger::error(const char *sfile, int sline, const char *fmt, ...)
 {
 	_assureIsInitialized();
-	va_list ap;
-	va_start(ap, fmt);
-	if(m_bEnabledScreen) clog_error(sfile,sline,0,fmt,ap);
-	if(m_bEnabledFile && m_bFileProvided) clog_error(sfile,sline,1,fmt,ap);
+	va_list ap, apf;
+	if(m_bEnabledScreen){
+        va_start(ap, fmt);
+        clog_error(sfile,sline,0,fmt,ap);
+        va_end(ap);
+    }
+	if(m_bEnabledFile && m_bFileProvided){
+        va_start(apf, fmt);
+        clog_error(sfile,sline,1,fmt,apf);
+        va_end(apf);
+    }
 }
 
 void CLogger::_setLevel(int id, log_level m_eLevel)
diff --git a/src/ParallelProjectionGeometry3D.cpp b/src/ParallelProjectionGeometry3D.cpp
index 1c87157..7b64fd9 100644
--- a/src/ParallelProjectionGeometry3D.cpp
+++ b/src/ParallelProjectionGeometry3D.cpp
@@ -27,8 +27,10 @@ $Id$
 */
 
 #include "astra/ParallelProjectionGeometry3D.h"
-#include <boost/lexical_cast.hpp>
 
+#include "astra/GeometryUtil3D.h"
+
+#include <boost/lexical_cast.hpp>
 #include <cstring>
 
 using namespace std;
@@ -185,9 +187,9 @@ CVector3D CParallelProjectionGeometry3D::getProjectionDirection(int _iProjection
 	return CVector3D(fDirX, fDirY, fDirZ);
 }
 
-void CParallelProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ,
+void CParallelProjectionGeometry3D::projectPoint(double fX, double fY, double fZ,
                                                  int iAngleIndex,
-                                                 float32 &fU, float32 &fV) const
+                                                 double &fU, double &fV) const
 {
 	ASTRA_ASSERT(iAngleIndex >= 0);
 	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
@@ -214,6 +216,79 @@ CParallelProjectionGeometry2D * CParallelProjectionGeometry3D::createProjectionG
 	return pOutput;
 }
 
+void CParallelProjectionGeometry3D::backprojectPointX(int iAngleIndex, double fU, double fV,
+	                               double fX, double &fY, double &fZ) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	SPar3DProjection *projs = genPar3DProjections(1, m_iDetectorColCount, m_iDetectorRowCount,
+	                                           m_fDetectorSpacingX, m_fDetectorSpacingY,
+	                                           &m_pfProjectionAngles[iAngleIndex]);
+
+	SPar3DProjection &proj = projs[0];
+
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fX - px) / proj.fRayX;
+
+	fY = py + a * proj.fRayY;
+	fZ = pz + a * proj.fRayZ;
+
+	delete[] projs;
+}
+
+void CParallelProjectionGeometry3D::backprojectPointY(int iAngleIndex, double fU, double fV,
+	                               double fY, double &fX, double &fZ) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	SPar3DProjection *projs = genPar3DProjections(1, m_iDetectorColCount, m_iDetectorRowCount,
+	                                           m_fDetectorSpacingX, m_fDetectorSpacingY,
+	                                           &m_pfProjectionAngles[iAngleIndex]);
+
+	SPar3DProjection &proj = projs[0];
+
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fY - py) / proj.fRayY;
+
+	fX = px + a * proj.fRayX;
+	fZ = pz + a * proj.fRayZ;
+
+	delete[] projs;
+}
+
+void CParallelProjectionGeometry3D::backprojectPointZ(int iAngleIndex, double fU, double fV,
+	                               double fZ, double &fX, double &fY) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	SPar3DProjection *projs = genPar3DProjections(1, m_iDetectorColCount, m_iDetectorRowCount,
+	                                           m_fDetectorSpacingX, m_fDetectorSpacingY,
+	                                           &m_pfProjectionAngles[iAngleIndex]);
+
+	SPar3DProjection &proj = projs[0];
+
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fZ - pz) / proj.fRayZ;
+
+	fX = px + a * proj.fRayX;
+	fY = py + a * proj.fRayY;
+
+	delete[] projs;
+}
+
+
 //----------------------------------------------------------------------------------------
 
 } // end namespace astra
diff --git a/src/ParallelVecProjectionGeometry3D.cpp b/src/ParallelVecProjectionGeometry3D.cpp
index ffad6d0..d04400b 100644
--- a/src/ParallelVecProjectionGeometry3D.cpp
+++ b/src/ParallelVecProjectionGeometry3D.cpp
@@ -239,9 +239,9 @@ CVector3D CParallelVecProjectionGeometry3D::getProjectionDirection(int _iProject
 	return CVector3D(p.fRayX, p.fRayY, p.fRayZ);
 }
 
-void CParallelVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ,
-                                                 int iAngleIndex,
-                                                 float32 &fU, float32 &fV) const
+void CParallelVecProjectionGeometry3D::projectPoint(double fX, double fY, double fZ,
+                                                    int iAngleIndex,
+                                                    double &fU, double &fV) const
 {
 	ASTRA_ASSERT(iAngleIndex >= 0);
 	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
@@ -258,6 +258,61 @@ void CParallelVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, floa
 
 }
 
+void CParallelVecProjectionGeometry3D::backprojectPointX(int iAngleIndex, double fU, double fV,
+	                               double fX, double &fY, double &fZ) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	SPar3DProjection &proj = m_pProjectionAngles[iAngleIndex];
+
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fX - px) / proj.fRayX;
+
+	fY = py + a * proj.fRayY;
+	fZ = pz + a * proj.fRayZ;
+}
+
+void CParallelVecProjectionGeometry3D::backprojectPointY(int iAngleIndex, double fU, double fV,
+	                               double fY, double &fX, double &fZ) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	SPar3DProjection &proj = m_pProjectionAngles[iAngleIndex];
+
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fY - py) / proj.fRayY;
+
+	fX = px + a * proj.fRayX;
+	fZ = pz + a * proj.fRayZ;
+}
+
+void CParallelVecProjectionGeometry3D::backprojectPointZ(int iAngleIndex, double fU, double fV,
+	                               double fZ, double &fX, double &fY) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	SPar3DProjection &proj = m_pProjectionAngles[iAngleIndex];
+
+	double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX;
+	double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY;
+	double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ;
+
+	double a = (fZ - pz) / proj.fRayZ;
+
+	fX = px + a * proj.fRayX;
+	fY = py + a * proj.fRayY;
+}
+
+
 //----------------------------------------------------------------------------------------
 
 bool CParallelVecProjectionGeometry3D::_check()
diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp
new file mode 100644
index 0000000..8f7dfc5
--- /dev/null
+++ b/src/PluginAlgorithm.cpp
@@ -0,0 +1,401 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+           2014-2015, CWI, Amsterdam
+
+Contact: astra@uantwerpen.be
+Website: http://sf.net/projects/astra-toolbox
+
+This file is part of the ASTRA Toolbox.
+
+
+The ASTRA Toolbox is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+The ASTRA Toolbox is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+$Id$
+*/
+
+#ifdef ASTRA_PYTHON
+
+#include "astra/PluginAlgorithm.h"
+#include "astra/Logging.h"
+#include <boost/algorithm/string.hpp>
+#include <boost/algorithm/string/split.hpp>
+#include <boost/lexical_cast.hpp>
+#include <iostream>
+#include <fstream>
+#include <string>
+
+#include <Python.h>
+#include "bytesobject.h"
+
+namespace astra {
+
+
+
+void logPythonError(){
+    if(PyErr_Occurred()){
+        PyObject *ptype, *pvalue, *ptraceback;
+        PyErr_Fetch(&ptype, &pvalue, &ptraceback);
+        PyErr_NormalizeException(&ptype, &pvalue, &ptraceback);
+        PyObject *traceback = PyImport_ImportModule("traceback");
+        if(traceback!=NULL){
+            PyObject *exc;
+            if(ptraceback==NULL){
+                exc = PyObject_CallMethod(traceback,"format_exception_only","OO",ptype, pvalue);
+            }else{
+                exc = PyObject_CallMethod(traceback,"format_exception","OOO",ptype, pvalue, ptraceback);
+            }
+            if(exc!=NULL){
+                PyObject *six = PyImport_ImportModule("six");
+                if(six!=NULL){
+                    PyObject *iter = PyObject_GetIter(exc);
+                    if(iter!=NULL){
+                        PyObject *line;
+                        std::string errStr = "";
+                        while(line = PyIter_Next(iter)){
+                            PyObject *retb = PyObject_CallMethod(six,"b","O",line);
+                            if(retb!=NULL){
+                                errStr += std::string(PyBytes_AsString(retb));
+                                Py_DECREF(retb);
+                            }
+                            Py_DECREF(line);
+                        }
+                        ASTRA_ERROR("%s",errStr.c_str());
+                        Py_DECREF(iter);
+                    }
+                    Py_DECREF(six);
+                }
+                Py_DECREF(exc);
+            }
+            Py_DECREF(traceback);
+        }
+        if(ptype!=NULL) Py_DECREF(ptype);
+        if(pvalue!=NULL) Py_DECREF(pvalue);
+        if(ptraceback!=NULL) Py_DECREF(ptraceback);
+    }
+}
+
+
+CPluginAlgorithm::CPluginAlgorithm(PyObject* pyclass){
+    instance = PyObject_CallObject(pyclass, NULL);
+    if(instance==NULL) logPythonError();
+}
+
+CPluginAlgorithm::~CPluginAlgorithm(){
+    if(instance!=NULL){
+        Py_DECREF(instance);
+        instance = NULL;
+    }
+}
+
+bool CPluginAlgorithm::initialize(const Config& _cfg){
+    if(instance==NULL) return false;
+    PyObject *cfgDict = XMLNode2dict(_cfg.self);
+    PyObject *retVal = PyObject_CallMethod(instance, "astra_init", "O",cfgDict);
+    Py_DECREF(cfgDict);
+    if(retVal==NULL){
+        logPythonError();
+        return false;
+    }
+    m_bIsInitialized = true;
+    Py_DECREF(retVal);
+    return m_bIsInitialized;
+}
+
+void CPluginAlgorithm::run(int _iNrIterations){
+    if(instance==NULL) return;
+    PyGILState_STATE state = PyGILState_Ensure();
+    PyObject *retVal = PyObject_CallMethod(instance, "run", "i",_iNrIterations);
+    if(retVal==NULL){
+        logPythonError();
+    }else{
+        Py_DECREF(retVal);
+    }
+    PyGILState_Release(state);
+}
+
+void fixLapackLoading(){
+    // When running in Matlab, we need to force numpy
+    // to use its internal lapack library instead of
+    // Matlab's MKL library to avoid errors. To do this,
+    // we set Python's dlopen flags to RTLD_NOW|RTLD_DEEPBIND
+    // and import 'numpy.linalg.lapack_lite' here. We reset
+    // Python's dlopen flags afterwards.
+    PyObject *sys = PyImport_ImportModule("sys");
+    if(sys!=NULL){
+        PyObject *curFlags = PyObject_CallMethod(sys,"getdlopenflags",NULL);
+        if(curFlags!=NULL){
+            PyObject *retVal = PyObject_CallMethod(sys, "setdlopenflags", "i",10);
+            if(retVal!=NULL){
+                PyObject *lapack = PyImport_ImportModule("numpy.linalg.lapack_lite");
+                if(lapack!=NULL){
+                    Py_DECREF(lapack);
+                }
+                PyObject_CallMethod(sys, "setdlopenflags", "O",curFlags);
+                Py_DECREF(retVal);
+            }
+            Py_DECREF(curFlags);
+        }
+        Py_DECREF(sys);
+    }
+}
+
+CPluginAlgorithmFactory::CPluginAlgorithmFactory(){
+    if(!Py_IsInitialized()){
+        Py_Initialize();
+        PyEval_InitThreads();
+    }
+#ifndef _MSC_VER
+    if(astra::running_in_matlab) fixLapackLoading();
+#endif
+    pluginDict = PyDict_New();
+    inspect = PyImport_ImportModule("inspect");
+    six = PyImport_ImportModule("six");
+}
+
+CPluginAlgorithmFactory::~CPluginAlgorithmFactory(){
+    if(pluginDict!=NULL){
+        Py_DECREF(pluginDict);
+    }
+    if(inspect!=NULL) Py_DECREF(inspect);
+    if(six!=NULL) Py_DECREF(six);
+}
+
+PyObject * getClassFromString(std::string str){
+    std::vector<std::string> items;
+    boost::split(items, str, boost::is_any_of("."));
+    PyObject *pyclass = PyImport_ImportModule(items[0].c_str());
+    if(pyclass==NULL){
+        logPythonError();
+        return NULL;
+    }
+    PyObject *submod = pyclass;
+    for(unsigned int i=1;i<items.size();i++){
+        submod = PyObject_GetAttrString(submod,items[i].c_str());
+        Py_DECREF(pyclass);
+        pyclass = submod;
+        if(pyclass==NULL){
+            logPythonError();
+            return NULL;
+        }
+    }
+    return pyclass;
+}
+
+bool CPluginAlgorithmFactory::registerPlugin(std::string name, std::string className){
+    PyObject *str = PyBytes_FromString(className.c_str());
+    PyDict_SetItemString(pluginDict, name.c_str(), str);
+    Py_DECREF(str);
+    return true;
+}
+
+bool CPluginAlgorithmFactory::registerPlugin(std::string className){
+    PyObject *pyclass = getClassFromString(className);
+    if(pyclass==NULL) return false;
+    bool ret = registerPluginClass(pyclass);
+    Py_DECREF(pyclass);
+    return ret;
+}
+
+bool CPluginAlgorithmFactory::registerPluginClass(std::string name, PyObject * className){
+    PyDict_SetItemString(pluginDict, name.c_str(), className);
+    return true;
+}
+
+bool CPluginAlgorithmFactory::registerPluginClass(PyObject * className){
+    PyObject *astra_name = PyObject_GetAttrString(className,"astra_name");
+    if(astra_name==NULL){
+        logPythonError();
+        return false;
+    }
+    PyObject *retb = PyObject_CallMethod(six,"b","O",astra_name);
+    if(retb!=NULL){
+        PyDict_SetItemString(pluginDict,PyBytes_AsString(retb),className);
+        Py_DECREF(retb);
+    }else{
+        logPythonError();
+    }
+    Py_DECREF(astra_name);
+    return true;
+}
+
+CPluginAlgorithm * CPluginAlgorithmFactory::getPlugin(std::string name){
+    PyObject *className = PyDict_GetItemString(pluginDict, name.c_str());
+    if(className==NULL) return NULL;
+    CPluginAlgorithm *alg = NULL;
+    if(PyBytes_Check(className)){
+        std::string str = std::string(PyBytes_AsString(className));
+    	PyObject *pyclass = getClassFromString(str);
+        if(pyclass!=NULL){
+            alg = new CPluginAlgorithm(pyclass);
+            Py_DECREF(pyclass);
+        }
+    }else{
+        alg = new CPluginAlgorithm(className);
+    }
+    return alg;
+}
+
+PyObject * CPluginAlgorithmFactory::getRegistered(){
+    Py_INCREF(pluginDict);
+    return pluginDict;
+}
+
+std::map<std::string, std::string> CPluginAlgorithmFactory::getRegisteredMap(){
+    std::map<std::string, std::string> ret;
+    PyObject *key, *value;
+    Py_ssize_t pos = 0;
+    while (PyDict_Next(pluginDict, &pos, &key, &value)) {
+        PyObject *keystr = PyObject_Str(key);
+        if(keystr!=NULL){
+            PyObject *valstr = PyObject_Str(value);
+            if(valstr!=NULL){
+                PyObject * keyb = PyObject_CallMethod(six,"b","O",keystr);
+                if(keyb!=NULL){
+                    PyObject * valb = PyObject_CallMethod(six,"b","O",valstr);
+                    if(valb!=NULL){
+                        ret[PyBytes_AsString(keyb)] = PyBytes_AsString(valb);
+                        Py_DECREF(valb);
+                    }
+                    Py_DECREF(keyb);
+                }
+                Py_DECREF(valstr);
+            }
+            Py_DECREF(keystr);
+        }
+        logPythonError();
+    }
+    return ret;
+}
+
+std::string CPluginAlgorithmFactory::getHelp(std::string name){
+    PyObject *className = PyDict_GetItemString(pluginDict, name.c_str());
+    if(className==NULL){
+        ASTRA_ERROR("Plugin %s not found!",name.c_str());
+        PyErr_Clear();
+        return "";
+    }
+    std::string ret = "";
+    PyObject *pyclass;
+    if(PyBytes_Check(className)){
+        std::string str = std::string(PyBytes_AsString(className));
+        pyclass = getClassFromString(str);
+    }else{
+        pyclass = className;
+    }
+    if(pyclass==NULL) return "";
+    if(inspect!=NULL && six!=NULL){
+        PyObject *retVal = PyObject_CallMethod(inspect,"getdoc","O",pyclass);
+        if(retVal!=NULL){
+            if(retVal!=Py_None){
+                PyObject *retb = PyObject_CallMethod(six,"b","O",retVal);
+                if(retb!=NULL){
+                    ret = std::string(PyBytes_AsString(retb));
+                    Py_DECREF(retb);
+                }
+            }
+            Py_DECREF(retVal);
+        }else{
+            logPythonError();
+        }
+    }
+    if(PyBytes_Check(className)){
+        Py_DECREF(pyclass);
+    }
+    return ret;
+}
+
+DEFINE_SINGLETON(CPluginAlgorithmFactory);
+
+#if PY_MAJOR_VERSION >= 3
+PyObject * pyStringFromString(std::string str){
+    return PyUnicode_FromString(str.c_str());
+}
+#else
+PyObject * pyStringFromString(std::string str){
+    return PyBytes_FromString(str.c_str());
+}
+#endif
+
+PyObject* stringToPythonValue(std::string str){
+    if(str.find(";")!=std::string::npos){
+        std::vector<std::string> rows, row;
+        boost::split(rows, str, boost::is_any_of(";"));
+        PyObject *mat = PyList_New(rows.size());
+        for(unsigned int i=0; i<rows.size(); i++){
+            boost::split(row, rows[i], boost::is_any_of(","));
+            PyObject *rowlist = PyList_New(row.size());
+            for(unsigned int j=0;j<row.size();j++){
+                PyList_SetItem(rowlist, j, PyFloat_FromDouble(boost::lexical_cast<double>(row[j])));
+            }
+            PyList_SetItem(mat, i, rowlist);
+        }
+        return mat;
+    }
+    if(str.find(",")!=std::string::npos){
+        std::vector<std::string> vec;
+        boost::split(vec, str, boost::is_any_of(","));
+        PyObject *veclist = PyList_New(vec.size());
+        for(unsigned int i=0;i<vec.size();i++){
+            PyList_SetItem(veclist, i, PyFloat_FromDouble(boost::lexical_cast<double>(vec[i])));
+        }
+        return veclist;
+    }
+    try{
+        return PyLong_FromLong(boost::lexical_cast<long>(str));
+    }catch(const boost::bad_lexical_cast &){
+        try{
+            return PyFloat_FromDouble(boost::lexical_cast<double>(str));
+        }catch(const boost::bad_lexical_cast &){
+            return pyStringFromString(str);
+        }
+    }
+}
+
+PyObject* XMLNode2dict(XMLNode node){
+    PyObject *dct = PyDict_New();
+    PyObject *opts = PyDict_New();
+    if(node.hasAttribute("type")){
+        PyObject *obj = pyStringFromString(node.getAttribute("type").c_str());
+        PyDict_SetItemString(dct, "type", obj);
+        Py_DECREF(obj);
+    }
+    std::list<XMLNode> nodes = node.getNodes();
+    std::list<XMLNode>::iterator it = nodes.begin();
+    while(it!=nodes.end()){
+        XMLNode subnode = *it;
+        if(subnode.getName()=="Option"){
+            PyObject *obj;
+            if(subnode.hasAttribute("value")){
+                obj = stringToPythonValue(subnode.getAttribute("value"));
+            }else{
+                obj = stringToPythonValue(subnode.getContent());
+            }
+            PyDict_SetItemString(opts, subnode.getAttribute("key").c_str(), obj);
+            Py_DECREF(obj);
+        }else{
+            PyObject *obj = stringToPythonValue(subnode.getContent());
+            PyDict_SetItemString(dct, subnode.getName().c_str(), obj);
+            Py_DECREF(obj);
+        }
+        ++it;
+    }
+    PyDict_SetItemString(dct, "options", opts);
+    Py_DECREF(opts);
+    return dct;
+}
+
+}
+#endif
diff --git a/src/ReconstructionAlgorithm2D.cpp b/src/ReconstructionAlgorithm2D.cpp
index 767efe6..4575ff7 100644
--- a/src/ReconstructionAlgorithm2D.cpp
+++ b/src/ReconstructionAlgorithm2D.cpp
@@ -85,9 +85,16 @@ bool CReconstructionAlgorithm2D::initialize(const Config& _cfg)
 	
 	// projector
 	XMLNode node = _cfg.self.getSingleNode("ProjectorId");
-	ASTRA_CONFIG_CHECK(node, "Reconstruction2D", "No ProjectorId tag specified.");
-	int id = boost::lexical_cast<int>(node.getContent());
-	m_pProjector = CProjector2DManager::getSingleton().get(id);
+	if (requiresProjector()) {
+		ASTRA_CONFIG_CHECK(node, "Reconstruction2D", "No ProjectorId tag specified.");
+	}
+	int id;
+	if (node) {
+		id = boost::lexical_cast<int>(node.getContent());
+		m_pProjector = CProjector2DManager::getSingleton().get(id);
+	} else {
+		m_pProjector = 0;
+	}
 	CC.markNodeParsed("ProjectorId");
 
 	// sinogram data
@@ -205,18 +212,22 @@ void CReconstructionAlgorithm2D::setSinogramMask(CFloat32ProjectionData2D* _pMas
 bool CReconstructionAlgorithm2D::_check() 
 {
 	// check pointers
-	ASTRA_CONFIG_CHECK(m_pProjector, "Reconstruction2D", "Invalid Projector Object.");
+	if (requiresProjector())
+		ASTRA_CONFIG_CHECK(m_pProjector, "Reconstruction2D", "Invalid Projector Object.");
 	ASTRA_CONFIG_CHECK(m_pSinogram, "Reconstruction2D", "Invalid Projection Data Object.");
 	ASTRA_CONFIG_CHECK(m_pReconstruction, "Reconstruction2D", "Invalid Reconstruction Data Object.");
 
 	// check initializations
-	ASTRA_CONFIG_CHECK(m_pProjector->isInitialized(), "Reconstruction2D", "Projector Object Not Initialized.");
+	if (requiresProjector())
+		ASTRA_CONFIG_CHECK(m_pProjector->isInitialized(), "Reconstruction2D", "Projector Object Not Initialized.");
 	ASTRA_CONFIG_CHECK(m_pSinogram->isInitialized(), "Reconstruction2D", "Projection Data Object Not Initialized.");
 	ASTRA_CONFIG_CHECK(m_pReconstruction->isInitialized(), "Reconstruction2D", "Reconstruction Data Object Not Initialized.");
 
 	// check compatibility between projector and data classes
-	ASTRA_CONFIG_CHECK(m_pSinogram->getGeometry()->isEqual(m_pProjector->getProjectionGeometry()), "Reconstruction2D", "Projection Data not compatible with the specified Projector.");
-	ASTRA_CONFIG_CHECK(m_pReconstruction->getGeometry()->isEqual(m_pProjector->getVolumeGeometry()), "Reconstruction2D", "Reconstruction Data not compatible with the specified Projector.");
+	if (requiresProjector()) {
+		ASTRA_CONFIG_CHECK(m_pSinogram->getGeometry()->isEqual(m_pProjector->getProjectionGeometry()), "Reconstruction2D", "Projection Data not compatible with the specified Projector.");
+		ASTRA_CONFIG_CHECK(m_pReconstruction->getGeometry()->isEqual(m_pProjector->getVolumeGeometry()), "Reconstruction2D", "Reconstruction Data not compatible with the specified Projector.");
+	}
 
 	// success
 	return true;
diff --git a/src/VolumeGeometry3D.cpp b/src/VolumeGeometry3D.cpp
index a1cf424..3de146f 100644
--- a/src/VolumeGeometry3D.cpp
+++ b/src/VolumeGeometry3D.cpp
@@ -45,6 +45,7 @@ bool CVolumeGeometry3D::_check()
 	ASTRA_CONFIG_CHECK(m_fWindowMinZ < m_fWindowMaxZ, "VolumeGeometry3D", "WindowMinZ should be lower than WindowMaxZ.");
 
 	ASTRA_CONFIG_CHECK(m_iGridTotCount == (m_iGridColCount * m_iGridRowCount * m_iGridSliceCount), "VolumeGeometry3D", "Internal configuration error.");
+#if 0
 	ASTRA_CONFIG_CHECK(m_fWindowLengthX == (m_fWindowMaxX - m_fWindowMinX), "VolumeGeometry3D", "Internal configuration error.");
 	ASTRA_CONFIG_CHECK(m_fWindowLengthY == (m_fWindowMaxY - m_fWindowMinY), "VolumeGeometry3D", "Internal configuration error.");
 	ASTRA_CONFIG_CHECK(m_fWindowLengthZ == (m_fWindowMaxZ - m_fWindowMinZ), "VolumeGeometry3D", "Internal configuration error.");
@@ -57,6 +58,7 @@ bool CVolumeGeometry3D::_check()
 	ASTRA_CONFIG_CHECK(m_fDivPixelLengthX == (1.0f / m_fPixelLengthX), "VolumeGeometry3D", "Internal configuration error.");
 	ASTRA_CONFIG_CHECK(m_fDivPixelLengthY == (1.0f / m_fPixelLengthY), "VolumeGeometry3D", "Internal configuration error.");
 	ASTRA_CONFIG_CHECK(m_fDivPixelLengthZ == (1.0f / m_fPixelLengthZ), "VolumeGeometry3D", "Internal configuration error.");
+#endif
 
 	return true;
 }