Merge pull request #101 from wjp/composite

 Add CompositeGeometryManager

28 files changed, 2306 insertions, 107 deletions

diff --git a/astra_vc09.vcproj b/astra_vc09.vcproj
index e5d7731..b928662 100644
--- a/astra_vc09.vcproj
+++ b/astra_vc09.vcproj
@@ -933,6 +933,10 @@
 					>
 				</File>
 				<File
+					RelativePath=".\include\astra\CompositeGeometryManager.h"
+					>
+				</File>
+				<File
 					RelativePath=".\include\astra\Config.h"
 					>
 				</File>
@@ -989,6 +993,10 @@
 					>
 				</File>
 				<File
+					RelativePath=".\src\CompositeGeometryManager.cpp"
+					>
+				</File>
+				<File
 					RelativePath=".\src\Config.cpp"
 					>
 				</File>
@@ -2229,6 +2237,10 @@
 					>
 				</File>
 				<File
+					RelativePath=".\cuda\3d\mem3d.h"
+					>
+				</File>
+				<File
 					RelativePath=".\cuda\3d\par3d_bp.h"
 					>
 				</File>
@@ -3041,6 +3053,42 @@
 					</FileConfiguration>
 				</File>
 				<File
+					RelativePath=".\cuda\3d\mem3d.cu"
+					>
+					<FileConfiguration
+						Name="Debug|Win32"
+						ExcludedFromBuild="true"
+						>
+						<Tool
+							Name="Cudart Build Rule"
+						/>
+					</FileConfiguration>
+					<FileConfiguration
+						Name="Debug|x64"
+						ExcludedFromBuild="true"
+						>
+						<Tool
+							Name="Cudart Build Rule"
+						/>
+					</FileConfiguration>
+					<FileConfiguration
+						Name="Release|Win32"
+						ExcludedFromBuild="true"
+						>
+						<Tool
+							Name="Cudart Build Rule"
+						/>
+					</FileConfiguration>
+					<FileConfiguration
+						Name="Release|x64"
+						ExcludedFromBuild="true"
+						>
+						<Tool
+							Name="Cudart Build Rule"
+						/>
+					</FileConfiguration>
+				</File>
+				<File
 					RelativePath=".\cuda\3d\par3d_bp.cu"
 					>
 					<FileConfiguration
diff --git a/astra_vc11.vcxproj b/astra_vc11.vcxproj
index bc11b23..fc8b9ce 100644
--- a/astra_vc11.vcxproj
+++ b/astra_vc11.vcxproj
@@ -380,6 +380,7 @@
     <ClCompile Include="src\AsyncAlgorithm.cpp" />
     <ClCompile Include="src\BackProjectionAlgorithm.cpp" />
     <ClCompile Include="src\CglsAlgorithm.cpp" />
+    <ClCompile Include="src\CompositeGeometryManager.cpp" />
     <ClCompile Include="src\ConeProjectionGeometry3D.cpp" />
     <ClCompile Include="src\ConeVecProjectionGeometry3D.cpp" />
     <ClCompile Include="src\Config.cpp" />
@@ -582,6 +583,7 @@
     <ClInclude Include="cuda\3d\darthelper3d.h" />
     <ClInclude Include="cuda\3d\dims3d.h" />
     <ClInclude Include="cuda\3d\fdk.h" />
+    <ClInclude Include="cuda\3d\mem3d.h" />
     <ClInclude Include="cuda\3d\par3d_bp.h" />
     <ClInclude Include="cuda\3d\par3d_fp.h" />
     <ClInclude Include="cuda\3d\sirt3d.h" />
@@ -594,6 +596,7 @@
     <ClInclude Include="include\astra\AsyncAlgorithm.h" />
     <ClInclude Include="include\astra\BackProjectionAlgorithm.h" />
     <ClInclude Include="include\astra\CglsAlgorithm.h" />
+    <ClInclude Include="include\astra\CompositeGeometryManager.h" />
     <ClInclude Include="include\astra\ConeProjectionGeometry3D.h" />
     <ClInclude Include="include\astra\ConeVecProjectionGeometry3D.h" />
     <ClInclude Include="include\astra\Config.h" />
@@ -804,6 +807,12 @@
       <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">true</ExcludedFromBuild>
       <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Release|x64'">true</ExcludedFromBuild>
     </CudaCompile>
+    <CudaCompile Include="cuda\3d\mem3d.cu">
+      <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">true</ExcludedFromBuild>
+      <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">true</ExcludedFromBuild>
+      <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Release|Win32'">true</ExcludedFromBuild>
+      <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Release|x64'">true</ExcludedFromBuild>
+    </CudaCompile>
     <CudaCompile Include="cuda\3d\par3d_bp.cu">
       <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Debug|Win32'">true</ExcludedFromBuild>
       <ExcludedFromBuild Condition="'$(Configuration)|$(Platform)'=='Debug|x64'">true</ExcludedFromBuild>
diff --git a/astra_vc11.vcxproj.filters b/astra_vc11.vcxproj.filters
index a597962..af8ca39 100644
--- a/astra_vc11.vcxproj.filters
+++ b/astra_vc11.vcxproj.filters
@@ -67,6 +67,9 @@
     <CudaCompile Include="cuda\3d\fdk.cu">
       <Filter>CUDA\cuda source</Filter>
     </CudaCompile>
+    <CudaCompile Include="cuda\3d\mem3d.cu">
+      <Filter>CUDA\cuda source</Filter>
+    </CudaCompile>
     <CudaCompile Include="cuda\3d\par3d_bp.cu">
       <Filter>CUDA\cuda source</Filter>
     </CudaCompile>
@@ -153,6 +156,9 @@
     <ClCompile Include="src\AstraObjectManager.cpp">
       <Filter>Global &amp; Other\source</Filter>
     </ClCompile>
+    <ClCompile Include="src\CompositeGeometryManager.cpp">
+      <Filter>Global &amp; Other\source</Filter>
+    </ClCompile>
     <ClCompile Include="src\Config.cpp">
       <Filter>Global &amp; Other\source</Filter>
     </ClCompile>
@@ -398,6 +404,9 @@
     <ClInclude Include="include\astra\clog.h">
       <Filter>Global &amp; Other\headers</Filter>
     </ClInclude>
+    <ClInclude Include="include\astra\CompositeGeometryManager.h">
+      <Filter>Global &amp; Other\headers</Filter>
+    </ClInclude>
     <ClInclude Include="include\astra\Config.h">
       <Filter>Global &amp; Other\headers</Filter>
     </ClInclude>
@@ -641,6 +650,9 @@
     <ClInclude Include="cuda\3d\fdk.h">
       <Filter>CUDA\cuda headers</Filter>
     </ClInclude>
+    <ClInclude Include="cuda\3d\mem3d.h">
+      <Filter>CUDA\cuda headers</Filter>
+    </ClInclude>
     <ClInclude Include="cuda\3d\par3d_bp.h">
       <Filter>CUDA\cuda headers</Filter>
     </ClInclude>
diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in
index bdffd4c..0361e97 100644
--- a/build/linux/Makefile.in
+++ b/build/linux/Makefile.in
@@ -105,6 +105,7 @@ BASE_OBJECTS=\
 	src/AstraObjectManager.lo \
 	src/BackProjectionAlgorithm.lo \
 	src/CglsAlgorithm.lo \
+	src/CompositeGeometryManager.lo \
 	src/ConeProjectionGeometry3D.lo \
 	src/ConeVecProjectionGeometry3D.lo \
 	src/Config.lo \
@@ -203,7 +204,8 @@ CUDA_OBJECTS=\
 	cuda/3d/sirt3d.lo \
 	cuda/3d/astra3d.lo \
 	cuda/3d/util3d.lo \
-	cuda/3d/arith3d.lo
+	cuda/3d/arith3d.lo \
+	cuda/3d/mem3d.lo
 
 ALL_OBJECTS=$(BASE_OBJECTS)
 ifeq ($(cuda),yes)
diff --git a/build/msvc/gen.py b/build/msvc/gen.py
index 72d4582..c18c1e8 100644
--- a/build/msvc/gen.py
+++ b/build/msvc/gen.py
@@ -168,6 +168,7 @@ P_astra["filters"]["CUDA\\cuda source"] = [
 "cuda\\3d\\cone_fp.cu",
 "cuda\\3d\\darthelper3d.cu",
 "cuda\\3d\\fdk.cu",
+"cuda\\3d\\mem3d.cu",
 "cuda\\3d\\par3d_bp.cu",
 "cuda\\3d\\par3d_fp.cu",
 "cuda\\3d\\sirt3d.cu",
@@ -205,6 +206,7 @@ P_astra["filters"]["Global &amp; Other\\source"] = [
 "1546cb47-7e5b-42c2-b695-ef172024c14b",
 "src\\AstraObjectFactory.cpp",
 "src\\AstraObjectManager.cpp",
+"src\\CompositeGeometryManager.cpp",
 "src\\Config.cpp",
 "src\\Fourier.cpp",
 "src\\Globals.cpp",
@@ -295,6 +297,7 @@ P_astra["filters"]["CUDA\\cuda headers"] = [
 "cuda\\3d\\darthelper3d.h",
 "cuda\\3d\\dims3d.h",
 "cuda\\3d\\fdk.h",
+"cuda\\3d\\mem3d.h",
 "cuda\\3d\\par3d_bp.h",
 "cuda\\3d\\par3d_fp.h",
 "cuda\\3d\\sirt3d.h",
@@ -336,6 +339,7 @@ P_astra["filters"]["Global &amp; Other\\headers"] = [
 "include\\astra\\AstraObjectFactory.h",
 "include\\astra\\AstraObjectManager.h",
 "include\\astra\\clog.h",
+"include\\astra\\CompositeGeometryManager.h",
 "include\\astra\\Config.h",
 "include\\astra\\Fourier.h",
 "include\\astra\\Globals.h",
diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu
index 3815a1a..8328229 100644
--- a/cuda/3d/astra3d.cu
+++ b/cuda/3d/astra3d.cu
@@ -58,88 +58,6 @@ enum CUDAProjectionType3d {
 };
 
 
-static 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;
-}
-
-static 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;
-}
-
 
 
 
diff --git a/cuda/3d/astra3d.h b/cuda/3d/astra3d.h
index 6c3fcfb..2782994 100644
--- a/cuda/3d/astra3d.h
+++ b/cuda/3d/astra3d.h
@@ -281,6 +281,15 @@ protected:
 	AstraCGLS3d_internal *pData;
 };
 
+bool convertAstraGeometry_dims(const CVolumeGeometry3D* pVolGeom,
+                               const CProjectionGeometry3D* pProjGeom,
+                               astraCUDA3d::SDimensions3D& dims);
+
+bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
+                          const CProjectionGeometry3D* pProjGeom,
+                          SPar3DProjection*& pParProjs,
+                          SConeProjection*& pConeProjs,
+                          float& fOutputScale);
 
 _AstraExport bool astraCudaFP(const float* pfVolume, float* pfProjections,
                       const CVolumeGeometry3D* pVolGeom,
diff --git a/cuda/3d/mem3d.cu b/cuda/3d/mem3d.cu
new file mode 100644
index 0000000..6d81dc0
--- /dev/null
+++ b/cuda/3d/mem3d.cu
@@ -0,0 +1,270 @@
+/*
+-----------------------------------------------------------------------
+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$
+*/
+
+#include <cstdio>
+#include <cassert>
+
+#include "util3d.h"
+
+#include "mem3d.h"
+
+#include "astra3d.h"
+#include "cone_fp.h"
+#include "cone_bp.h"
+#include "par3d_fp.h"
+#include "par3d_bp.h"
+
+#include "astra/Logging.h"
+
+
+namespace astraCUDA3d {
+
+
+struct SMemHandle3D_internal
+{
+	cudaPitchedPtr ptr;
+	unsigned int nx;
+	unsigned int ny;
+	unsigned int nz;
+};
+
+size_t availableGPUMemory()
+{
+	size_t free, total;
+	cudaError_t err = cudaMemGetInfo(&free, &total);
+	if (err != cudaSuccess)
+		return 0;
+	return free;
+}
+
+MemHandle3D allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, Mem3DZeroMode zero)
+{
+	SMemHandle3D_internal hnd;
+	hnd.nx = x;
+	hnd.ny = y;
+	hnd.nz = z;
+
+	size_t free = availableGPUMemory();
+
+	cudaError_t err;
+	err = cudaMalloc3D(&hnd.ptr, make_cudaExtent(sizeof(float)*x, y, z));
+
+	if (err != cudaSuccess) {
+		return MemHandle3D();
+	}
+
+	size_t free2 = availableGPUMemory();
+
+	ASTRA_DEBUG("Allocated %d x %d x %d on GPU. (Pre: %lu, post: %lu)", x, y, z, free, free2);
+
+
+
+	if (zero == INIT_ZERO) {
+		err = cudaMemset3D(hnd.ptr, 0, make_cudaExtent(sizeof(float)*x, y, z));
+		if (err != cudaSuccess) {
+			cudaFree(hnd.ptr.ptr);
+			return MemHandle3D();
+		}
+	}
+
+	MemHandle3D ret;
+	ret.d = boost::shared_ptr<SMemHandle3D_internal>(new SMemHandle3D_internal);
+	*ret.d = hnd;
+
+	return ret;
+}
+
+bool freeGPUMemory(MemHandle3D handle)
+{
+	size_t free = availableGPUMemory();
+	cudaError_t err = cudaFree(handle.d->ptr.ptr);
+	size_t free2 = availableGPUMemory();
+
+	ASTRA_DEBUG("Freeing memory. (Pre: %lu, post: %lu)", free, free2);
+
+	return err == cudaSuccess;
+}
+
+bool copyToGPUMemory(const float *src, MemHandle3D dst, const SSubDimensions3D &pos)
+{
+	ASTRA_DEBUG("Copying %d x %d x %d to GPU", pos.subnx, pos.subny, pos.subnz);
+	ASTRA_DEBUG("Offset %d,%d,%d", pos.subx, pos.suby, pos.subz);
+	cudaPitchedPtr s;
+	s.ptr = (void*)src; // const cast away
+	s.pitch = pos.pitch * sizeof(float);
+	s.xsize = pos.nx * sizeof(float);
+	s.ysize = pos.ny;
+	ASTRA_DEBUG("Pitch %d, xsize %d, ysize %d", s.pitch, s.xsize, s.ysize);
+
+	cudaMemcpy3DParms p;
+	p.srcArray = 0;
+	p.srcPos = make_cudaPos(pos.subx * sizeof(float), pos.suby, pos.subz);
+	p.srcPtr = s;
+
+	p.dstArray = 0;
+	p.dstPos = make_cudaPos(0, 0, 0);
+	p.dstPtr = dst.d->ptr;
+
+	p.extent = make_cudaExtent(pos.subnx * sizeof(float), pos.subny, pos.subnz);
+
+	p.kind = cudaMemcpyHostToDevice;
+
+	cudaError_t err = cudaMemcpy3D(&p);
+
+	return err == cudaSuccess;
+}
+
+
+bool copyFromGPUMemory(float *dst, MemHandle3D src, const SSubDimensions3D &pos)
+{
+	ASTRA_DEBUG("Copying %d x %d x %d from GPU", pos.subnx, pos.subny, pos.subnz);
+	ASTRA_DEBUG("Offset %d,%d,%d", pos.subx, pos.suby, pos.subz);
+	cudaPitchedPtr d;
+	d.ptr = (void*)dst;
+	d.pitch = pos.pitch * sizeof(float);
+	d.xsize = pos.nx * sizeof(float);
+	d.ysize = pos.ny;
+	ASTRA_DEBUG("Pitch %d, xsize %d, ysize %d", d.pitch, d.xsize, d.ysize);
+
+	cudaMemcpy3DParms p;
+	p.srcArray = 0;
+	p.srcPos = make_cudaPos(0, 0, 0);
+	p.srcPtr = src.d->ptr;
+
+	p.dstArray = 0;
+	p.dstPos = make_cudaPos(pos.subx * sizeof(float), pos.suby, pos.subz);
+	p.dstPtr = d;
+
+	p.extent = make_cudaExtent(pos.subnx * sizeof(float), pos.subny, pos.subnz);
+
+	p.kind = cudaMemcpyDeviceToHost;
+
+	cudaError_t err = cudaMemcpy3D(&p);
+
+	return err == cudaSuccess;
+
+}
+
+
+bool FP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, int iDetectorSuperSampling, astra::Cuda3DProjectionKernel projKernel)
+{
+	SDimensions3D dims;
+
+	bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims);
+	if (!ok)
+		return false;
+
+#if 1
+	dims.iRaysPerDetDim = iDetectorSuperSampling;
+	if (iDetectorSuperSampling == 0)
+		return false;
+#else
+	dims.iRaysPerDetDim = 1;
+	astra::Cuda3DProjectionKernel projKernel = astra::ker3d_default;
+#endif
+
+
+	SPar3DProjection* pParProjs;
+	SConeProjection* pConeProjs;
+
+	float outputScale = 1.0f;
+
+	ok = convertAstraGeometry(pVolGeom, pProjGeom,
+	                          pParProjs, pConeProjs,
+	                          outputScale);
+
+	if (pParProjs) {
+#if 0
+		for (int i = 0; i < dims.iProjAngles; ++i) {
+			ASTRA_DEBUG("Vec: %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\n",
+			    pParProjs[i].fRayX, pParProjs[i].fRayY, pParProjs[i].fRayZ,
+			    pParProjs[i].fDetSX, pParProjs[i].fDetSY, pParProjs[i].fDetSZ,
+			    pParProjs[i].fDetUX, pParProjs[i].fDetUY, pParProjs[i].fDetUZ,
+			    pParProjs[i].fDetVX, pParProjs[i].fDetVY, pParProjs[i].fDetVZ);
+		}
+#endif
+
+		switch (projKernel) {
+		case astra::ker3d_default:
+			ok &= Par3DFP(volData.d->ptr, projData.d->ptr, dims, pParProjs, outputScale);
+			break;
+		case astra::ker3d_sum_square_weights:
+			ok &= Par3DFP_SumSqW(volData.d->ptr, projData.d->ptr, dims, pParProjs, outputScale*outputScale);
+			break;
+		default:
+			ok = false;
+		}
+	} else {
+		switch (projKernel) {
+		case astra::ker3d_default:
+			ok &= ConeFP(volData.d->ptr, projData.d->ptr, dims, pConeProjs, outputScale);
+			break;
+		default:
+			ok = false;
+		}
+	}
+
+	return ok;
+}
+
+bool BP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, int iVoxelSuperSampling)
+{
+	SDimensions3D dims;
+
+	bool ok = convertAstraGeometry_dims(pVolGeom, pProjGeom, dims);
+	if (!ok)
+		return false;
+
+#if 1
+	dims.iRaysPerVoxelDim = iVoxelSuperSampling;
+#else
+	dims.iRaysPerVoxelDim = 1;
+#endif
+
+	SPar3DProjection* pParProjs;
+	SConeProjection* pConeProjs;
+
+	float outputScale = 1.0f;
+
+	ok = convertAstraGeometry(pVolGeom, pProjGeom,
+	                          pParProjs, pConeProjs,
+	                          outputScale);
+
+	if (pParProjs)
+		ok &= Par3DBP(volData.d->ptr, projData.d->ptr, dims, pParProjs, outputScale);
+	else
+		ok &= ConeBP(volData.d->ptr, projData.d->ptr, dims, pConeProjs, outputScale);
+
+	return ok;
+
+}
+
+
+
+
+}
diff --git a/cuda/3d/mem3d.h b/cuda/3d/mem3d.h
new file mode 100644
index 0000000..82bad19
--- /dev/null
+++ b/cuda/3d/mem3d.h
@@ -0,0 +1,99 @@
+/*
+-----------------------------------------------------------------------
+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/>.
+
+-----------------------------------------------------------------------
+*/
+
+#ifndef _CUDA_MEM3D_H
+#define _CUDA_MEM3D_H
+
+#include <boost/shared_ptr.hpp>
+
+#include "astra3d.h"
+
+namespace astra {
+class CVolumeGeometry3D;
+class CProjectionGeometry3D;	
+}
+
+namespace astraCUDA3d {
+
+// TODO: Make it possible to delete these handles when they're no longer
+// necessary inside the FP/BP
+//
+// TODO: Add functions for querying capacity
+
+struct SMemHandle3D_internal;
+
+struct MemHandle3D {
+	boost::shared_ptr<SMemHandle3D_internal> d;
+	operator bool() const { return (bool)d; }
+};
+
+struct SSubDimensions3D {
+	unsigned int nx;
+	unsigned int ny;
+	unsigned int nz;
+	unsigned int pitch;
+	unsigned int subnx;
+	unsigned int subny;
+	unsigned int subnz;
+	unsigned int subx;
+	unsigned int suby;
+	unsigned int subz;
+};
+
+/*
+// Useful or not?
+enum Mem3DCopyMode {
+	MODE_SET,
+	MODE_ADD
+};
+*/
+
+enum Mem3DZeroMode {
+	INIT_NO,
+	INIT_ZERO
+};
+
+size_t availableGPUMemory();
+
+MemHandle3D allocateGPUMemory(unsigned int x, unsigned int y, unsigned int z, Mem3DZeroMode zero);
+
+bool copyToGPUMemory(const float *src, MemHandle3D dst, const SSubDimensions3D &pos);
+
+bool copyFromGPUMemory(float *dst, MemHandle3D src, const SSubDimensions3D &pos);
+
+bool freeGPUMemory(MemHandle3D handle);
+
+
+bool FP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, int iDetectorSuperSampling, astra::Cuda3DProjectionKernel projKernel);
+
+bool BP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, int iVoxelSuperSampling);
+
+
+
+}
+
+#endif
diff --git a/include/astra/CompositeGeometryManager.h b/include/astra/CompositeGeometryManager.h
new file mode 100644
index 0000000..6610151
--- /dev/null
+++ b/include/astra/CompositeGeometryManager.h
@@ -0,0 +1,152 @@
+/*
+-----------------------------------------------------------------------
+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/>.
+
+-----------------------------------------------------------------------
+*/
+
+#ifndef _INC_ASTRA_COMPOSITEGEOMETRYMANAGER
+#define _INC_ASTRA_COMPOSITEGEOMETRYMANAGER
+
+#include "Globals.h"
+
+#ifdef ASTRA_CUDA
+
+#include <list>
+#include <map>
+#include <vector>
+#include <boost/shared_ptr.hpp>
+
+
+namespace astra {
+
+class CCompositeVolume;
+class CCompositeProjections;
+class CFloat32Data3DMemory;
+class CFloat32ProjectionData3DMemory;
+class CFloat32VolumeData3DMemory;
+class CVolumeGeometry3D;
+class CProjectionGeometry3D;
+class CProjector3D;
+
+
+
+class _AstraExport CCompositeGeometryManager {
+public:
+	class CPart;
+	typedef std::list<boost::shared_ptr<CPart> > TPartList;
+	class CPart {
+	public:
+		CPart() { }
+		CPart(const CPart& other);
+		virtual ~CPart() { }
+
+		enum {
+			PART_VOL, PART_PROJ
+		} eType;
+
+		CFloat32Data3DMemory* pData;
+		unsigned int subX;
+		unsigned int subY;
+		unsigned int subZ;
+
+		bool uploadToGPU();
+		bool downloadFromGPU(/*mode?*/);
+		virtual TPartList split(size_t maxSize, int div) = 0;
+		virtual CPart* reduce(const CPart *other) = 0;
+		virtual void getDims(size_t &x, size_t &y, size_t &z) = 0;
+		size_t getSize();
+	};
+
+	class CVolumePart : public CPart {
+	public:
+		CVolumePart() { eType = PART_VOL; }
+		CVolumePart(const CVolumePart& other);
+		virtual ~CVolumePart();
+
+		CVolumeGeometry3D* pGeom;
+
+		virtual TPartList split(size_t maxSize, int div);
+		virtual CPart* reduce(const CPart *other);
+		virtual void getDims(size_t &x, size_t &y, size_t &z);
+
+		CVolumePart* clone() const;
+	};
+	class CProjectionPart : public CPart {
+	public:
+		CProjectionPart() { eType = PART_PROJ; }
+		CProjectionPart(const CProjectionPart& other);
+		virtual ~CProjectionPart();
+
+		CProjectionGeometry3D* pGeom;
+
+		virtual TPartList split(size_t maxSize, int div);
+		virtual CPart* reduce(const CPart *other);
+		virtual void getDims(size_t &x, size_t &y, size_t &z);
+
+		CProjectionPart* clone() const;
+	};
+
+	struct SJob {
+	public:
+		boost::shared_ptr<CPart> pInput;
+		boost::shared_ptr<CPart> pOutput;
+		CProjector3D *pProjector; // For a `global' geometry. It will not match
+		                          // the geometries of the input and output.
+
+
+		enum {
+			JOB_FP, JOB_BP, JOB_NOP
+		} eType;
+		enum {
+			MODE_ADD, MODE_SET
+		} eMode;
+
+	};
+
+	typedef std::list<SJob> TJobList;
+	// output part -> list of jobs for that output
+	typedef std::map<CPart*, TJobList > TJobSet;
+
+	bool doJobs(TJobList &jobs);
+
+	// Convenience functions for creating and running a single FP or BP job
+	bool doFP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData,
+	          CFloat32ProjectionData3DMemory *pProjData);
+	bool doBP(CProjector3D *pProjector, CFloat32VolumeData3DMemory *pVolData,
+	          CFloat32ProjectionData3DMemory *pProjData);
+
+	bool doFP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData);
+	bool doBP(CProjector3D *pProjector, const std::vector<CFloat32VolumeData3DMemory *>& volData, const std::vector<CFloat32ProjectionData3DMemory *>& projData);
+
+protected:
+
+	bool splitJobs(TJobSet &jobs, size_t maxSize, int div, TJobSet &split);
+
+};
+
+}
+
+#endif
+
+#endif
diff --git a/include/astra/ConeProjectionGeometry3D.h b/include/astra/ConeProjectionGeometry3D.h
index 00e72ce..dede6e1 100644
--- a/include/astra/ConeProjectionGeometry3D.h
+++ b/include/astra/ConeProjectionGeometry3D.h
@@ -186,9 +186,15 @@ public:
 	 */
 	virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex) const;
 
-	virtual void projectPoint(float32 fX, float32 fY, float32 fZ,
+	virtual void projectPoint(double fX, double fY, double fZ,
 	                          int iAngleIndex,
-	                          float32 &fU, float32 &fV) const;
+	                          double &fU, double &fV) const;
+	virtual void backprojectPointX(int iAngleIndex, double fU, double fV,
+	                               double fX, double &fY, double &fZ) const;
+	virtual void backprojectPointY(int iAngleIndex, double fU, double fV,
+	                               double fY, double &fX, double &fZ) const;
+	virtual void backprojectPointZ(int iAngleIndex, double fU, double fV,
+	                               double fZ, double &fX, double &fY) const;
 
 };
 
diff --git a/include/astra/ConeVecProjectionGeometry3D.h b/include/astra/ConeVecProjectionGeometry3D.h
index 71e8010..f76f9dd 100644
--- a/include/astra/ConeVecProjectionGeometry3D.h
+++ b/include/astra/ConeVecProjectionGeometry3D.h
@@ -148,9 +148,16 @@ public:
 
 	const SConeProjection* getProjectionVectors() const { return m_pProjectionAngles; }
 
-	virtual void projectPoint(float32 fX, float32 fY, float32 fZ,
+	virtual void projectPoint(double fX, double fY, double fZ,
 	                          int iAngleIndex,
-	                          float32 &fU, float32 &fV) const;
+	                          double &fU, double &fV) const;
+	virtual void backprojectPointX(int iAngleIndex, double fU, double fV,
+	                               double fX, double &fY, double &fZ) const;
+	virtual void backprojectPointY(int iAngleIndex, double fU, double fV,
+	                               double fY, double &fX, double &fZ) const;
+	virtual void backprojectPointZ(int iAngleIndex, double fU, double fV,
+	                               double fZ, double &fX, double &fY) const;
+
 };
 
 } // namespace astra
diff --git a/include/astra/GeometryUtil3D.h b/include/astra/GeometryUtil3D.h
index 6ceac63..e4d73e4 100644
--- a/include/astra/GeometryUtil3D.h
+++ b/include/astra/GeometryUtil3D.h
@@ -119,6 +119,23 @@ void computeBP_UV_Coeffs(const SConeProjection& proj,
                          double &fDX, double &fDY, double &fDZ, double &fDC);
 
 
+SConeProjection* genConeProjections(unsigned int iProjAngles,
+                                    unsigned int iProjU,
+                                    unsigned int iProjV,
+                                    double fOriginSourceDistance,
+                                    double fOriginDetectorDistance,
+                                    double fDetUSize,
+                                    double fDetVSize,
+                                    const float *pfAngles);
+
+SPar3DProjection* genPar3DProjections(unsigned int iProjAngles,
+                                      unsigned int iProjU,
+                                      unsigned int iProjV,
+                                      double fDetUSize,
+                                      double fDetVSize,
+                                      const float *pfAngles);
+
+
 
 
 }
diff --git a/include/astra/ParallelProjectionGeometry3D.h b/include/astra/ParallelProjectionGeometry3D.h
index 72401e5..d95c050 100644
--- a/include/astra/ParallelProjectionGeometry3D.h
+++ b/include/astra/ParallelProjectionGeometry3D.h
@@ -147,9 +147,16 @@ public:
 	  */
 	virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex) const;
 
-	virtual void projectPoint(float32 fX, float32 fY, float32 fZ,
+	virtual void projectPoint(double fX, double fY, double fZ,
 	                          int iAngleIndex,
-	                          float32 &fU, float32 &fV) const;
+	                          double &fU, double &fV) const;
+	virtual void backprojectPointX(int iAngleIndex, double fU, double fV,
+	                               double fX, double &fY, double &fZ) const;
+	virtual void backprojectPointY(int iAngleIndex, double fU, double fV,
+	                               double fY, double &fX, double &fZ) const;
+	virtual void backprojectPointZ(int iAngleIndex, double fU, double fV,
+	                               double fZ, double &fX, double &fY) const;
+
 
 	 /**
 	  * Creates (= allocates) a 2D projection geometry used when projecting one slice using a 2D projector
diff --git a/include/astra/ParallelVecProjectionGeometry3D.h b/include/astra/ParallelVecProjectionGeometry3D.h
index 59238c8..ec91086 100644
--- a/include/astra/ParallelVecProjectionGeometry3D.h
+++ b/include/astra/ParallelVecProjectionGeometry3D.h
@@ -149,9 +149,15 @@ public:
 
 	const SPar3DProjection* getProjectionVectors() const { return m_pProjectionAngles; }
 
-	virtual void projectPoint(float32 fX, float32 fY, float32 fZ,
+	virtual void projectPoint(double fX, double fY, double fZ,
 	                          int iAngleIndex,
-	                          float32 &fU, float32 &fV) const;
+	                          double &fU, double &fV) const;
+	virtual void backprojectPointX(int iAngleIndex, double fU, double fV,
+	                               double fX, double &fY, double &fZ) const;
+	virtual void backprojectPointY(int iAngleIndex, double fU, double fV,
+	                               double fY, double &fX, double &fZ) const;
+	virtual void backprojectPointZ(int iAngleIndex, double fU, double fV,
+	                               double fZ, double &fX, double &fY) const;
 };
 
 } // namespace astra
diff --git a/include/astra/ProjectionGeometry3D.h b/include/astra/ProjectionGeometry3D.h
index 19ac3ab..0b60287 100644
--- a/include/astra/ProjectionGeometry3D.h
+++ b/include/astra/ProjectionGeometry3D.h
@@ -317,9 +317,24 @@ public:
 	 * @param iAngleIndex	the index of the angle to use
 	 * @param fU,fV		the projected point.
 	 */
-	virtual void projectPoint(float32 fX, float32 fY, float32 fZ,
+	virtual void projectPoint(double fX, double fY, double fZ,
 	                          int iAngleIndex,
-	                          float32 &fU, float32 &fV) const = 0;
+	                          double &fU, double &fV) const = 0;
+
+	/* Backproject a point onto a plane parallel to a coordinate plane.
+	 * The 2D point coordinates are the (unrounded) indices of the detector
+	 * column and row. The output is in 3D coordinates in units.
+	 * are in units. The output fU,fV are the (unrounded) indices of the
+	 * detector column and row.
+	 * This may fall outside of the actual detector.
+	 */
+	virtual void backprojectPointX(int iAngleIndex, double fU, double fV,
+	                               double fX, double &fY, double &fZ) const = 0;
+	virtual void backprojectPointY(int iAngleIndex, double fU, double fV,
+	                               double fY, double &fX, double &fZ) const = 0;
+	virtual void backprojectPointZ(int iAngleIndex, double fU, double fV,
+	                               double fZ, double &fX, double &fY) const = 0;
+
 
 	/** Returns true if the type of geometry defined in this class is the one specified in _sType.
 	 *
diff --git a/python/astra/PyIncludes.pxd b/python/astra/PyIncludes.pxd
index 77346b0..4a4ce43 100644
--- a/python/astra/PyIncludes.pxd
+++ b/python/astra/PyIncludes.pxd
@@ -226,6 +226,7 @@ cdef extern from "astra/Float32VolumeData3DMemory.h" namespace "astra":
 		int getRowCount()
 		int getColCount()
 		int getSliceCount()
+		bool isInitialized()
 
 
 
@@ -257,6 +258,7 @@ cdef extern from "astra/Float32ProjectionData3DMemory.h" namespace "astra":
 		int getDetectorColCount()
 		int getDetectorRowCount()
 		int getAngleCount()
+		bool isInitialized()
 
 cdef extern from "astra/Float32Data3D.h" namespace "astra":
 	cdef cppclass CFloat32Data3D:
diff --git a/python/astra/experimental.pyx b/python/astra/experimental.pyx
new file mode 100644
index 0000000..da27504
--- /dev/null
+++ b/python/astra/experimental.pyx
@@ -0,0 +1,84 @@
+#-----------------------------------------------------------------------
+# 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/>.
+#
+#-----------------------------------------------------------------------
+
+# distutils: language = c++
+# distutils: libraries = astra
+
+include "config.pxi"
+
+import six
+from .PyIncludes cimport *
+from libcpp.vector cimport vector
+
+cdef extern from "astra/CompositeGeometryManager.h" namespace "astra":
+    cdef cppclass CCompositeGeometryManager:
+        bool doFP(CProjector3D *, vector[CFloat32VolumeData3DMemory *], vector[CFloat32ProjectionData3DMemory *])
+        bool doBP(CProjector3D *, vector[CFloat32VolumeData3DMemory *], vector[CFloat32ProjectionData3DMemory *])
+
+cdef extern from *:
+    CFloat32VolumeData3DMemory * dynamic_cast_vol_mem "dynamic_cast<astra::CFloat32VolumeData3DMemory*>" (CFloat32Data3D * ) except NULL
+    CFloat32ProjectionData3DMemory * dynamic_cast_proj_mem "dynamic_cast<astra::CFloat32ProjectionData3DMemory*>" (CFloat32Data3D * ) except NULL
+
+cimport PyProjector3DManager
+from .PyProjector3DManager cimport CProjector3DManager
+cimport PyData3DManager
+from .PyData3DManager cimport CData3DManager
+
+cdef CProjector3DManager * manProj = <CProjector3DManager * >PyProjector3DManager.getSingletonPtr()
+cdef CData3DManager * man3d = <CData3DManager * >PyData3DManager.getSingletonPtr()
+
+def do_composite(projector_id, vol_ids, proj_ids, t):
+    cdef vector[CFloat32VolumeData3DMemory *] vol
+    cdef CFloat32VolumeData3DMemory * pVolObject
+    cdef CFloat32ProjectionData3DMemory * pProjObject
+    for v in vol_ids:
+        pVolObject = dynamic_cast_vol_mem(man3d.get(v))
+        if pVolObject == NULL:
+            raise Exception("Data object not found")
+        if not pVolObject.isInitialized():
+            raise Exception("Data object not initialized properly")
+        vol.push_back(pVolObject)
+    cdef vector[CFloat32ProjectionData3DMemory *] proj
+    for v in proj_ids:
+        pProjObject = dynamic_cast_proj_mem(man3d.get(v))
+        if pProjObject == NULL:
+            raise Exception("Data object not found")
+        if not pProjObject.isInitialized():
+            raise Exception("Data object not initialized properly")
+        proj.push_back(pProjObject)
+    cdef CCompositeGeometryManager m
+    cdef CProjector3D * projector = manProj.get(projector_id) # may be NULL
+    if t == "FP":
+        if not m.doFP(projector, vol, proj):
+            raise Exception("Failed to perform FP")
+    else:
+        if not m.doBP(projector, vol, proj):
+            raise Exception("Failed to perform BP")
+
+def do_composite_FP(projector_id, vol_ids, proj_ids):
+    do_composite(projector_id, vol_ids, proj_ids, "FP")
+
+def do_composite_BP(projector_id, vol_ids, proj_ids):
+    do_composite(projector_id, vol_ids, proj_ids, "BP")
diff --git a/samples/python/s018_experimental_multires.py b/samples/python/s018_experimental_multires.py
new file mode 100644
index 0000000..cf38e53
--- /dev/null
+++ b/samples/python/s018_experimental_multires.py
@@ -0,0 +1,84 @@
+#-----------------------------------------------------------------------
+#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam
+#
+#Author: Daniel M. Pelt
+#Contact: D.M.Pelt@cwi.nl
+#Website: http://dmpelt.github.io/pyastratoolbox/
+#
+#
+#This file is part of the Python interface to the
+#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox").
+#
+#The Python interface to the ASTRA Toolbox is free software: you can redistribute it and/or modify
+#it under the terms of the GNU General Public License as published by
+#the Free Software Foundation, either version 3 of the License, or
+#(at your option) any later version.
+#
+#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful,
+#but WITHOUT ANY WARRANTY; without even the implied warranty of
+#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#GNU General Public License for more details.
+#
+#You should have received a copy of the GNU General Public License
+#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+#
+#-----------------------------------------------------------------------
+
+import astra
+import numpy as np
+from astra.experimental import do_composite_FP
+
+astra.log.setOutputScreen(astra.log.STDERR, astra.log.DEBUG)
+
+# low res part (voxels of 4x4x4)
+vol_geom1 = astra.create_vol_geom(32, 16, 32, -64, 0, -64, 64, -64, 64)
+
+# high res part (voxels of 1x1x1)
+vol_geom2 = astra.create_vol_geom(128, 64, 128, 0, 64, -64, 64, -64, 64)
+
+
+# Split the output in two parts as well, for demonstration purposes
+angles1 = np.linspace(0, np.pi/2, 90, False)
+angles2 = np.linspace(np.pi/2, np.pi, 90, False)
+proj_geom1 = astra.create_proj_geom('parallel3d', 1.0, 1.0, 128, 192, angles1)
+proj_geom2 = astra.create_proj_geom('parallel3d', 1.0, 1.0, 128, 192, angles2)
+
+# Create a simple hollow cube phantom
+cube1 = np.zeros((32,32,16))
+cube1[4:28,4:28,4:16] = 1
+
+cube2 = np.zeros((128,128,64))
+cube2[16:112,16:112,0:112] = 1
+cube2[33:97,33:97,4:28] = 0
+
+vol1 = astra.data3d.create('-vol', vol_geom1, cube1)
+vol2 = astra.data3d.create('-vol', vol_geom2, cube2)
+
+proj1 = astra.data3d.create('-proj3d', proj_geom1, 0)
+proj2 = astra.data3d.create('-proj3d', proj_geom2, 0)
+
+# The actual geometries don't matter for this composite FP/BP case
+projector = astra.create_projector('cuda3d', proj_geom1, vol_geom1)
+
+do_composite_FP(projector, [vol1, vol2], [proj1, proj2])
+
+proj_data1 = astra.data3d.get(proj1)
+proj_data2 = astra.data3d.get(proj2)
+
+# Display a single projection image
+import pylab
+pylab.gray()
+pylab.figure(1)
+pylab.imshow(proj_data1[:,0,:])
+pylab.figure(2)
+pylab.imshow(proj_data2[:,0,:])
+pylab.show()
+
+
+# Clean up. Note that GPU memory is tied up in the algorithm object,
+# and main RAM in the data objects.
+astra.data3d.delete(vol1)
+astra.data3d.delete(vol2)
+astra.data3d.delete(proj1)
+astra.data3d.delete(proj2)
+astra.projector3d.delete(projector)
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/CudaBackProjectionAlgorithm3D.cpp b/src/CudaBackProjectionAlgorithm3D.cpp
index 8cf4c3b..ce8e111 100644
--- a/src/CudaBackProjectionAlgorithm3D.cpp
+++ b/src/CudaBackProjectionAlgorithm3D.cpp
@@ -37,6 +37,7 @@ $Id$
 #include "astra/ParallelProjectionGeometry3D.h"
 #include "astra/ParallelVecProjectionGeometry3D.h"
 #include "astra/ConeVecProjectionGeometry3D.h"
+#include "astra/CompositeGeometryManager.h"
 
 #include "astra/Logging.h"
 
@@ -203,9 +204,16 @@ void CCudaBackProjectionAlgorithm3D::run(int _iNrIterations)
 		                         &volgeom, projgeom,
 		                         m_iGPUIndex, m_iVoxelSuperSampling);
 	} else {
+
+#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/CudaForwardProjectionAlgorithm3D.cpp b/src/CudaForwardProjectionAlgorithm3D.cpp
index e57e077..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"
@@ -263,6 +265,12 @@ 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 CVolumeGeometry3D& volgeom = *m_pVolume->getGeometry();
 
@@ -294,6 +302,7 @@ void CCudaForwardProjectionAlgorithm3D::run(int)
 	astraCudaFP(m_pVolume->getDataConst(), m_pProjections->getData(),
 	            &volgeom, projgeom,
 	            m_iGPUIndex, m_iDetectorSuperSampling, projKernel);
+#endif
 }
 
 
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/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/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;
 }