Add CompositeGeometryManager

This handles FP and BP operations on multiple data objects at once,
splitting them to fit in GPU memory where necessary.

24 files changed, 2023 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 abbebe2..c555bca 100644
--- a/build/linux/Makefile.in
+++ b/build/linux/Makefile.in
@@ -99,6 +99,7 @@ BASE_OBJECTS=\
 	src/AstraObjectManager.lo \
 	src/BackProjectionAlgorithm.lo \
 	src/CglsAlgorithm.lo \
+	src/CompositeGeometryManager.lo \
 	src/ConeProjectionGeometry3D.lo \
 	src/ConeVecProjectionGeometry3D.lo \
 	src/Config.lo \
@@ -197,7 +198,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..a6e57f1
--- /dev/null
+++ b/include/astra/CompositeGeometryManager.h
@@ -0,0 +1,150 @@
+/*
+-----------------------------------------------------------------------
+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);
+
+
+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/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp
new file mode 100644
index 0000000..fc8bc2e
--- /dev/null
+++ b/src/CompositeGeometryManager.cpp
@@ -0,0 +1,884 @@
+/*
+-----------------------------------------------------------------------
+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;
+				split[outputPart.get()] = TJobList();
+
+				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\n", 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::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()