Merge branch 'consistent_scaling'

48 files changed, 841 insertions, 2532 deletions

diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in
index 209206e..f478bb7 100644
--- a/build/linux/Makefile.in
+++ b/build/linux/Makefile.in
@@ -257,8 +257,6 @@ endif
 TEST_OBJECTS=\
 	tests/main.o \
 	tests/test_AstraObjectManager.o \
-	tests/test_ParallelBeamLineKernelProjector2D.o \
-	tests/test_ParallelBeamLinearKernelProjector2D.o \
 	tests/test_Float32Data2D.o \
 	tests/test_VolumeGeometry2D.o \
 	tests/test_ParallelProjectionGeometry2D.o \
diff --git a/build/linux/acinclude.m4 b/build/linux/acinclude.m4
index 92f263a..4d278c7 100644
--- a/build/linux/acinclude.m4
+++ b/build/linux/acinclude.m4
@@ -166,3 +166,20 @@ fi
 rm -f conftest.cu conftest.o conftest.nvcc.out
 ])
 
+dnl ASTRA_CHECK_CUDA_BOOST(action-if-ok, action-if-not-ok)
+dnl Check for a specific incompatibility between boost and cuda version
+dnl (See https://github.com/boostorg/config/pull/175 )
+AC_DEFUN([ASTRA_CHECK_CUDA_BOOST],[
+cat >conftest.cu <<_ACEOF
+#include <boost/shared_ptr.hpp>
+int main() {
+  return 0;
+}
+_ACEOF
+ASTRA_RUN_LOGOUTPUT([$NVCC -c -o conftest.o conftest.cu $NVCCFLAGS])
+AS_IF([test $? = 0],[$1],[
+  AS_ECHO(["$as_me: failed program was:"]) >&AS_MESSAGE_LOG_FD
+  sed 's/^/| /' conftest.cu >&AS_MESSAGE_LOG_FD
+  $2])
+])
+
diff --git a/build/linux/configure.ac b/build/linux/configure.ac
index 0a9024e..bb4d113 100644
--- a/build/linux/configure.ac
+++ b/build/linux/configure.ac
@@ -122,6 +122,14 @@ if test x"$HAVECUDA" = xyes; then
   AC_MSG_RESULT($HAVECUDA)
 fi
 
+if test x"$HAVECUDA" = xyes; then
+  AC_MSG_CHECKING([if boost and CUDA versions are compatible])
+  ASTRA_CHECK_CUDA_BOOST(AC_MSG_RESULT([yes]),[
+    AC_MSG_RESULT([no])
+    AC_MSG_ERROR([Boost and CUDA versions are incompatible. You probably have to upgrade boost.])
+  ])
+fi
+
 AC_ARG_WITH(cuda_compute, [[  --with-cuda-compute=archs  comma separated list of CUDA compute models (optional)]],,)
 if test x"$HAVECUDA" = xyes; then
   AC_MSG_CHECKING([for nvcc archs])
diff --git a/cuda/2d/algo.cu b/cuda/2d/algo.cu
index b4c2864..be15b25 100644
--- a/cuda/2d/algo.cu
+++ b/cuda/2d/algo.cu
@@ -134,8 +134,8 @@ bool ReconAlgo::setGeometry(const astra::CVolumeGeometry2D* pVolGeom,
 	delete[] fanProjs;
 	fanProjs = 0;
 
-	fOutputScale = 1.0f;
-	ok = convertAstraGeometry(pVolGeom, pProjGeom, parProjs, fanProjs, fOutputScale);
+	fProjectorScale = 1.0f;
+	ok = convertAstraGeometry(pVolGeom, pProjGeom, parProjs, fanProjs, fProjectorScale);
 	if (!ok)
 		return false;
 
@@ -242,7 +242,7 @@ bool ReconAlgo::allocateBuffers()
 	return true;
 }
 
-bool ReconAlgo::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale,
+bool ReconAlgo::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch,
                               const float* pfReconstruction, unsigned int iReconstructionPitch,
                               const float* pfVolMask, unsigned int iVolMaskPitch,
                               const float* pfSinoMask, unsigned int iSinoMaskPitch)
@@ -258,11 +258,6 @@ bool ReconAlgo::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPit
 	if (!ok)
 		return false;
 
-	// rescale sinogram to adjust for pixel size
-	processSino<opMul>(D_sinoData, fSinogramScale,
-	                       //1.0f/(fPixelSize*fPixelSize),
-	                       sinoPitch, dims);
-
 	ok = copyVolumeToDevice(pfReconstruction, iReconstructionPitch,
 	                        dims,
 	                        D_volumeData, volumePitch);
@@ -316,11 +311,11 @@ bool ReconAlgo::callFP(float* D_volumeData, unsigned int volumePitch,
 	if (parProjs) {
 		assert(!fanProjs);
 		return FP(D_volumeData, volumePitch, D_projData, projPitch,
-		          dims, parProjs, fOutputScale * outputScale);
+		          dims, parProjs, fProjectorScale * outputScale);
 	} else {
 		assert(fanProjs);
 		return FanFP(D_volumeData, volumePitch, D_projData, projPitch,
-		             dims, fanProjs, fOutputScale * outputScale);
+		             dims, fanProjs, fProjectorScale * outputScale);
 	}
 }
 
@@ -331,11 +326,11 @@ bool ReconAlgo::callBP(float* D_volumeData, unsigned int volumePitch,
 	if (parProjs) {
 		assert(!fanProjs);
 		return BP(D_volumeData, volumePitch, D_projData, projPitch,
-		          dims, parProjs, outputScale);
+		          dims, parProjs, fProjectorScale * outputScale);
 	} else {
 		assert(fanProjs);
 		return FanBP(D_volumeData, volumePitch, D_projData, projPitch,
-		             dims, fanProjs, outputScale);
+		             dims, fanProjs, fProjectorScale * outputScale);
 	}
 
 }
diff --git a/cuda/2d/astra.cu b/cuda/2d/astra.cu
index ec03517..7ff1c95 100644
--- a/cuda/2d/astra.cu
+++ b/cuda/2d/astra.cu
@@ -302,7 +302,8 @@ static bool convertAstraGeometry_internal(const CVolumeGeometry2D* pVolGeom,
 		pProjs[i].scale(factor);
 	}
 	// CHECKME: Check factor
-	fOutputScale *= pVolGeom->getPixelLengthX() * pVolGeom->getPixelLengthY();
+	// NB: Only valid for square pixels
+	fOutputScale *= pVolGeom->getPixelLengthX();
 
 	return true;
 }
diff --git a/cuda/2d/cgls.cu b/cuda/2d/cgls.cu
index b6a9fae..e7238b9 100644
--- a/cuda/2d/cgls.cu
+++ b/cuda/2d/cgls.cu
@@ -29,10 +29,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/2d/util.h"
 #include "astra/cuda/2d/arith.h"
 
-#ifdef STANDALONE
-#include "testutil.h"
-#endif
-
 #include <cstdio>
 #include <cassert>
 
@@ -102,14 +98,14 @@ bool CGLS::setBuffers(float* _D_volumeData, unsigned int _volumePitch,
 	return true;
 }
 
-bool CGLS::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale,
+bool CGLS::copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch,
                          const float* pfReconstruction, unsigned int iReconstructionPitch,
                          const float* pfVolMask, unsigned int iVolMaskPitch,
                          const float* pfSinoMask, unsigned int iSinoMaskPitch)
 {
 	sliceInitialized = false;
 
-	return ReconAlgo::copyDataToGPU(pfSinogram, iSinogramPitch, fSinogramScale, pfReconstruction, iReconstructionPitch, pfVolMask, iVolMaskPitch, pfSinoMask, iSinoMaskPitch);
+	return ReconAlgo::copyDataToGPU(pfSinogram, iSinogramPitch, pfReconstruction, iReconstructionPitch, pfVolMask, iVolMaskPitch, pfSinoMask, iSinoMaskPitch);
 }
 
 bool CGLS::iterate(unsigned int iterations)
@@ -206,60 +202,3 @@ float CGLS::computeDiffNorm()
 
 
 }
-
-#ifdef STANDALONE
-
-using namespace astraCUDA;
-
-int main()
-{
-	float* D_volumeData;
-	float* D_sinoData;
-
-	SDimensions dims;
-	dims.iVolWidth = 1024;
-	dims.iVolHeight = 1024;
-	dims.iProjAngles = 512;
-	dims.iProjDets = 1536;
-	dims.fDetScale = 1.0f;
-	dims.iRaysPerDet = 1;
-	unsigned int volumePitch, sinoPitch;
-
-	allocateVolume(D_volumeData, dims.iVolWidth, dims.iVolHeight, volumePitch);
-	zeroVolume(D_volumeData, volumePitch, dims.iVolWidth, dims.iVolHeight);
-	printf("pitch: %u\n", volumePitch);
-
-	allocateVolume(D_sinoData, dims.iProjDets, dims.iProjAngles, sinoPitch);
-	zeroVolume(D_sinoData, sinoPitch, dims.iProjDets, dims.iProjAngles);
-	printf("pitch: %u\n", sinoPitch);
-	
-	unsigned int y, x;
-	float* sino = loadImage("sino.png", y, x);
-
-	float* img = new float[dims.iVolWidth*dims.iVolHeight];
-
-	copySinogramToDevice(sino, dims.iProjDets, dims.iProjDets, dims.iProjAngles, D_sinoData, sinoPitch);
-
-	float* angle = new float[dims.iProjAngles];
-
-	for (unsigned int i = 0; i < dims.iProjAngles; ++i)
-		angle[i] = i*(M_PI/dims.iProjAngles);
-
-	CGLS cgls;
-
-	cgls.setGeometry(dims, angle);
-	cgls.init();
-
-	cgls.setBuffers(D_volumeData, volumePitch, D_sinoData, sinoPitch);
-
-	cgls.iterate(25);
-
-	delete[] angle;
-
-	copyVolumeFromDevice(img, dims.iVolWidth, dims.iVolWidth, dims.iVolHeight, D_volumeData, volumePitch);
-
-	saveImage("vol.png",dims.iVolHeight,dims.iVolWidth,img);
-
-	return 0;
-}
-#endif
diff --git a/cuda/2d/em.cu b/cuda/2d/em.cu
index aa272d8..df140ec 100644
--- a/cuda/2d/em.cu
+++ b/cuda/2d/em.cu
@@ -29,10 +29,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/2d/util.h"
 #include "astra/cuda/2d/arith.h"
 
-#ifdef STANDALONE
-#include "testutil.h"
-#endif
-
 #include <cstdio>
 #include <cassert>
 
@@ -168,64 +164,3 @@ float EM::computeDiffNorm()
 
 
 }
-
-#ifdef STANDALONE
-
-using namespace astraCUDA;
-
-int main()
-{
-	float* D_volumeData;
-	float* D_sinoData;
-
-	SDimensions dims;
-	dims.iVolWidth = 1024;
-	dims.iVolHeight = 1024;
-	dims.iProjAngles = 512;
-	dims.iProjDets = 1536;
-	dims.fDetScale = 1.0f;
-	dims.iRaysPerDet = 1;
-	unsigned int volumePitch, sinoPitch;
-
-	allocateVolume(D_volumeData, dims.iVolWidth, dims.iVolHeight, volumePitch);
-	zeroVolume(D_volumeData, volumePitch, dims.iVolWidth, dims.iVolHeight);
-	printf("pitch: %u\n", volumePitch);
-
-	allocateVolume(D_sinoData, dims.iProjDets, dims.iProjAngles, sinoPitch);
-	zeroVolume(D_sinoData, sinoPitch, dims.iProjDets, dims.iProjAngles);
-	printf("pitch: %u\n", sinoPitch);
-	
-	unsigned int y, x;
-	float* sino = loadImage("sino.png", y, x);
-
-	float* img = new float[dims.iVolWidth*dims.iVolHeight];
-
-	copySinogramToDevice(sino, dims.iProjDets, dims.iProjDets, dims.iProjAngles, D_sinoData, sinoPitch);
-
-	float* angle = new float[dims.iProjAngles];
-
-	for (unsigned int i = 0; i < dims.iProjAngles; ++i)
-		angle[i] = i*(M_PI/dims.iProjAngles);
-
-	EM em;
-
-	em.setGeometry(dims, angle);
-	em.init();
-
-	// TODO: Initialize D_volumeData with an unfiltered backprojection
-
-	em.setBuffers(D_volumeData, volumePitch, D_sinoData, sinoPitch);
-
-	em.iterate(25);
-
-
-	delete[] angle;
-
-	copyVolumeFromDevice(img, dims.iVolWidth, dims.iVolWidth, dims.iVolHeight, D_volumeData, volumePitch);
-
-	saveImage("vol.png",dims.iVolHeight,dims.iVolWidth,img);
-
-	return 0;
-}
-
-#endif
diff --git a/cuda/2d/fan_bp.cu b/cuda/2d/fan_bp.cu
index dac3ac2..76d2fb9 100644
--- a/cuda/2d/fan_bp.cu
+++ b/cuda/2d/fan_bp.cu
@@ -28,10 +28,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/2d/util.h"
 #include "astra/cuda/2d/arith.h"
 
-#ifdef STANDALONE
-#include "testutil.h"
-#endif
-
 #include <cstdio>
 #include <cassert>
 #include <iostream>
@@ -50,12 +46,16 @@ const unsigned int g_blockSlices = 16;
 
 const unsigned int g_MaxAngles = 2560;
 
-__constant__ float gC_SrcX[g_MaxAngles];
-__constant__ float gC_SrcY[g_MaxAngles];
-__constant__ float gC_DetSX[g_MaxAngles];
-__constant__ float gC_DetSY[g_MaxAngles];
-__constant__ float gC_DetUX[g_MaxAngles];
-__constant__ float gC_DetUY[g_MaxAngles];
+struct DevFanParams {
+	float fNumC;
+	float fNumX;
+	float fNumY;
+	float fDenC;
+	float fDenX;
+	float fDenY;
+};
+
+__constant__ DevFanParams gC_C[g_MaxAngles];
 
 
 static bool bindProjDataTexture(float* data, unsigned int pitch, unsigned int width, unsigned int height, cudaTextureAddressMode mode = cudaAddressModeBorder)
@@ -74,6 +74,7 @@ static bool bindProjDataTexture(float* data, unsigned int pitch, unsigned int wi
 	return true;
 }
 
+template<bool FBPWEIGHT>
 __global__ void devFanBP(float* D_volData, unsigned int volPitch, unsigned int startAngle, const SDimensions dims, float fOutputScale)
 {
 	const int relX = threadIdx.x;
@@ -96,25 +97,25 @@ __global__ void devFanBP(float* D_volData, unsigned int volPitch, unsigned int s
 	float fVal = 0.0f;
 	float fA = startAngle + 0.5f;
 
-	// TODO: Distance correction?
-
 	for (int angle = startAngle; angle < endAngle; ++angle)
 	{
-		const float fSrcX = gC_SrcX[angle];
-		const float fSrcY = gC_SrcY[angle];
-		const float fDetSX = gC_DetSX[angle];
-		const float fDetSY = gC_DetSY[angle];
-		const float fDetUX = gC_DetUX[angle];
-		const float fDetUY = gC_DetUY[angle];
-
-		const float fXD = fSrcX - fX;
-		const float fYD = fSrcY - fY;
-
-		const float fNum = fDetSY * fXD - fDetSX * fYD + fX*fSrcY - fY*fSrcX;
-		const float fDen = fDetUX * fYD - fDetUY * fXD;
-		
-		const float fT = fNum / fDen;
-		fVal += tex2D(gT_FanProjTexture, fT, fA);
+		const float fNumC = gC_C[angle].fNumC;
+		const float fNumX = gC_C[angle].fNumX;
+		const float fNumY = gC_C[angle].fNumY;
+		const float fDenX = gC_C[angle].fDenX;
+		const float fDenY = gC_C[angle].fDenY;
+
+		const float fNum = fNumC + fNumX * fX + fNumY * fY;
+		const float fDen = (FBPWEIGHT ? 1.0 : gC_C[angle].fDenC) + fDenX * fX + fDenY * fY;
+
+		// Scale factor is the approximate number of rays traversing this pixel,
+		// given by the inverse size of a detector pixel scaled by the magnification
+		// factor of this pixel.
+		// Magnification factor is || u (d-s) || / || u (x-s) ||
+
+		const float fr = __fdividef(1.0f, fDen);
+		const float fT = fNum * fr;
+		fVal += tex2D(gT_FanProjTexture, fT, fA) * (FBPWEIGHT ? fr * fr : fr);
 		fA += 1.0f;
 	}
 
@@ -148,30 +149,27 @@ __global__ void devFanBP_SS(float* D_volData, unsigned int volPitch, unsigned in
 	float fVal = 0.0f;
 	float fA = startAngle + 0.5f;
 
-	// TODO: Distance correction?
-
 	for (int angle = startAngle; angle < endAngle; ++angle)
 	{
-		const float fSrcX = gC_SrcX[angle];
-		const float fSrcY = gC_SrcY[angle];
-		const float fDetSX = gC_DetSX[angle];
-		const float fDetSY = gC_DetSY[angle];
-		const float fDetUX = gC_DetUX[angle];
-		const float fDetUY = gC_DetUY[angle];
+		const float fNumC = gC_C[angle].fNumC;
+		const float fNumX = gC_C[angle].fNumX;
+		const float fNumY = gC_C[angle].fNumY;
+		const float fDenC = gC_C[angle].fDenC;
+		const float fDenX = gC_C[angle].fDenX;
+		const float fDenY = gC_C[angle].fDenY;
 
 		// TODO: Optimize these loops...
 		float fX = fXb;
 		for (int iSubX = 0; iSubX < dims.iRaysPerPixelDim; ++iSubX) {
 			float fY = fYb;
 			for (int iSubY = 0; iSubY < dims.iRaysPerPixelDim; ++iSubY) {
-				const float fXD = fSrcX - fX;
-				const float fYD = fSrcY - fY;
-
-				const float fNum = fDetSY * fXD - fDetSX * fYD + fX*fSrcY - fY*fSrcX;
-				const float fDen = fDetUX * fYD - fDetUY * fXD;
-		
-				const float fT = fNum / fDen;
-				fVal += tex2D(gT_FanProjTexture, fT, fA);
+
+				const float fNum = fNumC + fNumX * fX + fNumY * fY;
+				const float fDen = fDenC + fDenX * fX + fDenY * fY;
+				const float fr = __fdividef(1.0f, fDen);
+
+				const float fT = fNum * fr;
+				fVal += tex2D(gT_FanProjTexture, fT, fA) * fr;
 				fY -= fSubStep;
 			}
 			fX += fSubStep;
@@ -202,77 +200,97 @@ __global__ void devFanBP_SART(float* D_volData, unsigned int volPitch, const SDi
 
 	float* volData = (float*)D_volData;
 
-	// TODO: Distance correction?
-
-	// TODO: Constant memory vs parameters.
-	const float fSrcX = gC_SrcX[0];
-	const float fSrcY = gC_SrcY[0];
-	const float fDetSX = gC_DetSX[0];
-	const float fDetSY = gC_DetSY[0];
-	const float fDetUX = gC_DetUX[0];
-	const float fDetUY = gC_DetUY[0];
+	const float fNumC = gC_C[0].fNumC;
+	const float fNumX = gC_C[0].fNumX;
+	const float fNumY = gC_C[0].fNumY;
+	const float fDenC = gC_C[0].fDenC;
+	const float fDenX = gC_C[0].fDenX;
+	const float fDenY = gC_C[0].fDenY;
 
-	const float fXD = fSrcX - fX;
-	const float fYD = fSrcY - fY;
+	const float fNum = fNumC + fNumX * fX + fNumY * fY;
+	const float fDen = fDenC + fDenX * fX + fDenY * fY;
 
-	const float fNum = fDetSY * fXD - fDetSX * fYD + fX*fSrcY - fY*fSrcX;
-	const float fDen = fDetUX * fYD - fDetUY * fXD;
-		
-	const float fT = fNum / fDen;
+	const float fr = __fdividef(1.0f, fDen);
+	const float fT = fNum * fr;
+	// NB: The scale constant in devBP is cancelled out by the SART weighting
 	const float fVal = tex2D(gT_FanProjTexture, fT, 0.5f);
 
 	volData[Y*volPitch+X] += fVal * fOutputScale;
 }
 
-// Weighted BP for use in fan beam FBP
-// Each pixel/ray is weighted by 1/L^2 where L is the distance to the source.
-__global__ void devFanBP_FBPWeighted(float* D_volData, unsigned int volPitch, unsigned int startAngle, const SDimensions dims, float fOutputScale)
+struct Vec2 {
+	double x;
+	double y;
+	Vec2(double x_, double y_) : x(x_), y(y_) { }
+	Vec2 operator+(const Vec2 &b) const {
+		return Vec2(x + b.x, y + b.y);
+	}
+	Vec2 operator-(const Vec2 &b) const {
+		return Vec2(x - b.x, y - b.y);
+	}
+	Vec2 operator-() const {
+		return Vec2(-x, -y);
+	}
+	double norm() const {
+		return sqrt(x*x + y*y);
+	}
+};
+
+double det2(const Vec2 &a, const Vec2 &b) {
+	return a.x * b.y - a.y * b.x;
+}
+
+
+bool transferConstants(const SFanProjection* angles, unsigned int iProjAngles, bool FBP)
 {
-	const int relX = threadIdx.x;
-	const int relY = threadIdx.y;
+	DevFanParams *p = new DevFanParams[iProjAngles];
 
-	int endAngle = startAngle + g_anglesPerBlock;
-	if (endAngle > dims.iProjAngles)
-		endAngle = dims.iProjAngles;
-	const int X = blockIdx.x * g_blockSlices + relX;
-	const int Y = blockIdx.y * g_blockSliceSize + relY;
+	// We need three values in the kernel:
+	// projected coordinates of pixels on the detector:
+	// || x (s-d) || + ||s d|| / || u (s-x) ||
 
-	if (X >= dims.iVolWidth || Y >= dims.iVolHeight)
-		return;
+	// ray density weighting factor for the adjoint
+	// || u (s-d) || / ( |u| * || u (s-x) || )
 
-	const float fX = ( X - 0.5f*dims.iVolWidth + 0.5f );
-	const float fY = - ( Y - 0.5f*dims.iVolHeight + 0.5f );
+	// fan-beam FBP weighting factor
+	// ( || u s || / || u (s-x) || ) ^ 2
 
-	float* volData = (float*)D_volData;
 
-	float fVal = 0.0f;
-	float fA = startAngle + 0.5f;
 
-	// TODO: Distance correction?
+	for (unsigned int i = 0; i < iProjAngles; ++i) {
+		Vec2 u(angles[i].fDetUX, angles[i].fDetUY);
+		Vec2 s(angles[i].fSrcX, angles[i].fSrcY);
+		Vec2 d(angles[i].fDetSX, angles[i].fDetSY);
 
-	for (int angle = startAngle; angle < endAngle; ++angle)
-	{
-		const float fSrcX = gC_SrcX[angle];
-		const float fSrcY = gC_SrcY[angle];
-		const float fDetSX = gC_DetSX[angle];
-		const float fDetSY = gC_DetSY[angle];
-		const float fDetUX = gC_DetUX[angle];
-		const float fDetUY = gC_DetUY[angle];
-
-		const float fXD = fSrcX - fX;
-		const float fYD = fSrcY - fY;
-
-		const float fNum = fDetSY * fXD - fDetSX * fYD + fX*fSrcY - fY*fSrcX;
-		const float fDen = fDetUX * fYD - fDetUY * fXD;
-
-		const float fWeight = fXD*fXD + fYD*fYD;
-		
-		const float fT = fNum / fDen;
-		fVal += tex2D(gT_FanProjTexture, fT, fA) / fWeight;
-		fA += 1.0f;
+
+
+		double fScale;
+		if (!FBP) {
+			// goal: 1/fDen = || u (s-d) || / ( |u| * || u (s-x) || )
+			// fDen = ( |u| * || u (s-x) || ) / || u (s-d) ||
+			// i.e. scale = |u| /  || u (s-d) ||
+
+			fScale = u.norm() / det2(u, s-d);
+		} else {
+			// goal: 1/fDen = || u s || / || u (s-x) ||
+			// fDen = || u (s-x) || / || u s ||
+			// i.e., scale = 1 / || u s ||
+
+			fScale = 1.0 / det2(u, s);
+		}
+
+		p[i].fNumC = fScale * det2(s,d);
+		p[i].fNumX = fScale * (s-d).y;
+		p[i].fNumY = -fScale * (s-d).x;
+		p[i].fDenC = fScale * det2(u, s); // == 1.0 for FBP
+		p[i].fDenX = fScale * u.y;
+		p[i].fDenY = -fScale * u.x;
 	}
 
-	volData[Y*volPitch+X] += fVal * fOutputScale;
+	// TODO: Check for errors
+	cudaMemcpyToSymbol(gC_C, p, iProjAngles*sizeof(DevFanParams), 0, cudaMemcpyHostToDevice);
+
+	return true;
 }
 
 
@@ -285,21 +303,9 @@ bool FanBP_internal(float* D_volumeData, unsigned int volumePitch,
 
 	bindProjDataTexture(D_projData, projPitch, dims.iProjDets, dims.iProjAngles);
 
-	// transfer angles to constant memory
-	float* tmp = new float[dims.iProjAngles];
-
-#define TRANSFER_TO_CONSTANT(name) do { for (unsigned int i = 0; i < dims.iProjAngles; ++i) tmp[i] = angles[i].f##name ; cudaMemcpyToSymbol(gC_##name, tmp, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); } while (0)
-
-	TRANSFER_TO_CONSTANT(SrcX);
-	TRANSFER_TO_CONSTANT(SrcY);
-	TRANSFER_TO_CONSTANT(DetSX);
-	TRANSFER_TO_CONSTANT(DetSY);
-	TRANSFER_TO_CONSTANT(DetUX);
-	TRANSFER_TO_CONSTANT(DetUY);
-
-#undef TRANSFER_TO_CONSTANT
-
-	delete[] tmp;
+	bool ok = transferConstants(angles, dims.iProjAngles, false);
+	if (!ok)
+		return false;
 
 	dim3 dimBlock(g_blockSlices, g_blockSliceSize);
 	dim3 dimGrid((dims.iVolWidth+g_blockSlices-1)/g_blockSlices,
@@ -312,7 +318,7 @@ bool FanBP_internal(float* D_volumeData, unsigned int volumePitch,
 		if (dims.iRaysPerPixelDim > 1)
 			devFanBP_SS<<<dimGrid, dimBlock, 0, stream>>>(D_volumeData, volumePitch, i, dims, fOutputScale);
 		else
-			devFanBP<<<dimGrid, dimBlock, 0, stream>>>(D_volumeData, volumePitch, i, dims, fOutputScale);
+			devFanBP<false><<<dimGrid, dimBlock, 0, stream>>>(D_volumeData, volumePitch, i, dims, fOutputScale);
 	}
 	cudaThreadSynchronize();
 
@@ -332,21 +338,9 @@ bool FanBP_FBPWeighted_internal(float* D_volumeData, unsigned int volumePitch,
 
 	bindProjDataTexture(D_projData, projPitch, dims.iProjDets, dims.iProjAngles);
 
-	// transfer angles to constant memory
-	float* tmp = new float[dims.iProjAngles];
-
-#define TRANSFER_TO_CONSTANT(name) do { for (unsigned int i = 0; i < dims.iProjAngles; ++i) tmp[i] = angles[i].f##name ; cudaMemcpyToSymbol(gC_##name, tmp, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice); } while (0)
-
-	TRANSFER_TO_CONSTANT(SrcX);
-	TRANSFER_TO_CONSTANT(SrcY);
-	TRANSFER_TO_CONSTANT(DetSX);
-	TRANSFER_TO_CONSTANT(DetSY);
-	TRANSFER_TO_CONSTANT(DetUX);
-	TRANSFER_TO_CONSTANT(DetUY);
-
-#undef TRANSFER_TO_CONSTANT
-
-	delete[] tmp;
+	bool ok = transferConstants(angles, dims.iProjAngles, true);
+	if (!ok)
+		return false;
 
 	dim3 dimBlock(g_blockSlices, g_blockSliceSize);
 	dim3 dimGrid((dims.iVolWidth+g_blockSlices-1)/g_blockSlices,
@@ -356,7 +350,7 @@ bool FanBP_FBPWeighted_internal(float* D_volumeData, unsigned int volumePitch,
 	cudaStreamCreate(&stream);
 
 	for (unsigned int i = 0; i < dims.iProjAngles; i += g_anglesPerBlock) {
-		devFanBP_FBPWeighted<<<dimGrid, dimBlock, 0, stream>>>(D_volumeData, volumePitch, i, dims, fOutputScale);
+		devFanBP<true><<<dimGrid, dimBlock, 0, stream>>>(D_volumeData, volumePitch, i, dims, fOutputScale);
 	}
 	cudaThreadSynchronize();
 
@@ -377,17 +371,9 @@ bool FanBP_SART(float* D_volumeData, unsigned int volumePitch,
 	// only one angle
 	bindProjDataTexture(D_projData, projPitch, dims.iProjDets, 1, cudaAddressModeClamp);
 
-	// transfer angle to constant memory
-#define TRANSFER_TO_CONSTANT(name) do { cudaMemcpyToSymbol(gC_##name, &(angles[angle].f##name), sizeof(float), 0, cudaMemcpyHostToDevice); } while (0)
-
-	TRANSFER_TO_CONSTANT(SrcX);
-	TRANSFER_TO_CONSTANT(SrcY);
-	TRANSFER_TO_CONSTANT(DetSX);
-	TRANSFER_TO_CONSTANT(DetSY);
-	TRANSFER_TO_CONSTANT(DetUX);
-	TRANSFER_TO_CONSTANT(DetUY);
-
-#undef TRANSFER_TO_CONSTANT
+	bool ok = transferConstants(angles + angle, 1, false);
+	if (!ok)
+		return false;
 
 	dim3 dimBlock(g_blockSlices, g_blockSliceSize);
 	dim3 dimGrid((dims.iVolWidth+g_blockSlices-1)/g_blockSlices,
@@ -448,66 +434,3 @@ bool FanBP_FBPWeighted(float* D_volumeData, unsigned int volumePitch,
 
 
 }
-
-#ifdef STANDALONE
-
-using namespace astraCUDA;
-
-int main()
-{
-	float* D_volumeData;
-	float* D_projData;
-
-	SDimensions dims;
-	dims.iVolWidth = 128;
-	dims.iVolHeight = 128;
-	dims.iProjAngles = 180;
-	dims.iProjDets = 256;
-	dims.fDetScale = 1.0f;
-	dims.iRaysPerDet = 1;
-	unsigned int volumePitch, projPitch;
-
-	SFanProjection projs[180];
-
-	projs[0].fSrcX = 0.0f;
-	projs[0].fSrcY = 1536.0f;
-	projs[0].fDetSX = 128.0f;
-	projs[0].fDetSY = -512.0f;
-	projs[0].fDetUX = -1.0f;
-	projs[0].fDetUY = 0.0f;
-
-#define ROTATE0(name,i,alpha) do { projs[i].f##name##X = projs[0].f##name##X * cos(alpha) - projs[0].f##name##Y * sin(alpha); projs[i].f##name##Y = projs[0].f##name##X * sin(alpha) + projs[0].f##name##Y * cos(alpha); } while(0)
-
-	for (int i = 1; i < 180; ++i) {
-		ROTATE0(Src, i, i*2*M_PI/180);
-		ROTATE0(DetS, i, i*2*M_PI/180);
-		ROTATE0(DetU, i, i*2*M_PI/180);
-	}
-
-#undef ROTATE0
-
-	allocateVolume(D_volumeData, dims.iVolWidth, dims.iVolHeight, volumePitch);
-	printf("pitch: %u\n", volumePitch);
-
-	allocateVolume(D_projData, dims.iProjDets, dims.iProjAngles, projPitch);
-	printf("pitch: %u\n", projPitch);
-
-	unsigned int y, x;
-	float* sino = loadImage("sino.png", y, x);
-
-	float* img = new float[dims.iVolWidth*dims.iVolHeight];
-
-	memset(img, 0, dims.iVolWidth*dims.iVolHeight*sizeof(float));
-
-	copyVolumeToDevice(img, dims.iVolWidth, dims.iVolWidth, dims.iVolHeight, D_volumeData, volumePitch);
-	copySinogramToDevice(sino, dims.iProjDets, dims.iProjDets, dims.iProjAngles, D_projData, projPitch);
-
-	FanBP(D_volumeData, volumePitch, D_projData, projPitch, dims, projs, 1.0f);
-
-	copyVolumeFromDevice(img, dims.iVolWidth, dims.iVolWidth, dims.iVolHeight, D_volumeData, volumePitch);
-
-	saveImage("vol.png",dims.iVolHeight,dims.iVolWidth,img);
-
-	return 0;
-}
-#endif
diff --git a/cuda/2d/fan_fp.cu b/cuda/2d/fan_fp.cu
index 3479650..60c02f8 100644
--- a/cuda/2d/fan_fp.cu
+++ b/cuda/2d/fan_fp.cu
@@ -28,10 +28,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/2d/util.h"
 #include "astra/cuda/2d/arith.h"
 
-#ifdef STANDALONE
-#include "testutil.h"
-#endif
-
 #include <cstdio>
 #include <cassert>
 #include <iostream>
@@ -308,84 +304,3 @@ bool FanFP(float* D_volumeData, unsigned int volumePitch,
 }
 
 }
-
-#ifdef STANDALONE
-
-using namespace astraCUDA;
-
-int main()
-{
-	float* D_volumeData;
-	float* D_projData;
-
-	SDimensions dims;
-	dims.iVolWidth = 128;
-	dims.iVolHeight = 128;
-	dims.iProjAngles = 180;
-	dims.iProjDets = 256;
-	dims.fDetScale = 1.0f;
-	dims.iRaysPerDet = 1;
-	unsigned int volumePitch, projPitch;
-
-	SFanProjection projs[180];
-
-	projs[0].fSrcX = 0.0f;
-	projs[0].fSrcY = 1536.0f;
-	projs[0].fDetSX = 128.0f;
-	projs[0].fDetSY = -512.0f;
-	projs[0].fDetUX = -1.0f;
-	projs[0].fDetUY = 0.0f;
-
-#define ROTATE0(name,i,alpha) do { projs[i].f##name##X = projs[0].f##name##X * cos(alpha) - projs[0].f##name##Y * sin(alpha); projs[i].f##name##Y = projs[0].f##name##X * sin(alpha) + projs[0].f##name##Y * cos(alpha); } while(0)
-
-	for (int i = 1; i < 180; ++i) {
-		ROTATE0(Src, i, i*2*M_PI/180);
-		ROTATE0(DetS, i, i*2*M_PI/180);
-		ROTATE0(DetU, i, i*2*M_PI/180);
-	}
-
-#undef ROTATE0
-
-	allocateVolume(D_volumeData, dims.iVolWidth, dims.iVolHeight, volumePitch);
-	printf("pitch: %u\n", volumePitch);
-
-	allocateVolume(D_projData, dims.iProjDets, dims.iProjAngles, projPitch);
-	printf("pitch: %u\n", projPitch);
-
-	unsigned int y, x;
-	float* img = loadImage("phantom128.png", y, x);
-
-	float* sino = new float[dims.iProjAngles * dims.iProjDets];
-
-	memset(sino, 0, dims.iProjAngles * dims.iProjDets * sizeof(float));
-
-	copyVolumeToDevice(img, dims.iVolWidth, dims.iVolWidth, dims.iVolHeight, D_volumeData, volumePitch);
-	copySinogramToDevice(sino, dims.iProjDets, dims.iProjDets, dims.iProjAngles, D_projData, projPitch);
-
-	float* angle = new float[dims.iProjAngles];
-
-	for (unsigned int i = 0; i < dims.iProjAngles; ++i)
-		angle[i] = i*(M_PI/dims.iProjAngles);
-
-	FanFP(D_volumeData, volumePitch, D_projData, projPitch, dims, projs, 1.0f);
-
-	delete[] angle;
-
-	copySinogramFromDevice(sino, dims.iProjDets, dims.iProjDets, dims.iProjAngles, D_projData, projPitch);
-
-	float s = 0.0f;
-	for (unsigned int y = 0; y < dims.iProjAngles; ++y)
-		for (unsigned int x = 0; x < dims.iProjDets; ++x)
-			s += sino[y*dims.iProjDets+x] * sino[y*dims.iProjDets+x];
-	printf("cpu norm: %f\n", s);
-
-	//zeroVolume(D_projData, projPitch, dims.iProjDets, dims.iProjAngles);
-	s = dotProduct2D(D_projData, projPitch, dims.iProjDets, dims.iProjAngles);
-	printf("gpu norm: %f\n", s);
-
-	saveImage("sino.png",dims.iProjAngles,dims.iProjDets,sino);
-
-
-	return 0;
-}
-#endif
diff --git a/cuda/2d/fbp.cu b/cuda/2d/fbp.cu
index a5b8a7a..4fc3983 100644
--- a/cuda/2d/fbp.cu
+++ b/cuda/2d/fbp.cu
@@ -58,7 +58,8 @@ int FBP::calcFourierFilterSize(int _iDetectorCount)
 FBP::FBP() : ReconAlgo()
 {
 	D_filter = 0;
-
+	m_bShortScan = false;
+	fReconstructionScale = 1.0f;
 }
 
 FBP::~FBP()
@@ -72,6 +73,8 @@ void FBP::reset()
 		freeComplexOnDevice((cufftComplex *)D_filter);
 		D_filter = 0;
 	}
+	m_bShortScan = false;
+	fReconstructionScale = 1.0f;
 }
 
 bool FBP::init()
@@ -79,6 +82,12 @@ bool FBP::init()
 	return true;
 }
 
+bool FBP::setReconstructionScale(float fScale)
+{
+	fReconstructionScale = fScale;
+	return true;
+}
+
 bool FBP::setFilter(const astra::SFilterConfig &_cfg)
 {
 	if (D_filter)
@@ -292,7 +301,7 @@ bool FBP::iterate(unsigned int iterations)
 
 		astraCUDA3d::FDK_PreWeight(tmp, fOriginSource,
 		              fOriginDetector, 0.0f,
-		              fFanDetSize, 1.0f, /* fPixelSize */ 1.0f,
+		              fFanDetSize, 1.0f,
 		              m_bShortScan, dims3d, pfAngles);
 	} else {
 		// TODO: How should different detector pixel size in different
@@ -319,17 +328,14 @@ bool FBP::iterate(unsigned int iterations)
 	}
 
 	if (fanProjs) {
-		float fOutputScale = 1.0 / (/*fPixelSize * fPixelSize * fPixelSize * */ fFanDetSize * fFanDetSize);
-
-		ok = FanBP_FBPWeighted(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, fanProjs, fOutputScale);
+		ok = FanBP_FBPWeighted(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, fanProjs, fProjectorScale * fReconstructionScale);
 
 	} else {
 		// scale by number of angles. For the fan-beam case, this is already
 		// handled by FDK_PreWeight
 		float fOutputScale = (M_PI / 2.0f) / (float)dims.iProjAngles;
-		//fOutputScale /= fDetSize * fDetSize;
 
-		ok = BP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, parProjs, fOutputScale);
+		ok = BP(D_volumeData, volumePitch, D_sinoData, sinoPitch, dims, parProjs, fOutputScale * fProjectorScale * fReconstructionScale);
 	}
 	if(!ok)
 	{
diff --git a/cuda/2d/fft.cu b/cuda/2d/fft.cu
index 2e94b79..8361ad2 100644
--- a/cuda/2d/fft.cu
+++ b/cuda/2d/fft.cu
@@ -314,210 +314,3 @@ void genCuFFTFilter(const SFilterConfig &_cfg, int _iProjectionCount,
 
 
 }
-
-
-#ifdef STANDALONE
-
-__global__ static void doubleFourierOutput_kernel(int _iProjectionCount,
-                                                  int _iDetectorCount,
-                                                  cufftComplex* _pFourierOutput)
-{
-	int iIndex = threadIdx.x + blockIdx.x * blockDim.x;
-	int iProjectionIndex = iIndex / _iDetectorCount;
-	int iDetectorIndex = iIndex % _iDetectorCount;
-
-	if(iProjectionIndex >= _iProjectionCount)
-	{
-		return;
-	}
-
-	if(iDetectorIndex <= (_iDetectorCount / 2))
-	{
-		return;
-	}
-
-	int iOtherDetectorIndex = _iDetectorCount - iDetectorIndex;
-
-	_pFourierOutput[iProjectionIndex * _iDetectorCount + iDetectorIndex].x = _pFourierOutput[iProjectionIndex * _iDetectorCount + iOtherDetectorIndex].x;
-	_pFourierOutput[iProjectionIndex * _iDetectorCount + iDetectorIndex].y = -_pFourierOutput[iProjectionIndex * _iDetectorCount + iOtherDetectorIndex].y;
-}
-
-static void doubleFourierOutput(int _iProjectionCount, int _iDetectorCount,
-                                cufftComplex * _pFourierOutput)
-{
-	const int iBlockSize = 256;
-	int iElementCount = _iProjectionCount * _iDetectorCount;
-	int iBlockCount = (iElementCount + iBlockSize - 1) / iBlockSize;
-
-	doubleFourierOutput_kernel<<< iBlockCount, iBlockSize >>>(_iProjectionCount,
-	                                                          _iDetectorCount,
-	                                                          _pFourierOutput);
-	CHECK_ERROR("doubleFourierOutput_kernel failed");
-}
-
-
-
-static void writeToMatlabFile(const char * _fileName, const float * _pfData,
-                              int _iRowCount, int _iColumnCount)
-{
-	std::ofstream out(_fileName);
-
-	for(int iRowIndex = 0; iRowIndex < _iRowCount; iRowIndex++)
-	{
-		for(int iColumnIndex = 0; iColumnIndex < _iColumnCount; iColumnIndex++)
-		{
-			out << _pfData[iColumnIndex + iRowIndex * _iColumnCount] << " ";
-		}
-
-		out << std::endl;
-	}
-}
-
-static void convertComplexToRealImg(const cufftComplex * _pComplex,
-                                    int _iElementCount,
-                                    float * _pfReal, float * _pfImaginary)
-{
-	for(int iIndex = 0; iIndex < _iElementCount; iIndex++)
-	{
-		_pfReal[iIndex] = _pComplex[iIndex].x;
-		_pfImaginary[iIndex] = _pComplex[iIndex].y;
-	}
-}
-
-void testCudaFFT()
-{
-	const int iProjectionCount = 2;
-	const int iDetectorCount = 1024;
-	const int iTotalElementCount = iProjectionCount * iDetectorCount;
-
-	float * pfHostProj = new float[iTotalElementCount];
-	memset(pfHostProj, 0, sizeof(float) * iTotalElementCount);
-
-	for(int iProjectionIndex = 0; iProjectionIndex < iProjectionCount; iProjectionIndex++)
-	{
-		for(int iDetectorIndex = 0; iDetectorIndex < iDetectorCount; iDetectorIndex++)
-		{
-//			int
-
-//			pfHostProj[iIndex] = (float)rand() / (float)RAND_MAX;
-		}
-	}
-
-	writeToMatlabFile("proj.mat", pfHostProj, iProjectionCount, iDetectorCount);
-
-	float * pfDevProj = NULL;
-	SAFE_CALL(cudaMalloc((void **)&pfDevProj, sizeof(float) * iTotalElementCount));
-	SAFE_CALL(cudaMemcpy(pfDevProj, pfHostProj, sizeof(float) * iTotalElementCount, cudaMemcpyHostToDevice));
-
-	cufftComplex * pDevFourProj = NULL;
-	SAFE_CALL(cudaMalloc((void **)&pDevFourProj, sizeof(cufftComplex) * iTotalElementCount));
-
-	cufftHandle plan;
-	cufftResult result;
-
-	result = cufftPlan1d(&plan, iDetectorCount, CUFFT_R2C, iProjectionCount);
-	if(result != CUFFT_SUCCESS)
-	{
-		ASTRA_ERROR("Failed to plan 1d r2c fft");
-	}
-
-	result = cufftExecR2C(plan, pfDevProj, pDevFourProj);
-	if(result != CUFFT_SUCCESS)
-	{
-		ASTRA_ERROR("Failed to exec 1d r2c fft");
-	}
-
-	cufftDestroy(plan);
-
-	doubleFourierOutput(iProjectionCount, iDetectorCount, pDevFourProj);
-
-	cufftComplex * pHostFourProj = new cufftComplex[iTotalElementCount];
-	SAFE_CALL(cudaMemcpy(pHostFourProj, pDevFourProj, sizeof(cufftComplex) * iTotalElementCount, cudaMemcpyDeviceToHost));
-
-	float * pfHostFourProjReal = new float[iTotalElementCount];
-	float * pfHostFourProjImaginary = new float[iTotalElementCount];
-
-	convertComplexToRealImg(pHostFourProj, iTotalElementCount, pfHostFourProjReal, pfHostFourProjImaginary);
-
-	writeToMatlabFile("proj_four_real.mat", pfHostFourProjReal, iProjectionCount, iDetectorCount);
-	writeToMatlabFile("proj_four_imaginary.mat", pfHostFourProjImaginary, iProjectionCount, iDetectorCount);
-
-	float * pfDevInFourProj = NULL;
-	SAFE_CALL(cudaMalloc((void **)&pfDevInFourProj, sizeof(float) * iTotalElementCount));
-
-	result = cufftPlan1d(&plan, iDetectorCount, CUFFT_C2R, iProjectionCount);
-	if(result != CUFFT_SUCCESS)
-	{
-		ASTRA_ERROR("Failed to plan 1d c2r fft");
-	}
-
-	result = cufftExecC2R(plan, pDevFourProj, pfDevInFourProj);
-	if(result != CUFFT_SUCCESS)
-	{
-		ASTRA_ERROR("Failed to exec 1d c2r fft");
-	}
-
-	cufftDestroy(plan);
-
-	rescaleInverseFourier(iProjectionCount, iDetectorCount, pfDevInFourProj);
-
-	float * pfHostInFourProj = new float[iTotalElementCount];
-	SAFE_CALL(cudaMemcpy(pfHostInFourProj, pfDevInFourProj, sizeof(float) * iTotalElementCount, cudaMemcpyDeviceToHost));
-
-	writeToMatlabFile("in_four.mat", pfHostInFourProj, iProjectionCount, iDetectorCount);
-
-	SAFE_CALL(cudaFree(pDevFourProj));
-	SAFE_CALL(cudaFree(pfDevProj));
-
-	delete [] pfHostInFourProj;
-	delete [] pfHostFourProjReal;
-	delete [] pfHostFourProjImaginary;
-	delete [] pfHostProj;
-	delete [] pHostFourProj;
-}
-
-void downloadDebugFilterComplex(float * _pfHostSinogram, int _iProjectionCount,
-                                int _iDetectorCount,
-                                cufftComplex * _pDevFilter,
-                                int _iFilterDetCount)
-{
-	cufftComplex * pHostFilter = NULL;
-	size_t complMemSize = sizeof(cufftComplex) * _iFilterDetCount * _iProjectionCount;
-	pHostFilter = (cufftComplex *)malloc(complMemSize);
-	SAFE_CALL(cudaMemcpy(pHostFilter, _pDevFilter, complMemSize, cudaMemcpyDeviceToHost));
-
-	for(int iTargetProjIndex = 0; iTargetProjIndex < _iProjectionCount; iTargetProjIndex++)
-	{
-		for(int iTargetDetIndex = 0; iTargetDetIndex < min(_iDetectorCount, _iFilterDetCount); iTargetDetIndex++)
-		{
-			cufftComplex source = pHostFilter[iTargetDetIndex + iTargetProjIndex * _iFilterDetCount];
-			float fReadValue = sqrtf(source.x * source.x + source.y * source.y);
-			_pfHostSinogram[iTargetDetIndex + iTargetProjIndex * _iDetectorCount] = fReadValue;
-		}
-	}
-
-	free(pHostFilter);
-}
-
-void downloadDebugFilterReal(float * _pfHostSinogram, int _iProjectionCount,
-                             int _iDetectorCount, float * _pfDevFilter,
-                             int _iFilterDetCount)
-{
-	float * pfHostFilter = NULL;
-	size_t memSize = sizeof(float) * _iFilterDetCount * _iProjectionCount;
-	pfHostFilter = (float *)malloc(memSize);
-	SAFE_CALL(cudaMemcpy(pfHostFilter, _pfDevFilter, memSize, cudaMemcpyDeviceToHost));
-
-	for(int iTargetProjIndex = 0; iTargetProjIndex < _iProjectionCount; iTargetProjIndex++)
-	{
-		for(int iTargetDetIndex = 0; iTargetDetIndex < min(_iDetectorCount, _iFilterDetCount); iTargetDetIndex++)
-		{
-			float fSource = pfHostFilter[iTargetDetIndex + iTargetProjIndex * _iFilterDetCount];
-			_pfHostSinogram[iTargetDetIndex + iTargetProjIndex * _iDetectorCount] = fSource;
-		}
-	}
-
-	free(pfHostFilter);
-}
-
-#endif
diff --git a/cuda/2d/par_bp.cu b/cuda/2d/par_bp.cu
index 09a6554..f080abb 100644
--- a/cuda/2d/par_bp.cu
+++ b/cuda/2d/par_bp.cu
@@ -28,10 +28,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/2d/util.h"
 #include "astra/cuda/2d/arith.h"
 
-#ifdef STANDALONE
-#include "testutil.h"
-#endif
-
 #include <cstdio>
 #include <cassert>
 #include <iostream>
@@ -53,6 +49,7 @@ const unsigned int g_MaxAngles = 2560;
 __constant__ float gC_angle_scaled_sin[g_MaxAngles];
 __constant__ float gC_angle_scaled_cos[g_MaxAngles];
 __constant__ float gC_angle_offset[g_MaxAngles];
+__constant__ float gC_angle_scale[g_MaxAngles];
 
 static bool bindProjDataTexture(float* data, unsigned int pitch, unsigned int width, unsigned int height, cudaTextureAddressMode mode = cudaAddressModeBorder)
 {
@@ -70,6 +67,7 @@ static bool bindProjDataTexture(float* data, unsigned int pitch, unsigned int wi
 	return true;
 }
 
+// TODO: Templated version with/without scale? (Or only the global outputscale)
 __global__ void devBP(float* D_volData, unsigned int volPitch, unsigned int startAngle, const SDimensions dims, float fOutputScale)
 {
 	const int relX = threadIdx.x;
@@ -97,9 +95,10 @@ __global__ void devBP(float* D_volData, unsigned int volPitch, unsigned int star
 		const float scaled_cos_theta = gC_angle_scaled_cos[angle];
 		const float scaled_sin_theta = gC_angle_scaled_sin[angle];
 		const float TOffset = gC_angle_offset[angle];
+		const float scale = gC_angle_scale[angle];
 
 		const float fT = fX * scaled_cos_theta - fY * scaled_sin_theta + TOffset;
-		fVal += tex2D(gT_projTexture, fT, fA);
+		fVal += tex2D(gT_projTexture, fT, fA) * scale;
 		fA += 1.0f;
 	}
 
@@ -138,6 +137,7 @@ __global__ void devBP_SS(float* D_volData, unsigned int volPitch, unsigned int s
 		const float cos_theta = gC_angle_scaled_cos[angle];
 		const float sin_theta = gC_angle_scaled_sin[angle];
 		const float TOffset = gC_angle_offset[angle];
+		const float scale = gC_angle_scale[angle];
 
 		float fT = fX * cos_theta - fY * sin_theta + TOffset;
 
@@ -145,7 +145,7 @@ __global__ void devBP_SS(float* D_volData, unsigned int volPitch, unsigned int s
 			float fTy = fT;
 			fT += fSubStep * cos_theta;
 			for (int iSubY = 0; iSubY < dims.iRaysPerPixelDim; ++iSubY) {
-				fVal += tex2D(gT_projTexture, fTy, fA);
+				fVal += tex2D(gT_projTexture, fTy, fA) * scale;
 				fTy -= fSubStep * sin_theta;
 			}
 		}
@@ -172,6 +172,8 @@ __global__ void devBP_SART(float* D_volData, unsigned int volPitch, float offset
 	const float fT = fX * angle_cos - fY * angle_sin + offset;
 	const float fVal = tex2D(gT_projTexture, fT, 0.5f);
 
+	// NB: The 'scale' constant in devBP is cancelled out by the SART weighting
+
 	D_volData[Y*volPitch+X] += fVal * fOutputScale;
 }
 
@@ -186,27 +188,34 @@ bool BP_internal(float* D_volumeData, unsigned int volumePitch,
 	float* angle_scaled_sin = new float[dims.iProjAngles];
 	float* angle_scaled_cos = new float[dims.iProjAngles];
 	float* angle_offset = new float[dims.iProjAngles];
+	float* angle_scale = new float[dims.iProjAngles];
 
 	bindProjDataTexture(D_projData, projPitch, dims.iProjDets, dims.iProjAngles);
 
 	for (unsigned int i = 0; i < dims.iProjAngles; ++i) {
 		double d = angles[i].fDetUX * angles[i].fRayY - angles[i].fDetUY * angles[i].fRayX;
 		angle_scaled_cos[i] = angles[i].fRayY / d;
-		angle_scaled_sin[i] = -angles[i].fRayX / d; // TODO: Check signs
+		angle_scaled_sin[i] = -angles[i].fRayX / d;
 		angle_offset[i] = (angles[i].fDetSY * angles[i].fRayX - angles[i].fDetSX * angles[i].fRayY) / d;
+		angle_scale[i] = sqrt(angles[i].fRayX * angles[i].fRayX + angles[i].fRayY * angles[i].fRayY) / abs(d);
 	}
+	//fprintf(stderr, "outputscale in BP_internal: %f, %f\n", fOutputScale, angle_scale[0]);
+	//fprintf(stderr, "ray in BP_internal: %f,%f (length %f)\n", angles[0].fRayX, angles[0].fRayY, sqrt(angles[0].fRayX * angles[0].fRayX + angles[0].fRayY * angles[0].fRayY));
 
 	cudaError_t e1 = cudaMemcpyToSymbol(gC_angle_scaled_sin, angle_scaled_sin, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice);
 	cudaError_t e2 = cudaMemcpyToSymbol(gC_angle_scaled_cos, angle_scaled_cos, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice);
 	cudaError_t e3 = cudaMemcpyToSymbol(gC_angle_offset, angle_offset, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice);
+	cudaError_t e4 = cudaMemcpyToSymbol(gC_angle_scale, angle_scale, dims.iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice);
 	assert(e1 == cudaSuccess);
 	assert(e2 == cudaSuccess);
 	assert(e3 == cudaSuccess);
+	assert(e4 == cudaSuccess);
 
 
 	delete[] angle_scaled_sin;
 	delete[] angle_scaled_cos;
 	delete[] angle_offset;
+	delete[] angle_scale;
 
 	dim3 dimBlock(g_blockSlices, g_blockSliceSize);
 	dim3 dimGrid((dims.iVolWidth+g_blockSlices-1)/g_blockSlices,
@@ -267,6 +276,8 @@ bool BP_SART(float* D_volumeData, unsigned int volumePitch,
 	float angle_scaled_cos = angles[angle].fRayY / d;
 	float angle_scaled_sin = -angles[angle].fRayX / d; // TODO: Check signs
 	float angle_offset = (angles[angle].fDetSY * angles[angle].fRayX - angles[angle].fDetSX * angles[angle].fRayY) / d;
+	// NB: The adjoint scaling factor from regular BP is cancelled out by the SART weighting
+	//fOutputScale *= sqrt(angles[angle].fRayX * angles[angle].fRayX + angles[angle].fRayY * angles[angle].fRayY) / abs(d);
 
 	dim3 dimBlock(g_blockSlices, g_blockSliceSize);
 	dim3 dimGrid((dims.iVolWidth+g_blockSlices-1)/g_blockSlices,
@@ -282,55 +293,3 @@ bool BP_SART(float* D_volumeData, unsigned int volumePitch,
 
 
 }
-
-#ifdef STANDALONE
-
-using namespace astraCUDA;
-
-int main()
-{
-	float* D_volumeData;
-	float* D_projData;
-
-	SDimensions dims;
-	dims.iVolWidth = 1024;
-	dims.iVolHeight = 1024;
-	dims.iProjAngles = 512;
-	dims.iProjDets = 1536;
-	dims.fDetScale = 1.0f;
-	dims.iRaysPerDet = 1;
-
-	unsigned int volumePitch, projPitch;
-
-	allocateVolume(D_volumeData, dims.iVolWidth, dims.iVolHeight, volumePitch);
-	printf("pitch: %u\n", volumePitch);
-
-	allocateVolume(D_projData, dims.iProjDets, dims.iProjAngles, projPitch);
-	printf("pitch: %u\n", projPitch);
-
-	unsigned int y, x;
-	float* sino = loadImage("sino.png", y, x);
-
-	float* img = new float[dims.iVolWidth*dims.iVolHeight];
-
-	memset(img, 0, dims.iVolWidth*dims.iVolHeight*sizeof(float));
-
-	copyVolumeToDevice(img, dims.iVolWidth, dims.iVolWidth, dims.iVolHeight, D_volumeData, volumePitch);
-	copySinogramToDevice(sino, dims.iProjDets, dims.iProjDets, dims.iProjAngles, D_projData, projPitch);
-
-	float* angle = new float[dims.iProjAngles];
-
-	for (unsigned int i = 0; i < dims.iProjAngles; ++i)
-		angle[i] = i*(M_PI/dims.iProjAngles);
-
-	BP(D_volumeData, volumePitch, D_projData, projPitch, dims, angle, 0, 1.0f);
-
-	delete[] angle;
-
-	copyVolumeFromDevice(img, dims.iVolWidth, dims.iVolWidth, dims.iVolHeight, D_volumeData, volumePitch);
-
-	saveImage("vol.png",dims.iVolHeight,dims.iVolWidth,img);
-
-	return 0;
-}
-#endif
diff --git a/cuda/2d/par_fp.cu b/cuda/2d/par_fp.cu
index 0835301..ea436c3 100644
--- a/cuda/2d/par_fp.cu
+++ b/cuda/2d/par_fp.cu
@@ -28,10 +28,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/2d/util.h"
 #include "astra/cuda/2d/arith.h"
 
-#ifdef STANDALONE
-#include "testutil.h"
-#endif
-
 #include <cstdio>
 #include <cassert>
 #include <iostream>
@@ -115,10 +111,9 @@ __global__ void FPhorizontal_simple(float* D_projData, unsigned int projPitch, u
 	float fSliceStep = cos_theta / sin_theta;
 	float fDistCorr;
 	if (sin_theta > 0.0f)
-		fDistCorr = -fDetStep;
+		fDistCorr = outputScale / sin_theta;
 	else
-		fDistCorr = fDetStep;
-	fDistCorr *= outputScale;
+		fDistCorr = -outputScale / sin_theta;
 
 	float fVal = 0.0f;
 	// project detector on slice
@@ -193,10 +188,9 @@ __global__ void FPvertical_simple(float* D_projData, unsigned int projPitch, uns
 	float fSliceStep = sin_theta / cos_theta;
 	float fDistCorr;
 	if (cos_theta < 0.0f)
-		fDistCorr = -fDetStep;
+		fDistCorr = -outputScale / cos_theta; 
 	else
-		fDistCorr = fDetStep;
-	fDistCorr *= outputScale;
+		fDistCorr = outputScale / cos_theta;
 
 	float fVal = 0.0f;
 	float fP = (detector - 0.5f*dims.iProjDets + 0.5f - gC_angle_offset[angle]) * fDetStep + (startSlice - 0.5f*dims.iVolHeight + 0.5f) * fSliceStep + 0.5f*dims.iVolWidth - 0.5f + 0.5f;
@@ -375,65 +369,3 @@ bool FP(float* D_volumeData, unsigned int volumePitch,
 
 
 }
-
-#ifdef STANDALONE
-
-using namespace astraCUDA;
-
-int main()
-{
-	float* D_volumeData;
-	float* D_projData;
-
-	SDimensions dims;
-	dims.iVolWidth = 1024;
-	dims.iVolHeight = 1024;
-	dims.iProjAngles = 512;
-	dims.iProjDets = 1536;
-	dims.fDetScale = 1.0f;
-	dims.iRaysPerDet = 1;
-	unsigned int volumePitch, projPitch;
-
-	allocateVolume(D_volumeData, dims.iVolWidth, dims.iVolHeight, volumePitch);
-	printf("pitch: %u\n", volumePitch);
-
-	allocateVolume(D_projData, dims.iProjDets, dims.iProjAngles, projPitch);
-	printf("pitch: %u\n", projPitch);
-
-	unsigned int y, x;
-	float* img = loadImage("phantom.png", y, x);
-
-	float* sino = new float[dims.iProjAngles * dims.iProjDets];
-
-	memset(sino, 0, dims.iProjAngles * dims.iProjDets * sizeof(float));
-
-	copyVolumeToDevice(img, dims.iVolWidth, dims.iVolWidth, dims.iVolHeight, D_volumeData, volumePitch);
-	copySinogramToDevice(sino, dims.iProjDets, dims.iProjDets, dims.iProjAngles, D_projData, projPitch);
-
-	float* angle = new float[dims.iProjAngles];
-
-	for (unsigned int i = 0; i < dims.iProjAngles; ++i)
-		angle[i] = i*(M_PI/dims.iProjAngles);
-
-	FP(D_volumeData, volumePitch, D_projData, projPitch, dims, angle, 0, 1.0f);
-
-	delete[] angle;
-
-	copySinogramFromDevice(sino, dims.iProjDets, dims.iProjDets, dims.iProjAngles, D_projData, projPitch);
-
-	float s = 0.0f;
-	for (unsigned int y = 0; y < dims.iProjAngles; ++y)
-		for (unsigned int x = 0; x < dims.iProjDets; ++x)
-			s += sino[y*dims.iProjDets+x] * sino[y*dims.iProjDets+x];
-	printf("cpu norm: %f\n", s);
-
-	//zeroVolume(D_projData, projPitch, dims.iProjDets, dims.iProjAngles);
-	s = dotProduct2D(D_projData, projPitch, dims.iProjDets, dims.iProjAngles, 1, 0);
-	printf("gpu norm: %f\n", s);
-
-	saveImage("sino.png",dims.iProjAngles,dims.iProjDets,sino);
-
-
-	return 0;
-}
-#endif
diff --git a/cuda/2d/sart.cu b/cuda/2d/sart.cu
index 64973ba..12ad6df 100644
--- a/cuda/2d/sart.cu
+++ b/cuda/2d/sart.cu
@@ -254,11 +254,11 @@ bool SART::callFP_SART(float* D_volumeData, unsigned int volumePitch,
 	if (parProjs) {
 		assert(!fanProjs);
 		return FP(D_volumeData, volumePitch, D_projData, projPitch,
-		          d, &parProjs[angle], outputScale);
+		          d, &parProjs[angle], outputScale * fProjectorScale);
 	} else {
 		assert(fanProjs);
 		return FanFP(D_volumeData, volumePitch, D_projData, projPitch,
-		             d, &fanProjs[angle], outputScale);
+		             d, &fanProjs[angle], outputScale * fProjectorScale);
 	}
 }
 
@@ -266,6 +266,7 @@ bool SART::callBP_SART(float* D_volumeData, unsigned int volumePitch,
                        float* D_projData, unsigned int projPitch,
                        unsigned int angle, float outputScale)
 {
+	// NB: No fProjectorScale here, as that it is cancelled out in the SART weighting
 	if (parProjs) {
 		assert(!fanProjs);
 		return BP_SART(D_volumeData, volumePitch, D_projData, projPitch,
diff --git a/cuda/2d/sirt.cu b/cuda/2d/sirt.cu
index 2621490..2c5fdc9 100644
--- a/cuda/2d/sirt.cu
+++ b/cuda/2d/sirt.cu
@@ -29,10 +29,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/2d/util.h"
 #include "astra/cuda/2d/arith.h"
 
-#ifdef STANDALONE
-#include "testutil.h"
-#endif
-
 #include <cstdio>
 #include <cassert>
 
@@ -302,62 +298,3 @@ float SIRT::computeDiffNorm()
 
 
 }
-
-#ifdef STANDALONE
-
-using namespace astraCUDA;
-
-int main()
-{
-	float* D_volumeData;
-	float* D_sinoData;
-
-	SDimensions dims;
-	dims.iVolWidth = 1024;
-	dims.iVolHeight = 1024;
-	dims.iProjAngles = 512;
-	dims.iProjDets = 1536;
-	dims.fDetScale = 1.0f;
-	dims.iRaysPerDet = 1;
-	unsigned int volumePitch, sinoPitch;
-
-	allocateVolume(D_volumeData, dims.iVolWidth, dims.iVolHeight, volumePitch);
-	zeroVolume(D_volumeData, volumePitch, dims.iVolWidth, dims.iVolHeight);
-	printf("pitch: %u\n", volumePitch);
-
-	allocateVolume(D_sinoData, dims.iProjDets, dims.iProjAngles, sinoPitch);
-	zeroVolume(D_sinoData, sinoPitch, dims.iProjDets, dims.iProjAngles);
-	printf("pitch: %u\n", sinoPitch);
-	
-	unsigned int y, x;
-	float* sino = loadImage("sino.png", y, x);
-
-	float* img = new float[dims.iVolWidth*dims.iVolHeight];
-
-	copySinogramToDevice(sino, dims.iProjDets, dims.iProjDets, dims.iProjAngles, D_sinoData, sinoPitch);
-
-	float* angle = new float[dims.iProjAngles];
-
-	for (unsigned int i = 0; i < dims.iProjAngles; ++i)
-		angle[i] = i*(M_PI/dims.iProjAngles);
-
-	SIRT sirt;
-
-	sirt.setGeometry(dims, angle);
-	sirt.init();
-
-	sirt.setBuffers(D_volumeData, volumePitch, D_sinoData, sinoPitch);
-
-	sirt.iterate(25);
-
-
-	delete[] angle;
-
-	copyVolumeFromDevice(img, dims.iVolWidth, dims, D_volumeData, volumePitch);
-
-	saveImage("vol.png",dims.iVolHeight,dims.iVolWidth,img);
-
-	return 0;
-}
-#endif
-
diff --git a/cuda/3d/cgls3d.cu b/cuda/3d/cgls3d.cu
index 10c5f1e..4829574 100644
--- a/cuda/3d/cgls3d.cu
+++ b/cuda/3d/cgls3d.cu
@@ -33,10 +33,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include <cstdio>
 #include <cassert>
 
-#ifdef STANDALONE
-#include "testutil.h"
-#endif
-
 namespace astraCUDA3d {
 
 CGLS::CGLS() : ReconAlgo3D()
@@ -263,161 +259,3 @@ bool doCGLS(cudaPitchedPtr& D_volumeData,
 }
 
 }
-
-#ifdef STANDALONE
-
-using namespace astraCUDA3d;
-
-int main()
-{
-	SDimensions3D dims;
-	dims.iVolX = 256;
-	dims.iVolY = 256;
-	dims.iVolZ = 256;
-	dims.iProjAngles = 100;
-	dims.iProjU = 512;
-	dims.iProjV = 512;
-	dims.iRaysPerDet = 1;
-
-	SConeProjection angle[100];
-	angle[0].fSrcX = -2905.6;
-	angle[0].fSrcY = 0;
-	angle[0].fSrcZ = 0;
-
-	angle[0].fDetSX = 694.4;
-	angle[0].fDetSY = -122.4704;
-	angle[0].fDetSZ = -122.4704;
-
-	angle[0].fDetUX = 0;
-	angle[0].fDetUY = .4784;
-	//angle[0].fDetUY = .5;
-	angle[0].fDetUZ = 0;
-
-	angle[0].fDetVX = 0;
-	angle[0].fDetVY = 0;
-	angle[0].fDetVZ = .4784;
-
-#define ROTATE0(name,i,alpha) do { angle[i].f##name##X = angle[0].f##name##X * cos(alpha) - angle[0].f##name##Y * sin(alpha); angle[i].f##name##Y = angle[0].f##name##X * sin(alpha) + angle[0].f##name##Y * cos(alpha); } while(0)
-	for (int i = 1; i < 100; ++i) {
-		angle[i] = angle[0];
-		ROTATE0(Src, i, i*2*M_PI/100);
-		ROTATE0(DetS, i, i*2*M_PI/100);
-		ROTATE0(DetU, i, i*2*M_PI/100);
-		ROTATE0(DetV, i, i*2*M_PI/100);
-	}
-#undef ROTATE0
-
-
-	cudaPitchedPtr volData = allocateVolumeData(dims);
-	cudaPitchedPtr projData = allocateProjectionData(dims);
-	zeroProjectionData(projData, dims);
-
-	float* pbuf = new float[100*512*512];
-	copyProjectionsFromDevice(pbuf, projData, dims);
-	copyProjectionsToDevice(pbuf, projData, dims);
-	delete[] pbuf;
-
-#if 0
-	float* slice = new float[256*256];
-	cudaPitchedPtr ptr;
-	ptr.ptr = slice;
-	ptr.pitch = 256*sizeof(float);
-	ptr.xsize = 256*sizeof(float);
-	ptr.ysize = 256;
-
-	for (unsigned int i = 0; i < 256; ++i) {
-		for (unsigned int y = 0; y < 256; ++y)
-			for (unsigned int x = 0; x < 256; ++x)
-				slice[y*256+x] = (i-127.5)*(i-127.5)+(y-127.5)*(y-127.5)+(x-127.5)*(x-127.5) < 4900 ? 1.0f : 0.0f;
-
-		cudaExtent extentS;
-		extentS.width = dims.iVolX*sizeof(float);
-		extentS.height = dims.iVolY;
-		extentS.depth = 1;
-		cudaPos sp = { 0, 0, 0 };
-		cudaPos dp = { 0, 0, i };
-		cudaMemcpy3DParms p;
-		p.srcArray = 0;
-		p.srcPos = sp;
-		p.srcPtr = ptr;
-		p.dstArray = 0;
-		p.dstPos = dp;
-		p.dstPtr = volData;
-		p.extent = extentS;
-		p.kind = cudaMemcpyHostToDevice;
-		cudaMemcpy3D(&p);
-	}
-	astraCUDA3d::ConeFP(volData, projData, dims, angle, 1.0f);
-
-#else
-
-	for (int i = 0; i < 100; ++i) {
-		char fname[32];
-		sprintf(fname, "Tiffs/%04d.png", 4*i);
-		unsigned int w,h;
-		float* bufp = loadImage(fname, w,h);
-
-		for (int j = 0; j < 512*512; ++j) {
-			float v = bufp[j];
-			if (v > 236.0f) v = 236.0f;
-			v = logf(236.0f / v);
-			bufp[j] = 256*v;
-		}
-
-		for (int j = 0; j < 512; ++j) {
-			cudaMemcpy(((float*)projData.ptr)+100*512*j+512*i, bufp+512*j, 512*sizeof(float), cudaMemcpyHostToDevice);
-		}
-
-		delete[] bufp;
-
-	}
-#endif
-
-#if 0
-	float* bufs = new float[100*512];
-
-	for (int i = 0; i < 512; ++i) {
-		cudaMemcpy(bufs, ((float*)projData.ptr)+100*512*i, 100*512*sizeof(float), cudaMemcpyDeviceToHost);
-
-		printf("%d %d %d\n", projData.pitch, projData.xsize, projData.ysize);
-
-		char fname[20];
-		sprintf(fname, "sino%03d.png", i);
-		saveImage(fname, 100, 512, bufs);
-	}
-
-	float* bufp = new float[512*512];
-
-	for (int i = 0; i < 100; ++i) {
-		for (int j = 0; j < 512; ++j) {
-			cudaMemcpy(bufp+512*j, ((float*)projData.ptr)+100*512*j+512*i, 512*sizeof(float), cudaMemcpyDeviceToHost);
-		}
-
-		char fname[20];
-		sprintf(fname, "proj%03d.png", i);
-		saveImage(fname, 512, 512, bufp);
-	}
-#endif
-
-	zeroVolumeData(volData, dims);
-
-	cudaPitchedPtr maskData;
-	maskData.ptr = 0;
-
-	astraCUDA3d::doCGLS(volData, projData, maskData, dims, angle, 50);
-#if 1
-	float* buf = new float[256*256];
-
-	for (int i = 0; i < 256; ++i) {
-		cudaMemcpy(buf, ((float*)volData.ptr)+256*256*i, 256*256*sizeof(float), cudaMemcpyDeviceToHost);
-
-		char fname[20];
-		sprintf(fname, "vol%03d.png", i);
-		saveImage(fname, 256, 256, buf);
-	}
-#endif
-
-	return 0;
-}
-#endif
-
diff --git a/cuda/3d/cone_bp.cu b/cuda/3d/cone_bp.cu
index feebda2..7312bbc 100644
--- a/cuda/3d/cone_bp.cu
+++ b/cuda/3d/cone_bp.cu
@@ -28,11 +28,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/3d/util3d.h"
 #include "astra/cuda/3d/dims3d.h"
 
-#ifdef STANDALONE
-#include "astra/cuda/3d/cone_fp.h"
-#include "testutil.h"
-#endif
-
 #include <cstdio>
 #include <cassert>
 #include <iostream>
@@ -55,7 +50,13 @@ static const unsigned int g_volBlockY = 32;
 
 static const unsigned g_MaxAngles = 1024;
 
-__constant__ float gC_C[12*g_MaxAngles];
+struct DevConeParams {
+	float4 fNumU;
+	float4 fNumV;
+	float4 fDen;
+};
+
+__constant__ DevConeParams gC_C[g_MaxAngles];
 
 bool bindProjDataTexture(const cudaArray* array)
 {
@@ -118,16 +119,13 @@ __global__ void dev_cone_BP(void* D_volData, unsigned int volPitch, int startAng
 
 		for (int angle = startAngle; angle < endAngle; ++angle, fAngle += 1.0f)
 		{
-			float4 fCu  = make_float4(gC_C[12*angle+0], gC_C[12*angle+1], gC_C[12*angle+2], gC_C[12*angle+3]);
-			float4 fCv  = make_float4(gC_C[12*angle+4], gC_C[12*angle+5], gC_C[12*angle+6], gC_C[12*angle+7]);
-			float4 fCd  = make_float4(gC_C[12*angle+8], gC_C[12*angle+9], gC_C[12*angle+10], gC_C[12*angle+11]);
+			float4 fCu  = gC_C[angle].fNumU;
+			float4 fCv  = gC_C[angle].fNumV;
+			float4 fCd  = gC_C[angle].fDen;
 
 			float fUNum = fCu.w + fX * fCu.x + fY * fCu.y + fZ * fCu.z;
 			float fVNum = fCv.w + fX * fCv.x + fY * fCv.y + fZ * fCv.z;
-			float fDen  = fCd.w + fX * fCd.x + fY * fCd.y + fZ * fCd.z;
-
-			// fCd.w = -|| u v s || (determinant of 3x3 matrix with cols u,v,s)
-			// fDen =  || u v (x-s) ||
+			float fDen  = (FDKWEIGHT ? 1.0f : fCd.w) + fX * fCd.x + fY * fCd.y + fZ * fCd.z;
 
 			float fU,fV, fr;
 
@@ -137,18 +135,7 @@ __global__ void dev_cone_BP(void* D_volData, unsigned int volPitch, int startAng
 				fU = fUNum * fr;
 				fV = fVNum * fr;
 				float fVal = tex3D(gT_coneProjTexture, fU, fAngle, fV);
-				if (FDKWEIGHT) {
-					// The correct factor here is this one:
-					// Z[idx] += (fr*fCd.w)*(fr*fCd.w)*fVal;
-					// This is the square of the inverse magnification factor
-					// from fX,fY,fZ to the detector.
-
-					// Since we are assuming we have a circular cone
-					// beam trajectory, fCd.w is constant, and we instead
-					// multiply by fCd.w*fCd.w in the FDK preweighting step.
-					Z[idx] += fr*fr*fVal;
-				} else
-					Z[idx] += fVal;
+				Z[idx] += fr*fr*fVal;
 
 				fUNum += fCu.z;
 				fVNum += fCv.z;
@@ -215,19 +202,9 @@ __global__ void dev_cone_BP_SS(void* D_volData, unsigned int volPitch, int start
 
 		for (int angle = startAngle; angle < endAngle; ++angle, fAngle += 1.0f)
 		{
-
-			const float fCux = gC_C[12*angle+0];
-			const float fCuy = gC_C[12*angle+1];
-			const float fCuz = gC_C[12*angle+2];
-			const float fCuc = gC_C[12*angle+3];
-			const float fCvx = gC_C[12*angle+4];
-			const float fCvy = gC_C[12*angle+5];
-			const float fCvz = gC_C[12*angle+6];
-			const float fCvc = gC_C[12*angle+7];
-			const float fCdx = gC_C[12*angle+8];
-			const float fCdy = gC_C[12*angle+9];
-			const float fCdz = gC_C[12*angle+10];
-			const float fCdc = gC_C[12*angle+11];
+			float4 fCu  = gC_C[angle].fNumU;
+			float4 fCv  = gC_C[angle].fNumV;
+			float4 fCd  = gC_C[angle].fDen;
 
 			float fXs = fX;
 			for (int iSubX = 0; iSubX < iRaysPerVoxelDim; ++iSubX) {
@@ -236,14 +213,15 @@ __global__ void dev_cone_BP_SS(void* D_volData, unsigned int volPitch, int start
 			float fZs = fZ;
 			for (int iSubZ = 0; iSubZ < iRaysPerVoxelDim; ++iSubZ) {
 
-				const float fUNum = fCuc + fXs * fCux + fYs * fCuy + fZs * fCuz;
-				const float fVNum = fCvc + fXs * fCvx + fYs * fCvy + fZs * fCvz;
-				const float fDen = fCdc + fXs * fCdx + fYs * fCdy + fZs * fCdz;
+				const float fUNum = fCu.w + fX * fCu.x + fY * fCu.y + fZ * fCu.z;
+				const float fVNum = fCv.w + fX * fCv.x + fY * fCv.y + fZ * fCv.z;
+				const float fDen  = fCd.w + fX * fCd.x + fY * fCd.y + fZ * fCd.z;
 
-				const float fU = fUNum / fDen;
-				const float fV = fVNum / fDen;
+				const float fr = __fdividef(1.0f, fDen);
+				const float fU = fUNum * fr;
+				const float fV = fVNum * fr;
 
-				fVal += tex3D(gT_coneProjTexture, fU, fV, fAngle);
+				fVal += tex3D(gT_coneProjTexture, fU, fV, fAngle) * fr;
 
 				fZs += fSubStep;
 			}
@@ -259,6 +237,76 @@ __global__ void dev_cone_BP_SS(void* D_volData, unsigned int volPitch, int start
 }
 
 
+bool transferConstants(const SConeProjection* angles, unsigned int iProjAngles, const SProjectorParams3D& params)
+{
+	DevConeParams *p = new DevConeParams[iProjAngles];
+
+	// We need three things in the kernel:
+	// projected coordinates of pixels on the detector:
+
+	// u: || (x-s) v (s-d) || / || u v (s-x) ||
+	// v: -|| u (x-s) (s-d) || / || u v (s-x) ||
+
+	// ray density weighting factor for the adjoint
+	// || u v (s-d) ||^2 / ( |cross(u,v)| * || u v (s-x) ||^2 )
+
+	// FDK weighting factor
+	// ( || u v s || / || u v (s-x) || ) ^ 2
+
+	// Since u and v are ratios with the same denominator, we have
+	// a degree of freedom to scale the denominator. We use that to make
+	// the square of the denominator equal to the relevant weighting factor.
+
+
+	for (unsigned int i = 0; i < iProjAngles; ++i) {
+		Vec3 u(angles[i].fDetUX, angles[i].fDetUY, angles[i].fDetUZ);
+		Vec3 v(angles[i].fDetVX, angles[i].fDetVY, angles[i].fDetVZ);
+		Vec3 s(angles[i].fSrcX, angles[i].fSrcY, angles[i].fSrcZ);
+		Vec3 d(angles[i].fDetSX, angles[i].fDetSY, angles[i].fDetSZ);
+
+
+
+		double fScale;
+		if (!params.bFDKWeighting) {
+			// goal: 1/fDen^2 = || u v (s-d) ||^2 / ( |cross(u,v)| * || u v (s-x) ||^2 )
+			// fDen = ( sqrt(|cross(u,v)|) * || u v (s-x) || ) / || u v (s-d) || 
+			// i.e. scale = sqrt(|cross(u,v)|) * / || u v (s-d) ||
+
+
+			// NB: for cross(u,v) we invert the volume scaling (for the voxel
+			// size normalization) to get the proper dimensions for
+			// the scaling of the adjoint
+
+			fScale = sqrt(scaled_cross3(u,v,Vec3(params.fVolScaleX,params.fVolScaleY,params.fVolScaleZ)).norm()) / det3(u, v, s-d);
+		} else {
+			// goal: 1/fDen = || u v s || / || u v (s-x) ||
+			// fDen = || u v (s-x) || / || u v s ||
+			// i.e., scale = 1 / || u v s ||
+
+			fScale = 1.0 / det3(u, v, s);
+		}
+
+		p[i].fNumU.w = fScale * det3(s,v,d);
+		p[i].fNumU.x = fScale * det3x(v,s-d);
+		p[i].fNumU.y = fScale * det3y(v,s-d);
+		p[i].fNumU.z = fScale * det3z(v,s-d);
+		p[i].fNumV.w = -fScale * det3(s,u,d);
+		p[i].fNumV.x = -fScale * det3x(u,s-d);
+		p[i].fNumV.y = -fScale * det3y(u,s-d);
+		p[i].fNumV.z = -fScale * det3z(u,s-d);
+		p[i].fDen.w = fScale * det3(u, v, s); // == 1.0 for FDK
+		p[i].fDen.x = -fScale * det3x(u, v);
+		p[i].fDen.y = -fScale * det3y(u, v);
+		p[i].fDen.z = -fScale * det3z(u, v);
+	}
+
+	// TODO: Check for errors
+	cudaMemcpyToSymbol(gC_C, p, iProjAngles*sizeof(DevConeParams), 0, cudaMemcpyHostToDevice);
+
+	return true;
+}
+
+
 bool ConeBP_Array(cudaPitchedPtr D_volumeData,
                   cudaArray *D_projArray,
                   const SDimensions3D& dims, const SConeProjection* angles,
@@ -267,44 +315,21 @@ bool ConeBP_Array(cudaPitchedPtr D_volumeData,
 	bindProjDataTexture(D_projArray);
 
 	float fOutputScale;
-	if (params.bFDKWeighting)
-		fOutputScale = params.fOutputScale / (params.fVolScaleX * params.fVolScaleY * params.fVolScaleZ);
-	else
+	if (params.bFDKWeighting) {
+		// NB: assuming cube voxels here
+		fOutputScale = params.fOutputScale / (params.fVolScaleX);
+	} else {
 		fOutputScale = params.fOutputScale * (params.fVolScaleX * params.fVolScaleY * params.fVolScaleZ);
+	}
 
 	for (unsigned int th = 0; th < dims.iProjAngles; th += g_MaxAngles) {
 		unsigned int angleCount = g_MaxAngles;
 		if (th + angleCount > dims.iProjAngles)
 			angleCount = dims.iProjAngles - th;
 
-		// transfer angles to constant memory
-		float* tmp = new float[12*angleCount];
-
-
-		// NB: We increment angles at the end of the loop body.
-
-
-#define TRANSFER_TO_CONSTANT(expr,name) do { for (unsigned int i = 0; i < angleCount; ++i) tmp[12*i+name] = (expr) ; } while (0)
-
-		TRANSFER_TO_CONSTANT( (angles[i].fDetSZ - angles[i].fSrcZ)*angles[i].fDetVY - (angles[i].fDetSY - angles[i].fSrcY)*angles[i].fDetVZ , 0 );
-		TRANSFER_TO_CONSTANT( (angles[i].fDetSX - angles[i].fSrcX)*angles[i].fDetVZ -(angles[i].fDetSZ - angles[i].fSrcZ)*angles[i].fDetVX , 1 );
-		TRANSFER_TO_CONSTANT( (angles[i].fDetSY - angles[i].fSrcY)*angles[i].fDetVX - (angles[i].fDetSX - angles[i].fSrcX)*angles[i].fDetVY , 2 );
-		TRANSFER_TO_CONSTANT( (angles[i].fDetSY*angles[i].fDetVZ - angles[i].fDetSZ*angles[i].fDetVY)*angles[i].fSrcX - (angles[i].fDetSX*angles[i].fDetVZ - angles[i].fDetSZ*angles[i].fDetVX)*angles[i].fSrcY + (angles[i].fDetSX*angles[i].fDetVY - angles[i].fDetSY*angles[i].fDetVX)*angles[i].fSrcZ , 3 );
-
-		TRANSFER_TO_CONSTANT( (angles[i].fDetSY - angles[i].fSrcY)*angles[i].fDetUZ-(angles[i].fDetSZ - angles[i].fSrcZ)*angles[i].fDetUY, 4 );
-		TRANSFER_TO_CONSTANT( (angles[i].fDetSZ - angles[i].fSrcZ)*angles[i].fDetUX - (angles[i].fDetSX - angles[i].fSrcX)*angles[i].fDetUZ , 5 );
-		TRANSFER_TO_CONSTANT((angles[i].fDetSX - angles[i].fSrcX)*angles[i].fDetUY-(angles[i].fDetSY - angles[i].fSrcY)*angles[i].fDetUX , 6 );
-		TRANSFER_TO_CONSTANT( -(angles[i].fDetSY*angles[i].fDetUZ - angles[i].fDetSZ*angles[i].fDetUY)*angles[i].fSrcX + (angles[i].fDetSX*angles[i].fDetUZ - angles[i].fDetSZ*angles[i].fDetUX)*angles[i].fSrcY - (angles[i].fDetSX*angles[i].fDetUY - angles[i].fDetSY*angles[i].fDetUX)*angles[i].fSrcZ , 7 );
-
-		TRANSFER_TO_CONSTANT( angles[i].fDetUY*angles[i].fDetVZ - angles[i].fDetUZ*angles[i].fDetVY , 8 );
-		TRANSFER_TO_CONSTANT( angles[i].fDetUZ*angles[i].fDetVX - angles[i].fDetUX*angles[i].fDetVZ , 9 );
-		TRANSFER_TO_CONSTANT( angles[i].fDetUX*angles[i].fDetVY - angles[i].fDetUY*angles[i].fDetVX , 10 );
-		TRANSFER_TO_CONSTANT( -angles[i].fSrcX * (angles[i].fDetUY*angles[i].fDetVZ - angles[i].fDetUZ*angles[i].fDetVY) - angles[i].fSrcY * (angles[i].fDetUZ*angles[i].fDetVX - angles[i].fDetUX*angles[i].fDetVZ) - angles[i].fSrcZ * (angles[i].fDetUX*angles[i].fDetVY - angles[i].fDetUY*angles[i].fDetVX) , 11 );
-
-#undef TRANSFER_TO_CONSTANT
-		cudaMemcpyToSymbol(gC_C, tmp, angleCount*12*sizeof(float), 0, cudaMemcpyHostToDevice); 
-
-		delete[] tmp;
+		bool ok = transferConstants(angles, angleCount, params);
+		if (!ok)
+			return false;
 
 		dim3 dimBlock(g_volBlockX, g_volBlockY);
 
@@ -353,168 +378,3 @@ bool ConeBP(cudaPitchedPtr D_volumeData,
 
 
 }
-
-#ifdef STANDALONE
-int main()
-{
-	astraCUDA3d::SDimensions3D dims;
-	dims.iVolX = 512;
-	dims.iVolY = 512;
-	dims.iVolZ = 512;
-	dims.iProjAngles = 496;
-	dims.iProjU = 512;
-	dims.iProjV = 512;
-	dims.iRaysPerDetDim = 1;
-	dims.iRaysPerVoxelDim = 1;
-
-	cudaExtent extentV;
-	extentV.width = dims.iVolX*sizeof(float);
-	extentV.height = dims.iVolY;
-	extentV.depth = dims.iVolZ;
-
-	cudaPitchedPtr volData; // pitch, ptr, xsize, ysize
-
-	cudaMalloc3D(&volData, extentV);
-
-	cudaExtent extentP;
-	extentP.width = dims.iProjU*sizeof(float);
-	extentP.height = dims.iProjAngles;
-	extentP.depth = dims.iProjV;
-
-	cudaPitchedPtr projData; // pitch, ptr, xsize, ysize
-
-	cudaMalloc3D(&projData, extentP);
-	cudaMemset3D(projData, 0, extentP);
-
-#if 0
-	float* slice = new float[256*256];
-	cudaPitchedPtr ptr;
-	ptr.ptr = slice;
-	ptr.pitch = 256*sizeof(float);
-	ptr.xsize = 256*sizeof(float);
-	ptr.ysize = 256;
-
-	for (unsigned int i = 0; i < 256*256; ++i)
-		slice[i] = 1.0f;
-	for (unsigned int i = 0; i < 256; ++i) {
-		cudaExtent extentS;
-		extentS.width = dims.iVolX*sizeof(float);
-		extentS.height = dims.iVolY;
-		extentS.depth = 1;
-		cudaPos sp = { 0, 0, 0 };
-		cudaPos dp = { 0, 0, i };
-		cudaMemcpy3DParms p;
-		p.srcArray = 0;
-		p.srcPos = sp;
-		p.srcPtr = ptr;
-		p.dstArray = 0;
-		p.dstPos = dp;
-		p.dstPtr = volData;
-		p.extent = extentS;
-		p.kind = cudaMemcpyHostToDevice;
-		cudaMemcpy3D(&p);
-#if 0
-		if (i == 128) {
-			for (unsigned int j = 0; j < 256*256; ++j)
-				slice[j] = 0.0f;
-		}
-#endif 
-	}
-#endif
-
-
-	astraCUDA3d::SConeProjection angle[512];
-	angle[0].fSrcX = -5120;
-	angle[0].fSrcY = 0;
-	angle[0].fSrcZ = 0;
-
-	angle[0].fDetSX = 512;
-	angle[0].fDetSY = -256;
-	angle[0].fDetSZ = -256;
-
-	angle[0].fDetUX = 0;
-	angle[0].fDetUY = 1;
-	angle[0].fDetUZ = 0;
-
-	angle[0].fDetVX = 0;
-	angle[0].fDetVY = 0;
-	angle[0].fDetVZ = 1;
-
-#define ROTATE0(name,i,alpha) do { angle[i].f##name##X = angle[0].f##name##X * cos(alpha) - angle[0].f##name##Y * sin(alpha); angle[i].f##name##Y = angle[0].f##name##X * sin(alpha) + angle[0].f##name##Y * cos(alpha); } while(0)
-	for (int i = 1; i < 512; ++i) {
-		angle[i] = angle[0];
-		ROTATE0(Src, i, i*2*M_PI/512);
-		ROTATE0(DetS, i, i*2*M_PI/512);
-		ROTATE0(DetU, i, i*2*M_PI/512);
-		ROTATE0(DetV, i, i*2*M_PI/512);
-	}
-#undef ROTATE0
-
-#if 0
-	astraCUDA3d::ConeFP(volData, projData, dims, angle, 1.0f);
-#endif
-#if 0
-	float* bufs = new float[180*512];
-
-	for (int i = 0; i < 512; ++i) {
-		cudaMemcpy(bufs, ((float*)projData.ptr)+180*512*i, 180*512*sizeof(float), cudaMemcpyDeviceToHost);
-
-		printf("%d %d %d\n", projData.pitch, projData.xsize, projData.ysize);
-
-		char fname[20];
-		sprintf(fname, "sino%03d.png", i);
-		saveImage(fname, 180, 512, bufs);
-	}
-
-	float* bufp = new float[512*512];
-
-	for (int i = 0; i < 180; ++i) {
-		for (int j = 0; j < 512; ++j) {
-			cudaMemcpy(bufp+512*j, ((float*)projData.ptr)+180*512*j+512*i, 512*sizeof(float), cudaMemcpyDeviceToHost);
-		}
-
-		char fname[20];
-		sprintf(fname, "proj%03d.png", i);
-		saveImage(fname, 512, 512, bufp);
-	}
-#endif		
-#if 0
-	for (unsigned int i = 0; i < 256*256; ++i)
-		slice[i] = 0.0f;
-	for (unsigned int i = 0; i < 256; ++i) {
-		cudaExtent extentS;
-		extentS.width = dims.iVolX*sizeof(float);
-		extentS.height = dims.iVolY;
-		extentS.depth = 1;
-		cudaPos sp = { 0, 0, 0 };
-		cudaPos dp = { 0, 0, i };
-		cudaMemcpy3DParms p;
-		p.srcArray = 0;
-		p.srcPos = sp;
-		p.srcPtr = ptr;
-		p.dstArray = 0;
-		p.dstPos = dp;
-		p.dstPtr = volData;
-		p.extent = extentS;
-		p.kind = cudaMemcpyHostToDevice;
-		cudaMemcpy3D(&p);
-	}
-#endif
-
-	astraCUDA3d::ConeBP(volData, projData, dims, angle, 1.0f);
-#if 0
-	float* buf = new float[256*256];
-
-	for (int i = 0; i < 256; ++i) {
-		cudaMemcpy(buf, ((float*)volData.ptr)+256*256*i, 256*256*sizeof(float), cudaMemcpyDeviceToHost);
-
-		printf("%d %d %d\n", volData.pitch, volData.xsize, volData.ysize);
-
-		char fname[20];
-		sprintf(fname, "vol%03d.png", i);
-		saveImage(fname, 256, 256, buf);
-	}
-#endif
-
-}
-#endif
diff --git a/cuda/3d/cone_fp.cu b/cuda/3d/cone_fp.cu
index 7e0fae8..bd607fa 100644
--- a/cuda/3d/cone_fp.cu
+++ b/cuda/3d/cone_fp.cu
@@ -28,10 +28,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/3d/util3d.h"
 #include "astra/cuda/3d/dims3d.h"
 
-#ifdef STANDALONE
-#include "testutil.h"
-#endif
-
 #include <cstdio>
 #include <cassert>
 #include <iostream>
@@ -368,7 +364,7 @@ bool ConeFP_Array_internal(cudaPitchedPtr D_projData,
 	SCALE_NONCUBE snoncubeY;
 	fS1 = params.fVolScaleX / params.fVolScaleY;
 	snoncubeY.fScale1 = fS1 * fS1;
-	fS2 = params.fVolScaleY / params.fVolScaleY;
+	fS2 = params.fVolScaleZ / params.fVolScaleY;
 	snoncubeY.fScale2 = fS2 * fS2;
 	snoncubeY.fOutputScale = params.fOutputScale * params.fVolScaleY;
 
@@ -498,105 +494,3 @@ bool ConeFP(cudaPitchedPtr D_volumeData,
 
 
 }
-
-#ifdef STANDALONE
-int main()
-{
-	SDimensions3D dims;
-	dims.iVolX = 256;
-	dims.iVolY = 256;
-	dims.iVolZ = 256;
-	dims.iProjAngles = 32;
-	dims.iProjU = 512;
-	dims.iProjV = 512;
-	dims.iRaysPerDet = 1;
-
-	cudaExtent extentV;
-	extentV.width = dims.iVolX*sizeof(float);
-	extentV.height = dims.iVolY;
-	extentV.depth = dims.iVolZ;
-
-	cudaPitchedPtr volData; // pitch, ptr, xsize, ysize
-
-	cudaMalloc3D(&volData, extentV);
-
-	cudaExtent extentP;
-	extentP.width = dims.iProjU*sizeof(float);
-	extentP.height = dims.iProjV;
-	extentP.depth = dims.iProjAngles;
-
-	cudaPitchedPtr projData; // pitch, ptr, xsize, ysize
-
-	cudaMalloc3D(&projData, extentP);
-	cudaMemset3D(projData, 0, extentP);
-
-	float* slice = new float[256*256];
-	cudaPitchedPtr ptr;
-	ptr.ptr = slice;
-	ptr.pitch = 256*sizeof(float);
-	ptr.xsize = 256*sizeof(float);
-	ptr.ysize = 256;
-
-	for (unsigned int i = 0; i < 256*256; ++i)
-		slice[i] = 1.0f;
-	for (unsigned int i = 0; i < 256; ++i) {
-		cudaExtent extentS;
-		extentS.width = dims.iVolX*sizeof(float);
-		extentS.height = dims.iVolY;
-		extentS.depth = 1;
-		cudaPos sp = { 0, 0, 0 };
-		cudaPos dp = { 0, 0, i };
-		cudaMemcpy3DParms p;
-		p.srcArray = 0;
-		p.srcPos = sp;
-		p.srcPtr = ptr;
-		p.dstArray = 0;
-		p.dstPos = dp;
-		p.dstPtr = volData;
-		p.extent = extentS;
-		p.kind = cudaMemcpyHostToDevice;
-		cudaError err = cudaMemcpy3D(&p);
-		assert(!err);
-	}
-
-
-	SConeProjection angle[32];
-	angle[0].fSrcX = -1536;
-	angle[0].fSrcY = 0;
-	angle[0].fSrcZ = 200;
-
-	angle[0].fDetSX = 512;
-	angle[0].fDetSY = -256;
-	angle[0].fDetSZ = -256;
-
-	angle[0].fDetUX = 0;
-	angle[0].fDetUY = 1;
-	angle[0].fDetUZ = 0;
-
-	angle[0].fDetVX = 0;
-	angle[0].fDetVY = 0;
-	angle[0].fDetVZ = 1;
-
-#define ROTATE0(name,i,alpha) do { angle[i].f##name##X = angle[0].f##name##X * cos(alpha) - angle[0].f##name##Y * sin(alpha); angle[i].f##name##Y = angle[0].f##name##X * sin(alpha) + angle[0].f##name##Y * cos(alpha); } while(0)
-	for (int i = 1; i < 32; ++i) {
-		angle[i] = angle[0];
-		ROTATE0(Src, i, i*1*M_PI/180);
-		ROTATE0(DetS, i, i*1*M_PI/180);
-		ROTATE0(DetU, i, i*1*M_PI/180);
-		ROTATE0(DetV, i, i*1*M_PI/180);
-	}
-#undef ROTATE0
-
-	astraCUDA3d::ConeFP(volData, projData, dims, angle, 1.0f);
-
-	float* buf = new float[512*512];
-
-	cudaMemcpy(buf, ((float*)projData.ptr)+512*512*8, 512*512*sizeof(float), cudaMemcpyDeviceToHost);
-
-	printf("%d %d %d\n", projData.pitch, projData.xsize, projData.ysize);
-
-	saveImage("proj.png", 512, 512, buf);
-	
-
-}
-#endif
diff --git a/cuda/3d/fdk.cu b/cuda/3d/fdk.cu
index 1294721..456694f 100644
--- a/cuda/3d/fdk.cu
+++ b/cuda/3d/fdk.cu
@@ -32,11 +32,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 
 #include "astra/cuda/2d/fft.h"
 
-#ifdef STANDALONE
-#include "astra/cuda/3d/cone_fp.h"
-#include "testutil.h"
-#endif
-
 #include "astra/Logging.h"
 
 #include <cstdio>
@@ -57,10 +52,13 @@ static const unsigned g_MaxAngles = 12000;
 __constant__ float gC_angle[g_MaxAngles];
 
 
-// per-detector u/v shifts?
 
+// TODO: To support non-cube voxels, preweighting needs per-view
+// parameters. NB: Need to properly take into account the
+// anisotropic volume normalization done for that too.
 
-__global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsigned int startAngle, unsigned int endAngle, float fSrcOrigin, float fDetOrigin, float fZShift, float fDetUSize, float fDetVSize, float fVoxSize, const SDimensions3D dims)
+
+__global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsigned int startAngle, unsigned int endAngle, float fSrcOrigin, float fDetOrigin, float fZShift, float fDetUSize, float fDetVSize, const SDimensions3D dims)
 {
 	float* projData = (float*)D_projData;
 	int angle = startAngle + blockIdx.y * g_anglesPerWeightBlock + threadIdx.y;
@@ -88,14 +86,10 @@ __global__ void devFDK_preweight(void* D_projData, unsigned int projPitch, unsig
 	// fCentralRayLength / fRayLength   : the main FDK preweighting factor
 	// fSrcOrigin / (fDetUSize * fCentralRayLength)
 	//                                  : to adjust the filter to the det width
-	// || u v s || ^ 2                  : see cone_bp.cu, FDKWEIGHT
 	// pi / (2 * iProjAngles)           : scaling of the integral over angles
-	// fVoxSize ^ 2                     : ...
 
-	const float fW1 = fSrcOrigin * fDetUSize * fDetVSize;
 	const float fW2 = fCentralRayLength / (fDetUSize * fSrcOrigin);
-	const float fW3 = fVoxSize * fVoxSize;
-	const float fW = fCentralRayLength * fW1 * fW1 * fW2 * fW3 * (M_PI / 2.0f) / (float)dims.iProjAngles;
+	const float fW = fCentralRayLength * fW2 * (M_PI / 2.0f) / (float)dims.iProjAngles;
 
 	for (int detectorV = startDetectorV; detectorV < endDetectorV; ++detectorV)
 	{
@@ -167,7 +161,7 @@ __global__ void devFDK_ParkerWeight(void* D_projData, unsigned int projPitch, un
 bool FDK_PreWeight(cudaPitchedPtr D_projData,
                 float fSrcOrigin, float fDetOrigin,
                 float fZShift,
-                float fDetUSize, float fDetVSize, float fVoxSize,
+                float fDetUSize, float fDetVSize,
 				bool bShortScan,
                 const SDimensions3D& dims, const float* angles)
 {
@@ -180,7 +174,7 @@ bool FDK_PreWeight(cudaPitchedPtr D_projData,
 
 	int projPitch = D_projData.pitch/sizeof(float);
 
-	devFDK_preweight<<<dimGrid, dimBlock>>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fZShift, fDetUSize, fDetVSize, fVoxSize, dims);
+	devFDK_preweight<<<dimGrid, dimBlock>>>(D_projData.ptr, projPitch, 0, dims.iProjAngles, fSrcOrigin, fDetOrigin, fZShift, fDetUSize, fDetVSize, dims);
 
 	cudaTextForceKernelsCompletion();
 
@@ -344,9 +338,8 @@ bool FDK(cudaPitchedPtr D_volumeData,
 
 
 #if 1
-	// NB: assuming cube voxels (params.fVolScaleX)
 	ok = FDK_PreWeight(D_projData, fSrcOrigin, fDetOrigin,
-	                fZShift, fDetUSize, fDetVSize, params.fVolScaleX,
+	                fZShift, fDetUSize, fDetVSize,
 	                bShortScan, dims, pfAngles);
 #else
 	ok = true;
@@ -379,220 +372,3 @@ bool FDK(cudaPitchedPtr D_volumeData,
 
 
 }
-
-#ifdef STANDALONE
-void dumpVolume(const char* filespec, const cudaPitchedPtr& data, const SDimensions3D& dims, float fMin, float fMax)
-{
-	float* buf = new float[dims.iVolX*dims.iVolY];
-	unsigned int pitch = data.pitch / sizeof(float);
-
-	for (int i = 0; i < dims.iVolZ; ++i) {
-		cudaMemcpy2D(buf, dims.iVolX*sizeof(float), ((float*)data.ptr)+pitch*dims.iVolY*i, data.pitch, dims.iVolX*sizeof(float), dims.iVolY, cudaMemcpyDeviceToHost);
-
-		char fname[512];
-		sprintf(fname, filespec, dims.iVolZ-i-1);
-		saveImage(fname, dims.iVolY, dims.iVolX, buf, fMin, fMax);
-	}
-}
-
-void dumpSinograms(const char* filespec, const cudaPitchedPtr& data, const SDimensions3D& dims, float fMin, float fMax)
-{
-	float* bufs = new float[dims.iProjAngles*dims.iProjU];
-	unsigned int pitch = data.pitch / sizeof(float);
-
-	for (int i = 0; i < dims.iProjV; ++i) {
-		cudaMemcpy2D(bufs, dims.iProjU*sizeof(float), ((float*)data.ptr)+pitch*dims.iProjAngles*i, data.pitch, dims.iProjU*sizeof(float), dims.iProjAngles, cudaMemcpyDeviceToHost);
-
-		char fname[512];
-		sprintf(fname, filespec, i);
-		saveImage(fname, dims.iProjAngles, dims.iProjU, bufs, fMin, fMax);
-	}
-}
-
-void dumpProjections(const char* filespec, const cudaPitchedPtr& data, const SDimensions3D& dims, float fMin, float fMax)
-{
-	float* bufp = new float[dims.iProjV*dims.iProjU];
-	unsigned int pitch = data.pitch / sizeof(float);
-
-	for (int i = 0; i < dims.iProjAngles; ++i) {
-		for (int j = 0; j < dims.iProjV; ++j) {
-			cudaMemcpy(bufp+dims.iProjU*j, ((float*)data.ptr)+pitch*dims.iProjAngles*j+pitch*i, dims.iProjU*sizeof(float), cudaMemcpyDeviceToHost);
-		}
-
-		char fname[512];
-		sprintf(fname, filespec, i);
-		saveImage(fname, dims.iProjV, dims.iProjU, bufp, fMin, fMax);
-	}
-}
-
-
-
-
-int main()
-{
-#if 0
-	SDimensions3D dims;
-	dims.iVolX = 512;
-	dims.iVolY = 512;
-	dims.iVolZ = 512;
-	dims.iProjAngles = 180;
-	dims.iProjU = 1024;
-	dims.iProjV = 1024;
-	dims.iRaysPerDet = 1;
-
-	cudaExtent extentV;
-	extentV.width = dims.iVolX*sizeof(float);
-	extentV.height = dims.iVolY;
-	extentV.depth = dims.iVolZ;
-
-	cudaPitchedPtr volData; // pitch, ptr, xsize, ysize
-
-	cudaMalloc3D(&volData, extentV);
-
-	cudaExtent extentP;
-	extentP.width = dims.iProjU*sizeof(float);
-	extentP.height = dims.iProjAngles;
-	extentP.depth = dims.iProjV;
-
-	cudaPitchedPtr projData; // pitch, ptr, xsize, ysize
-
-	cudaMalloc3D(&projData, extentP);
-	cudaMemset3D(projData, 0, extentP);
-
-#if 0
-	float* slice = new float[256*256];
-	cudaPitchedPtr ptr;
-	ptr.ptr = slice;
-	ptr.pitch = 256*sizeof(float);
-	ptr.xsize = 256*sizeof(float);
-	ptr.ysize = 256;
-
-	for (unsigned int i = 0; i < 256*256; ++i)
-		slice[i] = 1.0f;
-	for (unsigned int i = 0; i < 256; ++i) {
-		cudaExtent extentS;
-		extentS.width = dims.iVolX*sizeof(float);
-		extentS.height = dims.iVolY;
-		extentS.depth = 1;
-		cudaPos sp = { 0, 0, 0 };
-		cudaPos dp = { 0, 0, i };
-		cudaMemcpy3DParms p;
-		p.srcArray = 0;
-		p.srcPos = sp;
-		p.srcPtr = ptr;
-		p.dstArray = 0;
-		p.dstPos = dp;
-		p.dstPtr = volData;
-		p.extent = extentS;
-		p.kind = cudaMemcpyHostToDevice;
-		cudaMemcpy3D(&p);
-#if 0
-		if (i == 128) {
-			for (unsigned int j = 0; j < 256*256; ++j)
-				slice[j] = 0.0f;
-		}
-#endif 
-	}
-#endif
-
-	SConeProjection angle[180];
-	angle[0].fSrcX = -1536;
-	angle[0].fSrcY = 0;
-	angle[0].fSrcZ = 0;
-
-	angle[0].fDetSX = 1024;
-	angle[0].fDetSY = -512;
-	angle[0].fDetSZ = 512;
-
-	angle[0].fDetUX = 0;
-	angle[0].fDetUY = 1;
-	angle[0].fDetUZ = 0;
-
-	angle[0].fDetVX = 0;
-	angle[0].fDetVY = 0;
-	angle[0].fDetVZ = -1;
-
-#define ROTATE0(name,i,alpha) do { angle[i].f##name##X = angle[0].f##name##X * cos(alpha) - angle[0].f##name##Y * sin(alpha); angle[i].f##name##Y = angle[0].f##name##X * sin(alpha) + angle[0].f##name##Y * cos(alpha); } while(0)
-	for (int i = 1; i < 180; ++i) {
-		angle[i] = angle[0];
-		ROTATE0(Src, i, i*2*M_PI/180);
-		ROTATE0(DetS, i, i*2*M_PI/180);
-		ROTATE0(DetU, i, i*2*M_PI/180);
-		ROTATE0(DetV, i, i*2*M_PI/180);
-	}
-#undef ROTATE0
-
-	astraCUDA3d::ConeFP(volData, projData, dims, angle, 1.0f);
-
-	//dumpSinograms("sino%03d.png", projData, dims, 0, 512);
-	//dumpProjections("proj%03d.png", projData, dims, 0, 512);
-
-	astraCUDA3d::zeroVolumeData(volData, dims);
-
-	float* angles = new float[dims.iProjAngles];
-	for (int i = 0; i < 180; ++i)
-		angles[i] = i*2*M_PI/180;
-
-	astraCUDA3d::FDK(volData, projData, 1536, 512, 0, 0, dims, angles);
-
-	dumpVolume("vol%03d.png", volData, dims, -20, 100);
-
-
-#else
-
-	SDimensions3D dims;
-	dims.iVolX = 1000;
-	dims.iVolY = 999;
-	dims.iVolZ = 500;
-	dims.iProjAngles = 376;
-	dims.iProjU = 1024;
-	dims.iProjV = 524;
-	dims.iRaysPerDet = 1;
-
-	float* angles = new float[dims.iProjAngles];
-	for (int i = 0; i < dims.iProjAngles; ++i)
-		angles[i] = -i*(M_PI)/360;
-
-	cudaPitchedPtr volData = astraCUDA3d::allocateVolumeData(dims);
-	cudaPitchedPtr projData = astraCUDA3d::allocateProjectionData(dims);
-	astraCUDA3d::zeroProjectionData(projData, dims);
-	astraCUDA3d::zeroVolumeData(volData, dims);
-
-	timeval t;
-	tic(t);
-
-	for (int i = 0; i < dims.iProjAngles; ++i) {
-		char fname[256];
-		sprintf(fname, "/home/wpalenst/tmp/Elke/proj%04d.png", i);
-		unsigned int w,h;
-		float* bufp = loadImage(fname, w,h);
-
-		int pitch = projData.pitch / sizeof(float);
-		for (int j = 0; j < dims.iProjV; ++j) {
-			cudaMemcpy(((float*)projData.ptr)+dims.iProjAngles*pitch*j+pitch*i, bufp+dims.iProjU*j, dims.iProjU*sizeof(float), cudaMemcpyHostToDevice);
-		}
-
-		delete[] bufp;
-	}
-	printf("Load time: %f\n", toc(t));
-
-	//dumpSinograms("sino%03d.png", projData, dims, -8.0f, 256.0f);
-	//astraCUDA3d::FDK(volData, projData, 7350, 62355, 0, 10, dims, angles);
-	//astraCUDA3d::FDK(volData, projData, 7350, -380, 0, 10, dims, angles);
-
-	tic(t);
-
-	astraCUDA3d::FDK(volData, projData, 7383.29867, 0, 0, 10, dims, angles);
-
-	printf("FDK time: %f\n", toc(t));
-	tic(t);
-
-	dumpVolume("vol%03d.png", volData, dims, -65.9f, 200.0f);
-	//dumpVolume("vol%03d.png", volData, dims, 0.0f, 256.0f);
-	printf("Save time: %f\n", toc(t));
-
-#endif
-
-
-}
-#endif
diff --git a/cuda/3d/mem3d.cu b/cuda/3d/mem3d.cu
index 697d2d2..50cfe75 100644
--- a/cuda/3d/mem3d.cu
+++ b/cuda/3d/mem3d.cu
@@ -268,7 +268,7 @@ bool FP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, con
 	return ok;
 }
 
-bool BP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, int iVoxelSuperSampling, bool bFDKWeighting)
+bool BP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, int iVoxelSuperSampling)
 {
 	assert(!volData.d->arr);
 	SDimensions3D dims;
@@ -289,7 +289,7 @@ bool BP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, con
 	                          pParProjs, pConeProjs,
 	                          params);
 
-	params.bFDKWeighting = bFDKWeighting;
+	params.bFDKWeighting = false;
 
 	if (pParProjs) {
 		if (projData.d->arr)
diff --git a/cuda/3d/par3d_bp.cu b/cuda/3d/par3d_bp.cu
index 3656f78..602f209 100644
--- a/cuda/3d/par3d_bp.cu
+++ b/cuda/3d/par3d_bp.cu
@@ -28,11 +28,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/3d/util3d.h"
 #include "astra/cuda/3d/dims3d.h"
 
-#ifdef STANDALONE
-#include "astra/cuda/3d/par3d_fp.h"
-#include "testutil.h"
-#endif
-
 #include <cstdio>
 #include <cassert>
 #include <iostream>
@@ -55,7 +50,13 @@ static const unsigned int g_volBlockY = 32;
 
 static const unsigned g_MaxAngles = 1024;
 
-__constant__ float gC_C[8*g_MaxAngles];
+struct DevPar3DParams {
+	float4 fNumU;
+	float4 fNumV;
+};
+
+__constant__ DevPar3DParams gC_C[g_MaxAngles];
+__constant__ float gC_scale[g_MaxAngles];
 
 
 static bool bindProjDataTexture(const cudaArray* array)
@@ -115,8 +116,9 @@ __global__ void dev_par3D_BP(void* D_volData, unsigned int volPitch, int startAn
 		for (int angle = startAngle; angle < endAngle; ++angle, fAngle += 1.0f)
 		{
 
-			float4 fCu = make_float4(gC_C[8*angle+0], gC_C[8*angle+1], gC_C[8*angle+2], gC_C[8*angle+3]);
-			float4 fCv = make_float4(gC_C[8*angle+4], gC_C[8*angle+5], gC_C[8*angle+6], gC_C[8*angle+7]);
+			float4 fCu = gC_C[angle].fNumU;
+			float4 fCv = gC_C[angle].fNumV;
+			float fS = gC_scale[angle];
 
 			float fU = fCu.w + fX * fCu.x + fY * fCu.y + fZ * fCu.z;
 			float fV = fCv.w + fX * fCv.x + fY * fCv.y + fZ * fCv.z;
@@ -124,7 +126,7 @@ __global__ void dev_par3D_BP(void* D_volData, unsigned int volPitch, int startAn
 			for (int idx = 0; idx < ZSIZE; ++idx) {
 
 				float fVal = tex3D(gT_par3DProjTexture, fU, fAngle, fV);
-				Z[idx] += fVal;
+				Z[idx] += fVal * fS;
 
 				fU += fCu.z;
 				fV += fCv.z;
@@ -190,14 +192,9 @@ __global__ void dev_par3D_BP_SS(void* D_volData, unsigned int volPitch, int star
 
 		for (int angle = startAngle; angle < endAngle; ++angle, fAngle += 1.0f)
 		{
-			const float fCux = gC_C[8*angle+0];
-			const float fCuy = gC_C[8*angle+1];
-			const float fCuz = gC_C[8*angle+2];
-			const float fCuc = gC_C[8*angle+3];
-			const float fCvx = gC_C[8*angle+4];
-			const float fCvy = gC_C[8*angle+5];
-			const float fCvz = gC_C[8*angle+6];
-			const float fCvc = gC_C[8*angle+7];
+			float4 fCu = gC_C[angle].fNumU;
+			float4 fCv = gC_C[angle].fNumV;
+			float fS = gC_scale[angle];
 
 			float fXs = fX;
 			for (int iSubX = 0; iSubX < iRaysPerVoxelDim; ++iSubX) {
@@ -206,10 +203,10 @@ __global__ void dev_par3D_BP_SS(void* D_volData, unsigned int volPitch, int star
 			float fZs = fZ;
 			for (int iSubZ = 0; iSubZ < iRaysPerVoxelDim; ++iSubZ) {
 
-				const float fU = fCuc + fXs * fCux + fYs * fCuy + fZs * fCuz;
-				const float fV = fCvc + fXs * fCvx + fYs * fCvy + fZs * fCvz;
+				const float fU = fCu.w + fXs * fCu.x + fYs * fCu.y + fZs * fCu.z;
+				const float fV = fCv.w + fXs * fCv.x + fYs * fCv.y + fZs * fCv.z;
 
-				fVal += tex3D(gT_par3DProjTexture, fU, fAngle, fV);
+				fVal += tex3D(gT_par3DProjTexture, fU, fAngle, fV) * fS;
 				fZs += fSubStep;
 			}
 			fYs += fSubStep;
@@ -224,6 +221,35 @@ __global__ void dev_par3D_BP_SS(void* D_volData, unsigned int volPitch, int star
 
 }
 
+bool transferConstants(const SPar3DProjection* angles, unsigned int iProjAngles, const SProjectorParams3D& params)
+{
+	DevPar3DParams *p = new DevPar3DParams[iProjAngles];
+	float *s = new float[iProjAngles];
+
+	for (unsigned int i = 0; i < iProjAngles; ++i) {
+		Vec3 u(angles[i].fDetUX, angles[i].fDetUY, angles[i].fDetUZ);
+		Vec3 v(angles[i].fDetVX, angles[i].fDetVY, angles[i].fDetVZ);
+		Vec3 r(angles[i].fRayX, angles[i].fRayY, angles[i].fRayZ);
+		Vec3 d(angles[i].fDetSX, angles[i].fDetSY, angles[i].fDetSZ);
+
+		double fDen = det3(r,u,v);
+		p[i].fNumU.x = -det3x(r,v) / fDen;
+		p[i].fNumU.y = -det3y(r,v) / fDen;
+		p[i].fNumU.z = -det3z(r,v) / fDen;
+		p[i].fNumU.w = -det3(r,d,v) / fDen;
+		p[i].fNumV.x = det3x(r,u) / fDen;
+		p[i].fNumV.y = det3y(r,u) / fDen;
+		p[i].fNumV.z = det3z(r,u) / fDen;
+		p[i].fNumV.w = det3(r,d,u) / fDen;
+
+		s[i] = 1.0 / scaled_cross3(u,v,Vec3(params.fVolScaleX,params.fVolScaleY,params.fVolScaleZ)).norm();
+	}
+
+	cudaMemcpyToSymbol(gC_C, p, iProjAngles*sizeof(DevPar3DParams), 0, cudaMemcpyHostToDevice);
+	cudaMemcpyToSymbol(gC_scale, s, iProjAngles*sizeof(float), 0, cudaMemcpyHostToDevice);
+	return true;
+}
+
 bool Par3DBP_Array(cudaPitchedPtr D_volumeData,
                    cudaArray *D_projArray,
                    const SDimensions3D& dims, const SPar3DProjection* angles,
@@ -238,33 +264,9 @@ bool Par3DBP_Array(cudaPitchedPtr D_volumeData,
 		if (th + angleCount > dims.iProjAngles)
 			angleCount = dims.iProjAngles - th;
 
-		// transfer angles to constant memory
-		float* tmp = new float[8*dims.iProjAngles];
-
-		// NB: We increment angles at the end of the loop body.
-
-
-		// TODO: Use functions from dims3d.cu for this:
-
-#define TRANSFER_TO_CONSTANT(expr,name) do { for (unsigned int i = 0; i < angleCount; ++i) tmp[8*i + name] = (expr) ; } while (0)
-
-#define DENOM (angles[i].fRayX*angles[i].fDetUY*angles[i].fDetVZ - angles[i].fRayX*angles[i].fDetUZ*angles[i].fDetVY - angles[i].fRayY*angles[i].fDetUX*angles[i].fDetVZ + angles[i].fRayY*angles[i].fDetUZ*angles[i].fDetVX + angles[i].fRayZ*angles[i].fDetUX*angles[i].fDetVY - angles[i].fRayZ*angles[i].fDetUY*angles[i].fDetVX)
-
-		TRANSFER_TO_CONSTANT( ( - (angles[i].fRayY*angles[i].fDetVZ - angles[i].fRayZ*angles[i].fDetVY)) / DENOM , 0 );
-		TRANSFER_TO_CONSTANT( ( (angles[i].fRayX*angles[i].fDetVZ - angles[i].fRayZ*angles[i].fDetVX)) / DENOM , 1 );
-		TRANSFER_TO_CONSTANT( (- (angles[i].fRayX*angles[i].fDetVY - angles[i].fRayY*angles[i].fDetVX) ) / DENOM , 2 );
-		TRANSFER_TO_CONSTANT( (-(angles[i].fDetSY*angles[i].fDetVZ - angles[i].fDetSZ*angles[i].fDetVY)*angles[i].fRayX + (angles[i].fRayY*angles[i].fDetVZ - angles[i].fRayZ*angles[i].fDetVY)*angles[i].fDetSX - (angles[i].fRayY*angles[i].fDetSZ - angles[i].fRayZ*angles[i].fDetSY)*angles[i].fDetVX) / DENOM , 3 );
-
-		TRANSFER_TO_CONSTANT( ((angles[i].fRayY*angles[i].fDetUZ - angles[i].fRayZ*angles[i].fDetUY) ) / DENOM , 4 );
-		TRANSFER_TO_CONSTANT( (- (angles[i].fRayX*angles[i].fDetUZ - angles[i].fRayZ*angles[i].fDetUX) ) / DENOM , 5 );
-		TRANSFER_TO_CONSTANT( ((angles[i].fRayX*angles[i].fDetUY - angles[i].fRayY*angles[i].fDetUX) ) / DENOM , 6 );
-		TRANSFER_TO_CONSTANT( ((angles[i].fDetSY*angles[i].fDetUZ - angles[i].fDetSZ*angles[i].fDetUY)*angles[i].fRayX - (angles[i].fRayY*angles[i].fDetUZ - angles[i].fRayZ*angles[i].fDetUY)*angles[i].fDetSX + (angles[i].fRayY*angles[i].fDetSZ - angles[i].fRayZ*angles[i].fDetSY)*angles[i].fDetUX ) / DENOM , 7 );
-
-#undef TRANSFER_TO_CONSTANT
-#undef DENOM
-		cudaMemcpyToSymbol(gC_C, tmp, angleCount*8*sizeof(float), 0, cudaMemcpyHostToDevice); 
-
-		delete[] tmp;
+		bool ok = transferConstants(angles, dims.iProjAngles, params);
+		if (!ok)
+			return false;
 
 		dim3 dimBlock(g_volBlockX, g_volBlockY);
 
@@ -310,161 +312,3 @@ bool Par3DBP(cudaPitchedPtr D_volumeData,
 
 
 }
-
-#ifdef STANDALONE
-int main()
-{
-	SDimensions3D dims;
-	dims.iVolX = 256;
-	dims.iVolY = 256;
-	dims.iVolZ = 256;
-	dims.iProjAngles = 180;
-	dims.iProjU = 512;
-	dims.iProjV = 512;
-	dims.iRaysPerDet = 1;
-
-	cudaExtent extentV;
-	extentV.width = dims.iVolX*sizeof(float);
-	extentV.height = dims.iVolY;
-	extentV.depth = dims.iVolZ;
-
-	cudaPitchedPtr volData; // pitch, ptr, xsize, ysize
-
-	cudaMalloc3D(&volData, extentV);
-
-	cudaExtent extentP;
-	extentP.width = dims.iProjU*sizeof(float);
-	extentP.height = dims.iProjAngles;
-	extentP.depth = dims.iProjV;
-
-	cudaPitchedPtr projData; // pitch, ptr, xsize, ysize
-
-	cudaMalloc3D(&projData, extentP);
-	cudaMemset3D(projData, 0, extentP);
-
-	float* slice = new float[256*256];
-	cudaPitchedPtr ptr;
-	ptr.ptr = slice;
-	ptr.pitch = 256*sizeof(float);
-	ptr.xsize = 256*sizeof(float);
-	ptr.ysize = 256;
-
-	for (unsigned int i = 0; i < 256*256; ++i)
-		slice[i] = 1.0f;
-	for (unsigned int i = 0; i < 256; ++i) {
-		cudaExtent extentS;
-		extentS.width = dims.iVolX*sizeof(float);
-		extentS.height = dims.iVolY;
-		extentS.depth = 1;
-		cudaPos sp = { 0, 0, 0 };
-		cudaPos dp = { 0, 0, i };
-		cudaMemcpy3DParms p;
-		p.srcArray = 0;
-		p.srcPos = sp;
-		p.srcPtr = ptr;
-		p.dstArray = 0;
-		p.dstPos = dp;
-		p.dstPtr = volData;
-		p.extent = extentS;
-		p.kind = cudaMemcpyHostToDevice;
-		cudaMemcpy3D(&p);
-#if 0
-		if (i == 128) {
-			for (unsigned int j = 0; j < 256*256; ++j)
-				slice[j] = 0.0f;
-		}
-#endif 
-	}
-
-
-	SPar3DProjection angle[180];
-	angle[0].fRayX = 1;
-	angle[0].fRayY = 0;
-	angle[0].fRayZ = 0;
-
-	angle[0].fDetSX = 512;
-	angle[0].fDetSY = -256;
-	angle[0].fDetSZ = -256;
-
-	angle[0].fDetUX = 0;
-	angle[0].fDetUY = 1;
-	angle[0].fDetUZ = 0;
-
-	angle[0].fDetVX = 0;
-	angle[0].fDetVY = 0;
-	angle[0].fDetVZ = 1;
-
-#define ROTATE0(name,i,alpha) do { angle[i].f##name##X = angle[0].f##name##X * cos(alpha) - angle[0].f##name##Y * sin(alpha); angle[i].f##name##Y = angle[0].f##name##X * sin(alpha) + angle[0].f##name##Y * cos(alpha); } while(0)
-	for (int i = 1; i < 180; ++i) {
-		angle[i] = angle[0];
-		ROTATE0(Ray, i, i*2*M_PI/180);
-		ROTATE0(DetS, i, i*2*M_PI/180);
-		ROTATE0(DetU, i, i*2*M_PI/180);
-		ROTATE0(DetV, i, i*2*M_PI/180);
-	}
-#undef ROTATE0
-
-	astraCUDA3d::Par3DFP(volData, projData, dims, angle, 1.0f);
-#if 1
-	float* bufs = new float[180*512];
-
-	for (int i = 0; i < 512; ++i) {
-		cudaMemcpy(bufs, ((float*)projData.ptr)+180*512*i, 180*512*sizeof(float), cudaMemcpyDeviceToHost);
-
-		printf("%d %d %d\n", projData.pitch, projData.xsize, projData.ysize);
-
-		char fname[20];
-		sprintf(fname, "sino%03d.png", i);
-		saveImage(fname, 180, 512, bufs, 0, 512);
-	}
-
-	float* bufp = new float[512*512];
-
-	for (int i = 0; i < 180; ++i) {
-		for (int j = 0; j < 512; ++j) {
-			cudaMemcpy(bufp+512*j, ((float*)projData.ptr)+180*512*j+512*i, 512*sizeof(float), cudaMemcpyDeviceToHost);
-		}
-
-		char fname[20];
-		sprintf(fname, "proj%03d.png", i);
-		saveImage(fname, 512, 512, bufp, 0, 512);
-	}
-#endif		
-	for (unsigned int i = 0; i < 256*256; ++i)
-		slice[i] = 0.0f;
-	for (unsigned int i = 0; i < 256; ++i) {
-		cudaExtent extentS;
-		extentS.width = dims.iVolX*sizeof(float);
-		extentS.height = dims.iVolY;
-		extentS.depth = 1;
-		cudaPos sp = { 0, 0, 0 };
-		cudaPos dp = { 0, 0, i };
-		cudaMemcpy3DParms p;
-		p.srcArray = 0;
-		p.srcPos = sp;
-		p.srcPtr = ptr;
-		p.dstArray = 0;
-		p.dstPos = dp;
-		p.dstPtr = volData;
-		p.extent = extentS;
-		p.kind = cudaMemcpyHostToDevice;
-		cudaMemcpy3D(&p);
-	}
-
-	astraCUDA3d::Par3DBP(volData, projData, dims, angle, 1.0f);
-#if 1
-	float* buf = new float[256*256];
-
-	for (int i = 0; i < 256; ++i) {
-		cudaMemcpy(buf, ((float*)volData.ptr)+256*256*i, 256*256*sizeof(float), cudaMemcpyDeviceToHost);
-
-		printf("%d %d %d\n", volData.pitch, volData.xsize, volData.ysize);
-
-		char fname[20];
-		sprintf(fname, "vol%03d.png", i);
-		saveImage(fname, 256, 256, buf, 0, 60000);
-	}
-#endif
-
-}
-#endif
diff --git a/cuda/3d/par3d_fp.cu b/cuda/3d/par3d_fp.cu
index 515b1ba..0a4a5cc 100644
--- a/cuda/3d/par3d_fp.cu
+++ b/cuda/3d/par3d_fp.cu
@@ -28,11 +28,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/3d/util3d.h"
 #include "astra/cuda/3d/dims3d.h"
 
-#ifdef STANDALONE
-#include "testutil.h"
-#endif
-
-
 #include <cstdio>
 #include <cassert>
 #include <iostream>
@@ -751,166 +746,3 @@ bool Par3DFP_SumSqW(cudaPitchedPtr D_volumeData,
 
 
 }
-
-#ifdef STANDALONE
-
-using namespace astraCUDA3d;
-
-int main()
-{
-	cudaSetDevice(1);
-
-
-	SDimensions3D dims;
-	dims.iVolX = 500;
-	dims.iVolY = 500;
-	dims.iVolZ = 81;
-	dims.iProjAngles = 241;
-	dims.iProjU = 600;
-	dims.iProjV = 100;
-	dims.iRaysPerDet = 1;
-
-	SPar3DProjection base;
-	base.fRayX = 1.0f;
-	base.fRayY = 0.0f;
-	base.fRayZ = 0.1f;
-
-	base.fDetSX = 0.0f;
-	base.fDetSY = -300.0f;
-	base.fDetSZ = -50.0f;
-
-	base.fDetUX = 0.0f;
-	base.fDetUY = 1.0f;
-	base.fDetUZ = 0.0f;
-
-	base.fDetVX = 0.0f;
-	base.fDetVY = 0.0f;
-	base.fDetVZ = 1.0f;
-
-	SPar3DProjection angle[dims.iProjAngles];
-
-	cudaPitchedPtr volData; // pitch, ptr, xsize, ysize
-
-	volData = allocateVolumeData(dims);
-
-	cudaPitchedPtr projData; // pitch, ptr, xsize, ysize
-
-	projData = allocateProjectionData(dims);
-
-	unsigned int ix = 500,iy = 500;
-
-	float* buf = new float[dims.iProjU*dims.iProjV];
-
-	float* slice = new float[dims.iVolX*dims.iVolY];
-	for (int i = 0; i < dims.iVolX*dims.iVolY; ++i)
-		slice[i] = 1.0f;
-
-	for (unsigned int a = 0; a < 241; a += dims.iProjAngles) {
-
-		zeroProjectionData(projData, dims);
-
-		for (int y = 0; y < iy; y += dims.iVolY) {
-			for (int x = 0; x < ix; x += dims.iVolX) { 
-
-				timeval st;
-				tic(st);
-
-				for (int z = 0; z < dims.iVolZ; ++z) {
-//					char sfn[256];
-//					sprintf(sfn, "/home/wpalenst/projects/cone_simulation/phantom_4096/mouse_fem_phantom_%04d.png", 30+z);
-//					float* slice = loadSubImage(sfn, x, y, dims.iVolX, dims.iVolY);
-
-					cudaPitchedPtr ptr;
-					ptr.ptr = slice;
-					ptr.pitch = dims.iVolX*sizeof(float);
-					ptr.xsize = dims.iVolX*sizeof(float);
-					ptr.ysize = dims.iVolY;
-					cudaExtent extentS;
-					extentS.width = dims.iVolX*sizeof(float);
-					extentS.height = dims.iVolY;
-					extentS.depth = 1;
-
-					cudaPos sp = { 0, 0, 0 };
-					cudaPos dp = { 0, 0, z };
-					cudaMemcpy3DParms p;
-					p.srcArray = 0;
-					p.srcPos = sp;
-					p.srcPtr = ptr;
-					p.dstArray = 0;
-					p.dstPos = dp;
-					p.dstPtr = volData;
-					p.extent = extentS;
-					p.kind = cudaMemcpyHostToDevice;
-					cudaError err = cudaMemcpy3D(&p);
-					assert(!err);
-//					delete[] slice;
-				}
-
-				printf("Load: %f\n", toc(st));
-
-#if 0
-
-	cudaPos zp = { 0, 0, 0 };
-
-	cudaPitchedPtr t;
-	t.ptr = new float[1024*1024];
-	t.pitch = 1024*4;
-	t.xsize = 1024*4;
-	t.ysize = 1024;
-
-	cudaMemcpy3DParms p;
-	p.srcArray = 0;
-	p.srcPos = zp;
-	p.srcPtr = volData;
-	p.extent = extentS;
-	p.dstArray = 0;
-	p.dstPtr = t;
-	p.dstPos = zp;
-	p.kind = cudaMemcpyDeviceToHost;
-	cudaError err = cudaMemcpy3D(&p);
-	assert(!err);
-
-	char fn[32];
-	sprintf(fn, "t%d%d.png", x / dims.iVolX, y / dims.iVolY);
-	saveImage(fn, 1024, 1024, (float*)t.ptr);
-	saveImage("s.png", 4096, 4096, slice);
-	delete[] (float*)t.ptr;
-#endif
-
-
-#define ROTATE0(name,i,alpha) do { angle[i].f##name##X = base.f##name##X * cos(alpha) - base.f##name##Y * sin(alpha); angle[i].f##name##Y = base.f##name##X * sin(alpha) + base.f##name##Y * cos(alpha); angle[i].f##name##Z = base.f##name##Z; } while(0)
-#define SHIFT(name,i,x,y) do { angle[i].f##name##X += x; angle[i].f##name##Y += y; } while(0)
-				for (int i = 0; i < dims.iProjAngles; ++i) {
-					ROTATE0(Ray, i, (a+i)*.8*M_PI/180);
-					ROTATE0(DetS, i, (a+i)*.8*M_PI/180);
-					ROTATE0(DetU, i, (a+i)*.8*M_PI/180);
-					ROTATE0(DetV, i, (a+i)*.8*M_PI/180);
-
-
-//					SHIFT(Src, i, (-x+1536), (-y+1536));
-//					SHIFT(DetS, i, (-x+1536), (-y+1536));
-				}
-#undef ROTATE0
-#undef SHIFT
-				tic(st);
-
-				astraCUDA3d::Par3DFP(volData, projData, dims, angle, 1.0f);
-
-				printf("FP: %f\n", toc(st));
-
-			}
-		}
-		for (unsigned int aa = 0; aa < dims.iProjAngles; ++aa) {
-			for (unsigned int v = 0; v < dims.iProjV; ++v)
-				cudaMemcpy(buf+v*dims.iProjU, ((float*)projData.ptr)+(v*dims.iProjAngles+aa)*(projData.pitch/sizeof(float)), dims.iProjU*sizeof(float), cudaMemcpyDeviceToHost);
-
-			char fname[32];
-			sprintf(fname, "proj%03d.png", a+aa);
-			saveImage(fname, dims.iProjV, dims.iProjU, buf, 0.0f, 1000.0f);
-		}
-	}
-
-	delete[] buf;
-
-}
-#endif
diff --git a/cuda/3d/sirt3d.cu b/cuda/3d/sirt3d.cu
index 869b2fd..e68bde8 100644
--- a/cuda/3d/sirt3d.cu
+++ b/cuda/3d/sirt3d.cu
@@ -30,10 +30,6 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/cuda/3d/arith3d.h"
 #include "astra/cuda/3d/cone_fp.h"
 
-#ifdef STANDALONE
-#include "testutil.h"
-#endif
-
 #include <cstdio>
 #include <cassert>
 
@@ -375,160 +371,3 @@ bool doSIRT(cudaPitchedPtr& D_volumeData,
 
 }
 
-#ifdef STANDALONE
-
-using namespace astraCUDA3d;
-
-int main()
-{
-	SDimensions3D dims;
-	dims.iVolX = 256;
-	dims.iVolY = 256;
-	dims.iVolZ = 256;
-	dims.iProjAngles = 100;
-	dims.iProjU = 512;
-	dims.iProjV = 512;
-	dims.iRaysPerDet = 1;
-
-	SConeProjection angle[100];
-	angle[0].fSrcX = -2905.6;
-	angle[0].fSrcY = 0;
-	angle[0].fSrcZ = 0;
-
-	angle[0].fDetSX = 694.4;
-	angle[0].fDetSY = -122.4704;
-	angle[0].fDetSZ = -122.4704;
-
-	angle[0].fDetUX = 0;
-	angle[0].fDetUY = .4784;
-	//angle[0].fDetUY = .5;
-	angle[0].fDetUZ = 0;
-
-	angle[0].fDetVX = 0;
-	angle[0].fDetVY = 0;
-	angle[0].fDetVZ = .4784;
-
-#define ROTATE0(name,i,alpha) do { angle[i].f##name##X = angle[0].f##name##X * cos(alpha) - angle[0].f##name##Y * sin(alpha); angle[i].f##name##Y = angle[0].f##name##X * sin(alpha) + angle[0].f##name##Y * cos(alpha); } while(0)
-	for (int i = 1; i < 100; ++i) {
-		angle[i] = angle[0];
-		ROTATE0(Src, i, i*2*M_PI/100);
-		ROTATE0(DetS, i, i*2*M_PI/100);
-		ROTATE0(DetU, i, i*2*M_PI/100);
-		ROTATE0(DetV, i, i*2*M_PI/100);
-	}
-#undef ROTATE0
-
-
-	cudaPitchedPtr volData = allocateVolumeData(dims);
-	cudaPitchedPtr projData = allocateProjectionData(dims);
-	zeroProjectionData(projData, dims);
-
-	float* pbuf = new float[100*512*512];
-	copyProjectionsFromDevice(pbuf, projData, dims);
-	copyProjectionsToDevice(pbuf, projData, dims);
-	delete[] pbuf;
-
-#if 0
-	float* slice = new float[256*256];
-	cudaPitchedPtr ptr;
-	ptr.ptr = slice;
-	ptr.pitch = 256*sizeof(float);
-	ptr.xsize = 256*sizeof(float);
-	ptr.ysize = 256;
-
-	for (unsigned int i = 0; i < 256; ++i) {
-		for (unsigned int y = 0; y < 256; ++y)
-			for (unsigned int x = 0; x < 256; ++x)
-				slice[y*256+x] = (i-127.5)*(i-127.5)+(y-127.5)*(y-127.5)+(x-127.5)*(x-127.5) < 4900 ? 1.0f : 0.0f;
-
-		cudaExtent extentS;
-		extentS.width = dims.iVolX*sizeof(float);
-		extentS.height = dims.iVolY;
-		extentS.depth = 1;
-		cudaPos sp = { 0, 0, 0 };
-		cudaPos dp = { 0, 0, i };
-		cudaMemcpy3DParms p;
-		p.srcArray = 0;
-		p.srcPos = sp;
-		p.srcPtr = ptr;
-		p.dstArray = 0;
-		p.dstPos = dp;
-		p.dstPtr = volData;
-		p.extent = extentS;
-		p.kind = cudaMemcpyHostToDevice;
-		cudaMemcpy3D(&p);
-	}
-	astraCUDA3d::ConeFP(volData, projData, dims, angle, 1.0f);
-
-#else
-
-	for (int i = 0; i < 100; ++i) {
-		char fname[32];
-		sprintf(fname, "Tiffs/%04d.png", 4*i);
-		unsigned int w,h;
-		float* bufp = loadImage(fname, w,h);
-
-		for (int j = 0; j < 512*512; ++j) {
-			float v = bufp[j];
-			if (v > 236.0f) v = 236.0f;
-			v = logf(236.0f / v);
-			bufp[j] = 256*v;
-		}
-
-		for (int j = 0; j < 512; ++j) {
-			cudaMemcpy(((float*)projData.ptr)+100*512*j+512*i, bufp+512*j, 512*sizeof(float), cudaMemcpyHostToDevice);
-		}
-
-		delete[] bufp;
-
-	}
-#endif
-
-#if 0
-	float* bufs = new float[100*512];
-
-	for (int i = 0; i < 512; ++i) {
-		cudaMemcpy(bufs, ((float*)projData.ptr)+100*512*i, 100*512*sizeof(float), cudaMemcpyDeviceToHost);
-
-		printf("%d %d %d\n", projData.pitch, projData.xsize, projData.ysize);
-
-		char fname[20];
-		sprintf(fname, "sino%03d.png", i);
-		saveImage(fname, 100, 512, bufs);
-	}
-
-	float* bufp = new float[512*512];
-
-	for (int i = 0; i < 100; ++i) {
-		for (int j = 0; j < 512; ++j) {
-			cudaMemcpy(bufp+512*j, ((float*)projData.ptr)+100*512*j+512*i, 512*sizeof(float), cudaMemcpyDeviceToHost);
-		}
-
-		char fname[20];
-		sprintf(fname, "proj%03d.png", i);
-		saveImage(fname, 512, 512, bufp);
-	}
-#endif
-
-	zeroVolumeData(volData, dims);
-
-	cudaPitchedPtr maskData;
-	maskData.ptr = 0;
-
-	astraCUDA3d::doSIRT(volData, projData, maskData, dims, angle, 50);
-#if 1
-	float* buf = new float[256*256];
-
-	for (int i = 0; i < 256; ++i) {
-		cudaMemcpy(buf, ((float*)volData.ptr)+256*256*i, 256*256*sizeof(float), cudaMemcpyDeviceToHost);
-
-		char fname[20];
-		sprintf(fname, "vol%03d.png", i);
-		saveImage(fname, 256, 256, buf);
-	}
-#endif
-
-	return 0;
-}
-#endif
-
diff --git a/include/astra/CudaProjector3D.h b/include/astra/CudaProjector3D.h
index 60df7bb..9b4ff1f 100644
--- a/include/astra/CudaProjector3D.h
+++ b/include/astra/CudaProjector3D.h
@@ -117,7 +117,6 @@ public:
 	int getVoxelSuperSampling() const { return m_iVoxelSuperSampling; }
 	int getDetectorSuperSampling() const { return m_iDetectorSuperSampling; }
 	int getGPUIndex() const { return m_iGPUIndex; }
-	bool getDensityWeighting() const { return m_bDensityWeighting; }
 
 protected:
 
@@ -125,7 +124,6 @@ protected:
 	int m_iVoxelSuperSampling;
 	int m_iDetectorSuperSampling;
 	int m_iGPUIndex;
-	bool m_bDensityWeighting;
 
 };
 
diff --git a/include/astra/FanFlatBeamLineKernelProjector2D.inl b/include/astra/FanFlatBeamLineKernelProjector2D.inl
index 58dec61..8f2e673 100644
--- a/include/astra/FanFlatBeamLineKernelProjector2D.inl
+++ b/include/astra/FanFlatBeamLineKernelProjector2D.inl
@@ -82,8 +82,6 @@ void CFanFlatBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, in
 
 		const SFanProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
 
-		float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY);
-
 		// loop detectors
 		for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
 			
@@ -104,7 +102,7 @@ void CFanFlatBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, in
 			// vertically
 			if (vertical) {
 				RxOverRy = Rx/Ry;
-				lengthPerRow = detSize * pixelLengthX * sqrt(Rx*Rx + Ry*Ry) / abs(Ry);
+				lengthPerRow = pixelLengthX * sqrt(Rx*Rx + Ry*Ry) / abs(Ry);
 				deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX;
 				S = 0.5f - 0.5f*fabs(RxOverRy);
 				T = 0.5f + 0.5f*fabs(RxOverRy);
@@ -154,7 +152,7 @@ void CFanFlatBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, in
 			// horizontally
 			else {
 				RyOverRx = Ry/Rx;
-				lengthPerCol = detSize * pixelLengthY * sqrt(Rx*Rx + Ry*Ry) / abs(Rx);
+				lengthPerCol = pixelLengthY * sqrt(Rx*Rx + Ry*Ry) / abs(Rx);
 				deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
 				S = 0.5f - 0.5f*fabs(RyOverRx);
 				T = 0.5f + 0.5f*fabs(RyOverRx);
diff --git a/include/astra/FanFlatBeamStripKernelProjector2D.inl b/include/astra/FanFlatBeamStripKernelProjector2D.inl
index 889d0a3..f5a688c 100644
--- a/include/astra/FanFlatBeamStripKernelProjector2D.inl
+++ b/include/astra/FanFlatBeamStripKernelProjector2D.inl
@@ -109,8 +109,12 @@ void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, i
 				// POLICY: RAY PRIOR
 				if (!p.rayPrior(iRayIndex)) continue;
 
-				float32 dist_srcDetPixSquared = projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() +
-				                                (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW * DW;
+				float32 dist_srcDetPixSquared = projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() + (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW * DW;
+				dist_srcDetPixSquared = dist_srcDetPixSquared * dist_srcDetPixSquared / (projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance()  * DW * DW);
+
+
+
+				//float32 InvRayWidthSquared = (projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance()) / dist_srcDetPixSquared;
 
 				float32 sin_theta_left, cos_theta_left;
 				float32 sin_theta_right, cos_theta_right;
@@ -257,8 +261,8 @@ void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, i
 				// POLICY: RAY PRIOR
 				if (!p.rayPrior(iRayIndex)) continue;
 
-				float32 dist_srcDetPixSquared = projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() +
-				                                (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW * DW;
+				float32 dist_srcDetPixSquared = projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance() + (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * (iDetector + 0.5f - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW * DW;
+				dist_srcDetPixSquared = dist_srcDetPixSquared * dist_srcDetPixSquared / (projgeom->getSourceDetectorDistance() * projgeom->getSourceDetectorDistance()  * DW * DW);
 
 				// get theta_l = alpha_left + theta and theta_r = alpha_right + theta
 				float32 sin_theta_left, cos_theta_left;
diff --git a/include/astra/Features.h b/include/astra/Features.h
index d88ae71..c7ef98c 100644
--- a/include/astra/Features.h
+++ b/include/astra/Features.h
@@ -38,10 +38,22 @@ _AstraExport bool hasFeature(const std::string &feature);
 
 FEATURES:
 
-cuda: is cuda support compiled in?
+cuda
+	is cuda support compiled in?
 	NB: To check if there is also actually a usable GPU, use cudaAvailable()
 
-mex_link: is there support for the matlab command astra_mex_data3d('link')?
+mex_link
+	is there support for the matlab command astra_mex_data3d('link')?
+
+projectors_scaled_as_line_integrals
+	This is set since all 2D and 3D, CPU and GPU projectors scale their outputs
+	to approximate line integrals. (Previously, some 2D projectors were scaled
+	as area integrals.)
+
+fan_cone_BP_density_weighting_by_default
+	This is set since fan beam and cone beam BP operations perform ray density
+	weighting by default to more closely approximate the true mathematical adjoint.
+	The DensityWeighting cuda3d projector option is removed.
 
 For future backward-incompatible changes, extra features will be added here
 
diff --git a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl
index aedcee9..3a18c6f 100644
--- a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl
+++ b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl
@@ -72,7 +72,7 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro
 	const int rowCount = m_pVolumeGeometry->getGridRowCount();
 
 	// Performance note:
-	// This is not a very well optimizated version of the distance driven
+	// This is not a very well optimized version of the distance driven
 	// projector. The CPU projector model in ASTRA requires ray-driven iteration,
 	// which limits re-use of intermediate computations.
 
@@ -86,6 +86,9 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro
 		const float32 Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f;
 		const float32 Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
 
+		const float32 rayWidth = fabs(proj->fDetUX * proj->fRayY - proj->fDetUY * proj->fRayX) /
+		                         sqrt(proj->fRayX * proj->fRayX + proj->fRayY * proj->fRayY);
+
 		// loop detectors
 		for (int iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
 
@@ -100,7 +103,7 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro
 			if (vertical) {
 
 				const float32 RxOverRy = proj->fRayX/proj->fRayY;
-				const float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY();
+				const float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY() / rayWidth;
 				const float32 deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX;
 				const float32 deltad = 0.5f * fabs((proj->fDetUX - proj->fDetUY * RxOverRy) * inv_pixelLengthX);
 
@@ -157,7 +160,7 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro
 			} else {
 
 				const float32 RyOverRx = proj->fRayY/proj->fRayX;
-				const float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY();
+				const float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY() / rayWidth;
 				const float32 deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
 				const float32 deltad = 0.5f * fabs((proj->fDetUY - proj->fDetUX * RyOverRx) * inv_pixelLengthY);
 
diff --git a/include/astra/ParallelBeamLineKernelProjector2D.inl b/include/astra/ParallelBeamLineKernelProjector2D.inl
index a9f1aa0..903ebb6 100644
--- a/include/astra/ParallelBeamLineKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLineKernelProjector2D.inl
@@ -166,24 +166,7 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i
 
 		const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
 
-		float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY);
-
 		bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY);
-		if (vertical) {
-			RxOverRy = proj->fRayX/proj->fRayY;
-			lengthPerRow = detSize * pixelLengthX * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY);
-			deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX;
-			S = 0.5f - 0.5f*fabs(RxOverRy);
-			T = 0.5f + 0.5f*fabs(RxOverRy);
-			invTminSTimesLengthPerRow = lengthPerRow / (T - S);
-		} else {
-			RyOverRx = proj->fRayY/proj->fRayX;
-			lengthPerCol = detSize * pixelLengthY * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX);
-			deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
-			S = 0.5f - 0.5f*fabs(RyOverRx);
-			T = 0.5f + 0.5f*fabs(RyOverRx);
-			invTminSTimesLengthPerCol = lengthPerCol / (T - S);
-		}
 
 		Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f;
 		Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
@@ -204,6 +187,13 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i
 			// vertically
 			if (vertical) {
 
+				RxOverRy = proj->fRayX/proj->fRayY;
+				lengthPerRow = pixelLengthX * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY);
+				deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX;
+				S = 0.5f - 0.5f*fabs(RxOverRy);
+				T = 0.5f + 0.5f*fabs(RxOverRy);
+				invTminSTimesLengthPerRow = lengthPerRow / (T - S);
+
 				// calculate c for row 0
 				c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX;
 
@@ -248,6 +238,13 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i
 			// horizontally
 			else {
 
+				RyOverRx = proj->fRayY/proj->fRayX;
+				lengthPerCol = pixelLengthY * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX);
+				deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
+				S = 0.5f - 0.5f*fabs(RyOverRx);
+				T = 0.5f + 0.5f*fabs(RyOverRx);
+				invTminSTimesLengthPerCol = lengthPerCol / (T - S);
+
 				// calculate r for col 0
 				r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY;
 
diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.inl b/include/astra/ParallelBeamLinearKernelProjector2D.inl
index 10d4892..53451e5 100644
--- a/include/astra/ParallelBeamLinearKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLinearKernelProjector2D.inl
@@ -154,18 +154,7 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom,
 
 		const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
 
-		float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY);
-
 		const bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY);
-		if (vertical) {
-			RxOverRy = proj->fRayX/proj->fRayY;
-			lengthPerRow = detSize * m_pVolumeGeometry->getPixelLengthX() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY);
-			deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX;
-		} else {
-			RyOverRx = proj->fRayY/proj->fRayX;
-			lengthPerCol = detSize * m_pVolumeGeometry->getPixelLengthY() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX);
-			deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
-		}
 
 		Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f;
 		Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
@@ -186,6 +175,10 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom,
 			// vertically
 			if (vertical) {
 
+				RxOverRy = proj->fRayX/proj->fRayY;
+				lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY);
+				deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX;
+
 				// calculate c for row 0
 				c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX;
 
@@ -209,6 +202,10 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom,
 			// horizontally
 			else {
 
+				RyOverRx = proj->fRayY/proj->fRayX;
+				lengthPerCol = m_pVolumeGeometry->getPixelLengthY() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX);
+				deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
+
 				// calculate r for col 0
 				r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY;
 
diff --git a/include/astra/ParallelBeamStripKernelProjector2D.inl b/include/astra/ParallelBeamStripKernelProjector2D.inl
index 0d775b3..2031560 100644
--- a/include/astra/ParallelBeamStripKernelProjector2D.inl
+++ b/include/astra/ParallelBeamStripKernelProjector2D.inl
@@ -142,20 +142,11 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom,
 
 		const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
 
+		const float32 rayWidth = fabs(proj->fDetUX * proj->fRayY - proj->fDetUY * proj->fRayX) /
+		                         sqrt(proj->fRayX * proj->fRayX + proj->fRayY * proj->fRayY);
+		const float32 relPixelArea = pixelArea / rayWidth;
+
 		bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY);
-		if (vertical) {
-			RxOverRy = proj->fRayX/proj->fRayY;
-			deltac = -m_pVolumeGeometry->getPixelLengthY() * RxOverRy * inv_pixelLengthX;
-			S = 0.5f - 0.5f*fabs(RxOverRy);
-			T = 0.5f + 0.5f*fabs(RxOverRy);
-			invTminS = 1.0f / (T-S);
-		} else {
-			RyOverRx = proj->fRayY/proj->fRayX;
-			deltar = -m_pVolumeGeometry->getPixelLengthX() * RyOverRx * inv_pixelLengthY;
-			S = 0.5f - 0.5f*fabs(RyOverRx);
-			T = 0.5f + 0.5f*fabs(RyOverRx);
-			invTminS = 1.0f / (T-S);
-		}
 
 		Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f;
 		Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
@@ -176,6 +167,12 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom,
 			// vertically
 			if (vertical) {
 
+				RxOverRy = proj->fRayX/proj->fRayY;
+				deltac = -m_pVolumeGeometry->getPixelLengthY() * RxOverRy * inv_pixelLengthX;
+				S = 0.5f - 0.5f*fabs(RxOverRy);
+				T = 0.5f + 0.5f*fabs(RxOverRy);
+				invTminS = 1.0f / (T-S);
+
 				// calculate cL and cR for row 0
 				cL = (DLx + (Ey - DLy)*RxOverRy - Ex) * inv_pixelLengthX;
 				cR = (DRx + (Ey - DRy)*RxOverRy - Ex) * inv_pixelLengthX;
@@ -219,7 +216,7 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom,
 							else if (-S < offsetL)  res -= 0.5f + offsetL;
 							else if (-T < offsetL)  res -= 0.5f*(offsetL+T)*(offsetL+T)*invTminS;
 
-							p.addWeight(iRayIndex, iVolumeIndex, pixelArea*res);
+							p.addWeight(iRayIndex, iVolumeIndex, relPixelArea*res);
 							p.pixelPosterior(iVolumeIndex);
 						}
 					}
@@ -229,6 +226,12 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom,
 			// horizontally
 			else {
 
+				RyOverRx = proj->fRayY/proj->fRayX;
+				deltar = -m_pVolumeGeometry->getPixelLengthX() * RyOverRx * inv_pixelLengthY;
+				S = 0.5f - 0.5f*fabs(RyOverRx);
+				T = 0.5f + 0.5f*fabs(RyOverRx);
+				invTminS = 1.0f / (T-S);
+
 				// calculate rL and rR for row 0
 				rL = -(DLy + (Ex - DLx)*RyOverRx - Ey) * inv_pixelLengthY;
 				rR = -(DRy + (Ex - DRx)*RyOverRx - Ey) * inv_pixelLengthY;
@@ -272,7 +275,7 @@ void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom,
 							else if (-S < offsetL)  res -= 0.5f + offsetL;
 							else if (-T < offsetL)  res -= 0.5f*(offsetL+T)*(offsetL+T)*invTminS;
 
-							p.addWeight(iRayIndex, iVolumeIndex, pixelArea*res);
+							p.addWeight(iRayIndex, iVolumeIndex, relPixelArea*res);
 							p.pixelPosterior(iVolumeIndex);
 						}
 					}
diff --git a/include/astra/cuda/2d/algo.h b/include/astra/cuda/2d/algo.h
index 3fccdef..b670b8b 100644
--- a/include/astra/cuda/2d/algo.h
+++ b/include/astra/cuda/2d/algo.h
@@ -56,6 +56,10 @@ public:
 
 	bool setSuperSampling(int raysPerDet, int raysPerPixelDim);
 
+	// Scale the final reconstruction.
+	// May be called at any time after setGeometry and before iterate(). Multiple calls stack.
+	bool setReconstructionScale(float fScale);
+
 	virtual bool enableVolumeMask();
 	virtual bool enableSinogramMask();
 
@@ -88,8 +92,7 @@ public:
 	// sinogram, reconstruction, volume mask and sinogram mask in system RAM,
 	// respectively. The corresponding pitch variables give the pitches
 	// of these buffers, measured in floats.
-	// The sinogram is multiplied by fSinogramScale after uploading it.
-	virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale,
+	virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch,
 	                           const float* pfReconstruction, unsigned int iReconstructionPitch,
 	                           const float* pfVolMask, unsigned int iVolMaskPitch,
 	                           const float* pfSinoMask, unsigned int iSinoMaskPitch);
@@ -133,7 +136,7 @@ protected:
 	SDimensions dims;
 	SParProjection* parProjs;
 	SFanProjection* fanProjs;
-	float fOutputScale;
+	float fProjectorScale;
 
 	bool freeGPUMemory;
 
diff --git a/include/astra/cuda/2d/cgls.h b/include/astra/cuda/2d/cgls.h
index 375a425..a854a74 100644
--- a/include/astra/cuda/2d/cgls.h
+++ b/include/astra/cuda/2d/cgls.h
@@ -47,7 +47,7 @@ public:
 	virtual bool setBuffers(float* D_volumeData, unsigned int volumePitch,
 	                        float* D_projData, unsigned int projPitch);
 
-	virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch, float fSinogramScale,
+	virtual bool copyDataToGPU(const float* pfSinogram, unsigned int iSinogramPitch,
 	                           const float* pfReconstruction, unsigned int iReconstructionPitch,
 	                           const float* pfVolMask, unsigned int iVolMaskPitch,
 	                           const float* pfSinoMask, unsigned int iSinoMaskPitch);
diff --git a/include/astra/cuda/2d/fbp.h b/include/astra/cuda/2d/fbp.h
index 1adf3b1..3aa4741 100644
--- a/include/astra/cuda/2d/fbp.h
+++ b/include/astra/cuda/2d/fbp.h
@@ -79,6 +79,11 @@ public:
 
 	bool setShortScan(bool ss) { m_bShortScan = ss; return true; }
 
+	// Scale the final reconstruction.
+	// May be called at any time before iterate().
+	bool setReconstructionScale(float fScale);
+
+
 	virtual bool init();
 
 	virtual bool iterate(unsigned int iterations);
@@ -90,6 +95,7 @@ protected:
 
 	void* D_filter; // cufftComplex*
 	bool m_bShortScan;
+	float fReconstructionScale;
 };
 
 }
diff --git a/include/astra/cuda/3d/fdk.h b/include/astra/cuda/3d/fdk.h
index 3b1a9e1..6817154 100644
--- a/include/astra/cuda/3d/fdk.h
+++ b/include/astra/cuda/3d/fdk.h
@@ -35,7 +35,7 @@ namespace astraCUDA3d {
 bool FDK_PreWeight(cudaPitchedPtr D_projData,
                 float fSrcOrigin, float fDetOrigin,
                 float fZShift,
-                float fDetUSize, float fDetVSize, float fVoxSize,
+                float fDetUSize, float fDetVSize,
                 bool bShortScan,
                 const SDimensions3D& dims, const float* angles);
 
diff --git a/include/astra/cuda/3d/mem3d.h b/include/astra/cuda/3d/mem3d.h
index 8c3956e..fff1490 100644
--- a/include/astra/cuda/3d/mem3d.h
+++ b/include/astra/cuda/3d/mem3d.h
@@ -110,7 +110,7 @@ bool copyIntoArray(MemHandle3D handle, MemHandle3D subdata, const SSubDimensions
 
 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, bool bFDKWeighting);
+bool BP(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, int iVoxelSuperSampling);
 
 bool FDK(const astra::CProjectionGeometry3D* pProjGeom, MemHandle3D projData, const astra::CVolumeGeometry3D* pVolGeom, MemHandle3D volData, bool bShortScan, const float *pfFilter = 0);
 
diff --git a/include/astra/cuda/3d/util3d.h b/include/astra/cuda/3d/util3d.h
index 0146d2d..200abfc 100644
--- a/include/astra/cuda/3d/util3d.h
+++ b/include/astra/cuda/3d/util3d.h
@@ -64,6 +64,53 @@ float dotProduct3D(cudaPitchedPtr data, unsigned int x, unsigned int y, unsigned
 
 int calcNextPowerOfTwo(int _iValue);
 
+struct Vec3 {
+	double x;
+	double y;
+	double z;
+	Vec3(double x_, double y_, double z_) : x(x_), y(y_), z(z_) { }
+	Vec3 operator+(const Vec3 &b) const {
+		return Vec3(x + b.x, y + b.y, z + b.z);
+	}
+	Vec3 operator-(const Vec3 &b) const {
+		return Vec3(x - b.x, y - b.y, z - b.z);
+	}
+	Vec3 operator-() const {
+		return Vec3(-x, -y, -z);
+	}
+	double norm() const {
+		return sqrt(x*x + y*y + z*z);
+	}
+};
+
+static double det3x(const Vec3 &b, const Vec3 &c) {
+	return (b.y * c.z - b.z * c.y);
+}
+static double det3y(const Vec3 &b, const Vec3 &c) {
+	return -(b.x * c.z - b.z * c.x);
+}
+
+static double det3z(const Vec3 &b, const Vec3 &c) {
+	return (b.x * c.y - b.y * c.x);
+}
+
+static double det3(const Vec3 &a, const Vec3 &b, const Vec3 &c) {
+	return a.x * det3x(b,c) + a.y * det3y(b,c) + a.z * det3z(b,c);
+}
+
+static Vec3 cross3(const Vec3 &a, const Vec3 &b) {
+	return Vec3(det3x(a,b), det3y(a,b), det3z(a,b));
+}
+
+static Vec3 scaled_cross3(const Vec3 &a, const Vec3 &b, const Vec3 &sc) {
+	Vec3 ret = cross3(a, b);
+	ret.x *= sc.y * sc.z;
+	ret.y *= sc.x * sc.z;
+	ret.z *= sc.x * sc.y;
+	return ret;
+}
+
+
 }
 
 #endif
diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp
index 1319a87..822f746 100644
--- a/src/CompositeGeometryManager.cpp
+++ b/src/CompositeGeometryManager.cpp
@@ -1462,12 +1462,10 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
 		Cuda3DProjectionKernel projKernel = ker3d_default;
 		int detectorSuperSampling = 1;
 		int voxelSuperSampling = 1;
-		bool densityWeighting = false;
 		if (projector) {
 			projKernel = projector->getProjectionKernel();
 			detectorSuperSampling = projector->getDetectorSuperSampling();
 			voxelSuperSampling = projector->getVoxelSuperSampling();
-			densityWeighting = projector->getDensityWeighting();
 		}
 
 		size_t inx, iny, inz;
@@ -1513,7 +1511,7 @@ static bool doJob(const CCompositeGeometryManager::TJobSet::const_iterator& iter
 
 			ASTRA_DEBUG("CCompositeGeometryManager::doJobs: doing BP");
 
-			ok = astraCUDA3d::BP(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, srcMem->hnd, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, dstMem->hnd, voxelSuperSampling, densityWeighting);
+			ok = astraCUDA3d::BP(((CCompositeGeometryManager::CProjectionPart*)j.pInput.get())->pGeom, srcMem->hnd, ((CCompositeGeometryManager::CVolumePart*)j.pOutput.get())->pGeom, dstMem->hnd, voxelSuperSampling);
 			if (!ok) ASTRA_ERROR("Error performing sub-BP");
 			ASTRA_DEBUG("CCompositeGeometryManager::doJobs: BP done");
 		}
diff --git a/src/CudaFilteredBackProjectionAlgorithm.cpp b/src/CudaFilteredBackProjectionAlgorithm.cpp
index 88e235b..c1d3dc8 100644
--- a/src/CudaFilteredBackProjectionAlgorithm.cpp
+++ b/src/CudaFilteredBackProjectionAlgorithm.cpp
@@ -151,6 +151,13 @@ void CCudaFilteredBackProjectionAlgorithm::initCUDAAlgorithm()
 	if (!ok) {
 		ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set short-scan mode");
 	}
+
+	const CVolumeGeometry2D& volGeom = *m_pReconstruction->getGeometry();
+	float fPixelArea = volGeom.getPixelArea();
+	ok &= pFBP->setReconstructionScale(1.0f/fPixelArea);
+	if (!ok) {
+		ASTRA_ERROR("CCudaFilteredBackProjectionAlgorithm: Failed to set reconstruction scale");
+	}
 }
 
 
diff --git a/src/CudaProjector3D.cpp b/src/CudaProjector3D.cpp
index 3ea7043..e5c55cc 100644
--- a/src/CudaProjector3D.cpp
+++ b/src/CudaProjector3D.cpp
@@ -67,7 +67,6 @@ void CCudaProjector3D::_clear()
 	m_iVoxelSuperSampling = 1;
 	m_iDetectorSuperSampling = 1;
 	m_iGPUIndex = -1;
-	m_bDensityWeighting = false;
 }
 
 //----------------------------------------------------------------------------------------
@@ -132,13 +131,6 @@ bool CCudaProjector3D::initialize(const Config& _cfg)
 	m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", 1);
 	CC.markOptionParsed("DetectorSuperSampling");
 
-	if (dynamic_cast<CConeProjectionGeometry3D*>(m_pProjectionGeometry) ||
-	    dynamic_cast<CConeVecProjectionGeometry3D*>(m_pProjectionGeometry))
-	{
-		m_bDensityWeighting = _cfg.self.getOptionBool("DensityWeighting", false);
-		CC.markOptionParsed("DensityWeighting");
-	}
-
 	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", -1);
 	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
 	CC.markOptionParsed("GPUIndex");
diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp
index 1e81390..6730cea 100644
--- a/src/CudaReconstructionAlgorithm2D.cpp
+++ b/src/CudaReconstructionAlgorithm2D.cpp
@@ -309,10 +309,7 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations)
 		m_bAlgoInit = true;
 	}
 
-	float fPixelSize = volgeom.getPixelLengthX();
-	float fSinogramScale = 1.0f/(fPixelSize*fPixelSize);
-
-	ok = m_pAlgo->copyDataToGPU(m_pSinogram->getDataConst(), m_pSinogram->getGeometry()->getDetectorCount(), fSinogramScale,
+	ok = m_pAlgo->copyDataToGPU(m_pSinogram->getDataConst(), m_pSinogram->getGeometry()->getDetectorCount(),
 	                            m_pReconstruction->getDataConst(), volgeom.getGridColCount(),
 	                            m_bUseReconstructionMask ? m_pReconstructionMask->getDataConst() : 0, volgeom.getGridColCount(),
 	                            m_bUseSinogramMask ? m_pSinogramMask->getDataConst() : 0, m_pSinogram->getGeometry()->getDetectorCount());
diff --git a/src/Features.cpp b/src/Features.cpp
index 9114131..09a3499 100644
--- a/src/Features.cpp
+++ b/src/Features.cpp
@@ -34,6 +34,12 @@ _AstraExport bool hasFeature(const std::string &flag) {
 	if (flag == "cuda") {
 		return cudaEnabled();
 	}
+	if (flag == "projectors_scaled_as_line_integrals") {
+		return true;
+	}
+	if (flag == "fan_cone_BP_density_weighting_by_default") {
+		return true;
+	}
 
 	return false;
 }
diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp
index 423dc6c..6b4093d 100644
--- a/src/FilteredBackProjectionAlgorithm.cpp
+++ b/src/FilteredBackProjectionAlgorithm.cpp
@@ -167,6 +167,11 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
 
 	m_filterConfig = getFilterConfigForAlgorithm(_cfg, this);
 
+	const CParallelProjectionGeometry2D* parprojgeom = dynamic_cast<CParallelProjectionGeometry2D*>(m_pSinogram->getGeometry());
+	if (!parprojgeom) {
+		ASTRA_ERROR("FBP currently only supports parallel projection geometries.");
+		return false;
+	}
 
 	// TODO: check that the angles are linearly spaced between 0 and pi
 
@@ -264,8 +269,12 @@ void CFilteredBackProjectionAlgorithm::run(int _iNrIterations)
 	            DefaultBPPolicy(m_pReconstruction, &filteredSinogram));
 
 	// Scale data
-	int iAngleCount = m_pProjector->getProjectionGeometry()->getProjectionAngleCount();
-	(*m_pReconstruction) *= (PI/2)/iAngleCount;
+	const CVolumeGeometry2D& volGeom = *m_pProjector->getVolumeGeometry();
+	const CProjectionGeometry2D& projGeom = *m_pProjector->getProjectionGeometry();
+
+	int iAngleCount = projGeom.getProjectionAngleCount();
+	float fPixelArea = volGeom.getPixelArea();
+	(*m_pReconstruction) *= PI/(2*iAngleCount*fPixelArea);
 
 	m_pReconstruction->updateStatistics();
 }
diff --git a/src/GeometryUtil2D.cpp b/src/GeometryUtil2D.cpp
index e09a3bc..806572f 100644
--- a/src/GeometryUtil2D.cpp
+++ b/src/GeometryUtil2D.cpp
@@ -28,6 +28,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/GeometryUtil2D.h"
 
 #include <cmath>
+#include <cstdio>
 
 namespace astra {
 
@@ -158,14 +159,16 @@ bool getFanParameters(const SFanProjection &proj, unsigned int iProjDets, float
 	// project origin on detector line ( == project source on detector line)
 
 	double t = (- proj.fDetSX) * proj.fDetUX + (- proj.fDetSY) * proj.fDetUY;
+	t /= (proj.fDetUX * proj.fDetUX + proj.fDetUY * proj.fDetUY);
 
 	fOffset = (float)t - 0.5*iProjDets;
 
-	// TODO: CHECKME
 	fOriginDetector = sqrt((proj.fDetSX + t * proj.fDetUX)*(proj.fDetSX + t * proj.fDetUX) + (proj.fDetSY + t * proj.fDetUY)*(proj.fDetSY + t * proj.fDetUY));
 
-	//float fAngle = atan2(proj.fDetSX + t * proj.fDetUX - proj.fSrcX, proj.fDetSY + t * proj.fDetUY); // TODO: Fix order + sign
-	fAngle = atan2(proj.fDetUY, proj.fDetUX); // TODO: Check order + sign
+	fAngle = atan2(proj.fDetUY, proj.fDetUX);
+
+	//fprintf(stderr, "getFanParams: s = (%f,%f) d = (%f,%f) u = (%f,%f)\n", proj.fSrcX, proj.fSrcY, proj.fDetSX, proj.fDetSY, proj.fDetUX, proj.fDetUY);
+	//fprintf(stderr, "getFanParams: fOS = %f, fOD = %f, detsize = %f, offset = %f (t = %f), angle = %f\n", fOriginSource, fOriginDetector, fDetSize, fOffset, t, fAngle);
 
 	return true;
 }
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
index d04ffb8..e5d8f2b 100644
--- a/tests/python/test_line2d.py
+++ b/tests/python/test_line2d.py
@@ -7,9 +7,9 @@ import pylab
 # Display sinograms with mismatch on test failure
 DISPLAY=False
 
-NONUNITDET=False
-OBLIQUE=False
-FLEXVOL=False
+NONUNITDET=True
+OBLIQUE=True
+FLEXVOL=True
 NONSQUARE=False  # non-square pixels not supported yet by most projectors
 
 # Round interpolation weight to 8 bits to emulate CUDA texture unit precision
@@ -20,15 +20,8 @@ nloops = 50
 seed = 123
 
 
-# FAILURES:
-# fan/cuda with flexible volume
-# detweight for fan/cuda
-# fan/strip relatively high numerical errors?
-# parvec/line+linear for oblique
-
-# INCONSISTENCY:
-# effective_detweight vs norm(detu) in line/linear (oblique)
-
+# KNOWN FAILURES:
+# fan/strip relatively high numerical errors around 45 degrees
 
 
 # return length of intersection of the line through points src = (x,y)
@@ -454,23 +447,15 @@ class Test2DKernel(unittest.TestCase):
         for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
           (src, det) = center
 
-          try:
-            detweight = pg['DetectorWidth']
-          except KeyError:
-            if 'fan' not in type:
-              detweight = effective_detweight(src, det, pg['Vectors'][i//pg['DetectorCount'],4:6])
-            else:
-              detweight = np.linalg.norm(pg['Vectors'][i//pg['DetectorCount'],4:6], ord=2)
-
           # We compute line intersections with slightly bigger (cw) and
           # smaller (aw) rectangles, and see if the kernel falls
           # between these two values.
           (aw,bw,cw) = intersect_line_rectangle_interval(src, det,
                         xmin, xmax, ymin, ymax,
                         1e-3)
-          a[i] = aw * detweight
-          b[i] = bw * detweight
-          c[i] = cw * detweight
+          a[i] = aw
+          b[i] = bw
+          c[i] = cw
         a = a.reshape(astra.functions.geom_size(pg))
         b = b.reshape(astra.functions.geom_size(pg))
         c = c.reshape(astra.functions.geom_size(pg))
@@ -494,17 +479,9 @@ class Test2DKernel(unittest.TestCase):
         for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
           (src, det) = center
           (xd, yd) = det - src
-          try:
-            detweight = pg['DetectorWidth']
-          except KeyError:
-            if 'fan' not in type:
-              detweight = effective_detweight(src, det, pg['Vectors'][i//pg['DetectorCount'],4:6])
-            else:
-              detweight = np.linalg.norm(pg['Vectors'][i//pg['DetectorCount'],4:6], ord=2)
-
           l = 0.0
           if np.abs(xd) > np.abs(yd): # horizontal ray
-            length = math.sqrt(1.0 + abs(yd/xd)**2)
+            length = math.sqrt(1.0 + abs(yd/xd)**2) * pixsize[0]
             y_seg = (ymin, ymax)
             for j in range(rect_min[0], rect_max[0]):
               x = origin[0] + (-0.5 * shape[0] + j + 0.5) * pixsize[0]
@@ -512,9 +489,9 @@ class Test2DKernel(unittest.TestCase):
               # limited interpolation precision with cuda
               if CUDA_8BIT_LINEAR and proj_type == 'cuda':
                 w = np.round(w * 256.0) / 256.0
-              l += w * length * pixsize[0] * detweight
+              l += w * length
           else:
-            length = math.sqrt(1.0 + abs(xd/yd)**2)
+            length = math.sqrt(1.0 + abs(xd/yd)**2) * pixsize[1]
             x_seg = (xmin, xmax)
             for j in range(rect_min[1], rect_max[1]):
               y = origin[1] + (+0.5 * shape[1] - j - 0.5) * pixsize[1]
@@ -522,7 +499,7 @@ class Test2DKernel(unittest.TestCase):
               # limited interpolation precision with cuda
               if CUDA_8BIT_LINEAR and proj_type == 'cuda':
                 w = np.round(w * 256.0) / 256.0
-              l += w * length * pixsize[1] * detweight
+              l += w * length
           a[i] = l
         a = a.reshape(astra.functions.geom_size(pg))
         if not np.all(np.isfinite(a)):
@@ -532,21 +509,26 @@ class Test2DKernel(unittest.TestCase):
         if DISPLAY and x > TOL:
           display_mismatch(data, sinogram, a)
         self.assertFalse(x > TOL)
-      elif proj_type == 'distance_driven':
+      elif proj_type == 'distance_driven' and 'par' in type:
         a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32)
         for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
-          (xd, yd) = center[1] - center[0]
+          (src, det) = center
+          try:
+            detweight = pg['DetectorWidth']
+          except KeyError:
+            detweight = effective_detweight(src, det, pg['Vectors'][i//pg['DetectorCount'],4:6])
+          (xd, yd) = det - src
           l = 0.0
           if np.abs(xd) > np.abs(yd): # horizontal ray
             y_seg = (ymin, ymax)
             for j in range(rect_min[0], rect_max[0]):
               x = origin[0] + (-0.5 * shape[0] + j + 0.5) * pixsize[0]
-              l += intersect_ray_vertical_segment(edge1, edge2, x, y_seg) * pixsize[0]
+              l += intersect_ray_vertical_segment(edge1, edge2, x, y_seg) * pixsize[0] / detweight
           else:
             x_seg = (xmin, xmax)
             for j in range(rect_min[1], rect_max[1]):
               y = origin[1] + (+0.5 * shape[1] - j - 0.5) * pixsize[1]
-              l += intersect_ray_horizontal_segment(edge1, edge2, y, x_seg) * pixsize[1]
+              l += intersect_ray_horizontal_segment(edge1, edge2, y, x_seg) * pixsize[1] / detweight
           a[i] = l
         a = a.reshape(astra.functions.geom_size(pg))
         if not np.all(np.isfinite(a)):
@@ -560,6 +542,7 @@ class Test2DKernel(unittest.TestCase):
         a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32)
         for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
           (src, det) = center
+          detweight = effective_detweight(src, det, edge2[1] - edge1[1])
           det_dist = np.linalg.norm(src-det, ord=2)
           l = 0.0
           for j in range(rect_min[0], rect_max[0]):
@@ -570,7 +553,7 @@ class Test2DKernel(unittest.TestCase):
               ymin = origin[1] + (+0.5 * shape[1] - k - 1) * pixsize[1]
               ymax = origin[1] + (+0.5 * shape[1] - k) * pixsize[1]
               ycen = 0.5 * (ymin + ymax)
-              scale = det_dist / np.linalg.norm( src - np.array((xcen,ycen)), ord=2 )
+              scale = det_dist / (np.linalg.norm( src - np.array((xcen,ycen)), ord=2 ) * detweight)
               w = intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax)
               l += w * scale
           a[i] = l
@@ -578,14 +561,20 @@ class Test2DKernel(unittest.TestCase):
         if not np.all(np.isfinite(a)):
           raise RuntimeError("Invalid value in reference sinogram")
         x = np.max(np.abs(sinogram-a))
-        TOL = 8e-3
+        # BUG: Known bug in fan/strip code around 45 degree projections causing larger errors than desirable
+        TOL = 4e-2
         if DISPLAY and x > TOL:
           display_mismatch(data, sinogram, a)
         self.assertFalse(x > TOL)
       elif proj_type == 'strip':
         a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32)
         for i, (center, edge1, edge2) in enumerate(gen_lines(pg)):
-          a[i] = intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax)
+          (src, det) = center
+          try:
+            detweight = pg['DetectorWidth']
+          except KeyError:
+            detweight = effective_detweight(src, det, pg['Vectors'][i//pg['DetectorCount'],4:6])
+          a[i] = intersect_ray_rect(edge1, edge2, xmin, xmax, ymin, ymax) / detweight
         a = a.reshape(astra.functions.geom_size(pg))
         if not np.all(np.isfinite(a)):
           raise RuntimeError("Invalid value in reference sinogram")
@@ -594,46 +583,83 @@ class Test2DKernel(unittest.TestCase):
         if DISPLAY and x > TOL:
           display_mismatch(data, sinogram, a)
         self.assertFalse(x > TOL)
+      else:
+        raise RuntimeError("Unsupported projector")
 
-  def multi_test(self, type, proj_type):
-    np.random.seed(seed)
-    for _ in range(nloops):
-      self.single_test(type, proj_type)
-
-  def test_par(self):
-    self.multi_test('parallel', 'line')
-  def test_par_linear(self):
-    self.multi_test('parallel', 'linear')
-  def test_par_cuda(self):
-    self.multi_test('parallel', 'cuda')
-  def test_par_dd(self):
-    self.multi_test('parallel', 'distance_driven')
-  def test_par_strip(self):
-    self.multi_test('parallel', 'strip')
-  def test_fan(self):
-    self.multi_test('fanflat', 'line')
-  def test_fan_strip(self):
-    self.multi_test('fanflat', 'strip')
-  def test_fan_cuda(self):
-    self.multi_test('fanflat', 'cuda')
-  def test_parvec(self):
-    self.multi_test('parallel_vec', 'line')
-  def test_parvec_linear(self):
-    self.multi_test('parallel_vec', 'linear')
-  def test_parvec_dd(self):
-    self.multi_test('parallel_vec', 'distance_driven')
-  def test_parvec_strip(self):
-    self.multi_test('parallel_vec', 'strip')
-  def test_parvec_cuda(self):
-    self.multi_test('parallel_vec', 'cuda')
-  def test_fanvec(self):
-    self.multi_test('fanflat_vec', 'line')
-  def test_fanvec_cuda(self):
-    self.multi_test('fanflat_vec', 'cuda')
+  def single_test_adjoint(self, type, proj_type):
+      shape = np.random.randint(*range2d, size=2)
+      if FLEXVOL:
+          if not NONSQUARE:
+            pixsize = np.array([0.5, 0.5]) + np.random.random()
+          else:
+            pixsize = 0.5 + np.random.random(size=2)
+          origin = 10 * np.random.random(size=2)
+      else:
+          pixsize = (1.,1.)
+          origin = (0.,0.)
+      vg = astra.create_vol_geom(shape[1], shape[0],
+                                 origin[0] - 0.5 * shape[0] * pixsize[0],
+                                 origin[0] + 0.5 * shape[0] * pixsize[0],
+                                 origin[1] - 0.5 * shape[1] * pixsize[1],
+                                 origin[1] + 0.5 * shape[1] * pixsize[1])
 
+      if type == 'parallel':
+        pg = gen_random_geometry_parallel()
+        projector_id = astra.create_projector(proj_type, pg, vg)
+      elif type == 'parallel_vec':
+        pg = gen_random_geometry_parallel_vec()
+        projector_id = astra.create_projector(proj_type, pg, vg)
+      elif type == 'fanflat':
+        pg = gen_random_geometry_fanflat()
+        projector_id = astra.create_projector(proj_type_to_fan(proj_type), pg, vg)
+      elif type == 'fanflat_vec':
+        pg = gen_random_geometry_fanflat_vec()
+        projector_id = astra.create_projector(proj_type_to_fan(proj_type), pg, vg)
 
+      for i in range(5):
+        X = np.random.random((shape[1], shape[0]))
+        Y = np.random.random(astra.geom_size(pg))
+
+        sinogram_id, fX = astra.create_sino(X, projector_id)
+        bp_id, fTY = astra.create_backprojection(Y, projector_id)
+
+        astra.data2d.delete(sinogram_id)
+        astra.data2d.delete(bp_id)
+
+        da = np.dot(fX.ravel(), Y.ravel())
+        db = np.dot(X.ravel(), fTY.ravel())
+        m = np.abs(da - db)
+        TOL = 1e-3 if 'cuda' not in proj_type else 1e-1
+        if m / da >= TOL:
+          print(vg)
+          print(pg)
+          print(m/da, da/db, da, db)
+        self.assertTrue(m / da < TOL)
+      astra.projector.delete(projector_id)
 
 
+  def multi_test(self, type, proj_type):
+    np.random.seed(seed)
+    for _ in range(nloops):
+      self.single_test(type, proj_type)
+  def multi_test_adjoint(self, type, proj_type):
+    np.random.seed(seed)
+    for _ in range(nloops):
+      self.single_test_adjoint(type, proj_type)
+
+__combinations = { 'parallel': [ 'line', 'linear', 'distance_driven', 'strip', 'cuda' ],
+                   'parallel_vec': [ 'line', 'linear', 'distance_driven', 'strip', 'cuda' ],
+                   'fanflat': [ 'line', 'strip', 'cuda' ],
+                   'fanflat_vec': [ 'line', 'cuda' ] }
+
+for k, l in __combinations.items():
+  for v in l:
+    def f(k,v):
+      return lambda self: self.multi_test(k, v)
+    def f_adj(k,v):
+      return lambda self: self.multi_test_adjoint(k, v)
+    setattr(Test2DKernel, 'test_' + k + '_' + v, f(k,v))
+    setattr(Test2DKernel, 'test_' + k + '_' + v + '_adjoint', f_adj(k,v))
 
 if __name__ == '__main__':
   unittest.main()
diff --git a/tests/python/test_rec_scaling.py b/tests/python/test_rec_scaling.py
new file mode 100644
index 0000000..621fd8a
--- /dev/null
+++ b/tests/python/test_rec_scaling.py
@@ -0,0 +1,213 @@
+import numpy as np
+import unittest
+import astra
+import math
+import pylab
+
+DISPLAY=False
+
+def VolumeGeometries(is3D,noncube):
+  if not is3D:
+    for s in [0.8, 1.0, 1.25]:
+      yield astra.create_vol_geom(128, 128, -64*s, 64*s, -64*s, 64*s)
+  elif noncube:
+    for sx in [0.8, 1.0]:
+      for sy in [0.8, 1.0]:
+        for sz in [0.8, 1.0]:
+          yield astra.create_vol_geom(64, 64, 64, -32*sx, 32*sx, -32*sy, 32*sy, -32*sz, 32*sz)
+  else:
+    for s in [0.8, 1.0]:
+      yield astra.create_vol_geom(64, 64, 64, -32*s, 32*s, -32*s, 32*s, -32*s, 32*s)
+
+
+def ProjectionGeometries(type):
+  if type == 'parallel':
+    for dU in [0.8, 1.0, 1.25]:
+      yield astra.create_proj_geom('parallel', dU, 256, np.linspace(0,np.pi,180,False))
+  elif type == 'fanflat':
+    for dU in [0.8, 1.0, 1.25]:
+      for src in [500, 1000]:
+        for det in [0, 250, 500]:
+          yield astra.create_proj_geom('fanflat', dU, 256, np.linspace(0,2*np.pi,180,False), src, det)
+  elif type == 'parallel3d':
+    for dU in [0.8, 1.0]:
+      for dV in [0.8, 1.0]:
+        yield astra.create_proj_geom('parallel3d', dU, dV, 128, 128, np.linspace(0,np.pi,180,False))
+  elif type == 'parallel3d_vec':
+    for j in range(10):
+       Vectors = np.zeros([180,12])
+       wu = 0.6 + 0.8 * np.random.random()
+       wv = 0.6 + 0.8 * np.random.random()
+       for i in range(Vectors.shape[0]):
+         l = 0.6 + 0.8 * np.random.random()
+         angle1 = 2*np.pi*np.random.random()
+         angle2 = angle1 + 0.5 * np.random.random()
+         angle3 = 0.1*np.pi*np.random.random()
+         detc = 10 * np.random.random(size=3)
+         detu = [ math.cos(angle1) * wu, math.sin(angle1) * wu, 0 ]
+         detv = [ -math.sin(angle1) * math.sin(angle3) * wv, math.cos(angle1) * math.sin(angle3) * wv, math.cos(angle3) * wv ]
+         ray = [ math.sin(angle2) * l, -math.cos(angle2) * l, 0 ]
+         Vectors[i, :] = [ ray[0], ray[1], ray[2], detc[0], detc[1], detc[2], detu[0], detu[1], detu[2], detv[0], detv[1], detv[2] ]
+       pg = astra.create_proj_geom('parallel3d_vec', 128, 128, Vectors)
+       yield pg
+  elif type == 'cone':
+    for dU in [0.8, 1.0]:
+      for dV in [0.8, 1.0]:
+        for src in [500, 1000]:
+          for det in [0, 250]:
+            yield astra.create_proj_geom('cone', dU, dV, 128, 128, np.linspace(0,2*np.pi,180,False), src, det)
+  elif type == 'cone_vec':
+    for j in range(10):
+       Vectors = np.zeros([180,12])
+       wu = 0.6 + 0.8 * np.random.random()
+       wv = 0.6 + 0.8 * np.random.random()
+       for i in range(Vectors.shape[0]):
+         l = 256 * (0.5 * np.random.random())
+         angle1 = 2*np.pi*np.random.random()
+         angle2 = angle1 + 0.5 * np.random.random()
+         angle3 = 0.1*np.pi*np.random.random()
+         detc = 10 * np.random.random(size=3)
+         detu = [ math.cos(angle1) * wu, math.sin(angle1) * wu, 0 ]
+         detv = [ -math.sin(angle1) * math.sin(angle3) * wv, math.cos(angle1) * math.sin(angle3) * wv, math.cos(angle3) * wv ]
+         src = [ math.sin(angle2) * l, -math.cos(angle2) * l, 0 ]
+         Vectors[i, :] = [ src[0], src[1], src[2], detc[0], detc[1], detc[2], detu[0], detu[1], detu[2], detv[0], detv[1], detv[2] ]
+       pg = astra.create_proj_geom('parallel3d_vec', 128, 128, Vectors)
+       yield pg
+
+
+class TestRecScale(unittest.TestCase):
+  def single_test(self, geom_type, proj_type, alg, iters):
+    if alg == 'FBP' and 'fanflat' in geom_type:
+      self.skipTest('CPU FBP is parallel-beam only')
+    is3D = (geom_type in ['parallel3d', 'cone'])
+    for vg in VolumeGeometries(is3D, 'FDK' not in alg):
+      for pg in ProjectionGeometries(geom_type):
+        if not is3D:
+          vol = np.zeros((128,128),dtype=np.float32)
+          vol[50:70,50:70] = 1
+        else:
+          vol = np.zeros((64,64,64),dtype=np.float32)
+          vol[25:35,25:35,25:35] = 1
+        proj_id = astra.create_projector(proj_type, pg, vg)
+        if not is3D:
+          sino_id, sinogram = astra.create_sino(vol, proj_id)
+        else:
+          sino_id, sinogram = astra.create_sino3d_gpu(vol, pg, vg)
+        if not is3D:
+          DATA = astra.data2d
+        else:
+          DATA = astra.data3d
+
+        rec_id = DATA.create('-vol', vg, 0.0 if 'EM' not in alg else 1.0)
+
+        cfg = astra.astra_dict(alg)
+        cfg['ReconstructionDataId'] = rec_id
+        cfg['ProjectionDataId'] = sino_id
+        cfg['ProjectorId'] = proj_id
+        alg_id = astra.algorithm.create(cfg)
+
+        for i in range(iters):
+            astra.algorithm.run(alg_id, 1)
+        rec = DATA.get(rec_id)
+        astra.astra.delete([sino_id, alg_id, alg_id, proj_id])
+        if not is3D:
+          val = np.sum(rec[55:65,55:65]) / 100.
+        else:
+          val = np.sum(rec[27:32,27:32,27:32]) / 125.
+        TOL = 5e-2
+        if DISPLAY and abs(val-1.0) >= TOL:
+          print(geom_type, proj_type, alg, vg, pg)
+          print(val)
+          pylab.gray()
+          if not is3D:
+            pylab.imshow(rec)
+          else:
+            pylab.imshow(rec[:,32,:])
+          pylab.show()
+        self.assertTrue(abs(val-1.0) < TOL)
+
+  def single_test_adjoint3D(self, geom_type, proj_type):
+    for vg in VolumeGeometries(True, True):
+      for pg in ProjectionGeometries(geom_type):
+        for i in range(5):
+          X = np.random.random(astra.geom_size(vg))
+          Y = np.random.random(astra.geom_size(pg))
+          proj_id, fX = astra.create_sino3d_gpu(X, pg, vg)
+          bp_id, fTY = astra.create_backprojection3d_gpu(Y, pg, vg)
+
+          astra.data3d.delete([proj_id, bp_id])
+
+          da = np.dot(fX.ravel(), Y.ravel())
+          db = np.dot(X.ravel(), fTY.ravel())
+          m = np.abs(da - db)
+          TOL = 1e-1
+          if m / da >= TOL:
+            print(vg)
+            print(pg)
+            print(m/da, da/db, da, db)
+          self.assertTrue(m / da < TOL)
+
+
+
+
+
+__combinations = {
+  'parallel': [ 'line', 'linear', 'distance_driven', 'strip', 'cuda' ],
+  'fanflat': [ 'line_fanflat', 'strip_fanflat', 'cuda' ],
+  'parallel3d': [ 'cuda3d' ],
+  'cone': [ 'cuda3d' ],
+}
+
+__combinations_adjoint = {
+  'parallel3d': [ 'cuda3d' ],
+  'cone': [ 'cuda3d' ],
+  'parallel3d_vec': [ 'cuda3d' ],
+  'cone_vec': [ 'cuda3d' ],
+}
+
+__algs = {
+   'SIRT': 50, 'SART': 10*180, 'CGLS': 30,
+   'FBP': 1
+}
+
+__algs_CUDA = {
+  'SIRT_CUDA': 50, 'SART_CUDA': 10*180, 'CGLS_CUDA': 30, 'EM_CUDA': 50,
+  'FBP_CUDA': 1
+}
+
+__algs_parallel3d = {
+  'SIRT3D_CUDA': 200, 'CGLS3D_CUDA': 20,
+}
+
+__algs_cone = {
+  'SIRT3D_CUDA': 200, 'CGLS3D_CUDA': 20,
+  'FDK_CUDA': 1
+}
+
+
+
+for k, l in __combinations.items():
+  for v in l:
+    is3D = (k in ['parallel3d', 'cone'])
+    if k == 'parallel3d':
+      A = __algs_parallel3d
+    elif k == 'cone':
+      A = __algs_cone
+    elif v == 'cuda':
+      A = __algs_CUDA
+    else:
+      A = __algs
+    for a, i in A.items():
+      def f(k, v, a, i):
+        return lambda self: self.single_test(k, v, a, i)
+      setattr(TestRecScale, 'test_' + a + '_' + k + '_' + v, f(k,v,a,i))
+
+for k, l in __combinations_adjoint.items():
+  for v in l:
+    def g(k, v):
+      return lambda self: self.single_test_adjoint3D(k, v)
+    setattr(TestRecScale, 'test_adjoint_' + k + '_' + v, g(k,v))
+
+if __name__ == '__main__':
+  unittest.main()
+
diff --git a/tests/test_ParallelBeamLineKernelProjector2D.cpp b/tests/test_ParallelBeamLineKernelProjector2D.cpp
deleted file mode 100644
index 71130c1..0000000
--- a/tests/test_ParallelBeamLineKernelProjector2D.cpp
+++ /dev/null
@@ -1,82 +0,0 @@
-/*
------------------------------------------------------------------------
-Copyright: 2010-2018, imec Vision Lab, University of Antwerp
-           2014-2018, CWI, Amsterdam
-
-Contact: astra@astra-toolbox.com
-Website: http://www.astra-toolbox.com/
-
-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/>.
-
------------------------------------------------------------------------
-*/
-
-
-#define BOOST_TEST_DYN_LINK
-#include <boost/test/unit_test.hpp>
-#include <boost/test/auto_unit_test.hpp>
-#include <boost/test/floating_point_comparison.hpp>
-
-#include "astra/ParallelBeamLineKernelProjector2D.h"
-#include "astra/ParallelProjectionGeometry2D.h"
-#include "astra/VolumeGeometry2D.h"
-
-struct TestParallelBeamLineKernelProjector2D {
-        TestParallelBeamLineKernelProjector2D()
-	{
-		astra::float32 angles[] = { 1.0f };
-		BOOST_REQUIRE( projGeom.initialize(1, 9, 1.0f, angles) );
-		BOOST_REQUIRE( volGeom.initialize(6, 4) );
-		BOOST_REQUIRE( proj.initialize(&projGeom, &volGeom) );
-	}
-        ~TestParallelBeamLineKernelProjector2D()
-	{
-
-	}
-
-	astra::CParallelBeamLineKernelProjector2D proj;
-	astra::CParallelProjectionGeometry2D projGeom;
-	astra::CVolumeGeometry2D volGeom;
-};
-
-BOOST_FIXTURE_TEST_CASE( testParallelBeamLineKernelProjector2D_General, TestParallelBeamLineKernelProjector2D )
-{
-
-}
-
-BOOST_FIXTURE_TEST_CASE( testParallelBeamLineKernelProjector2D_Rectangle, TestParallelBeamLineKernelProjector2D )
-{
-	int iMax = proj.getProjectionWeightsCount(0);
-	BOOST_REQUIRE(iMax > 0);
-
-	astra::SPixelWeight* pPix = new astra::SPixelWeight[iMax];
-	BOOST_REQUIRE(pPix);
-
-	int iCount;
-	proj.computeSingleRayWeights(0, 4, pPix, iMax, iCount); 
-	BOOST_REQUIRE(iCount <= iMax);
-
-	astra::float32 fWeight = 0;
-	for (int i = 0; i < iCount; ++i)
-		fWeight += pPix[i].m_fWeight;
-
-	BOOST_CHECK_SMALL(fWeight - 7.13037f, 0.00001f); // 6 / sin(1)
-
-	delete[] pPix;
-}
-
-
diff --git a/tests/test_ParallelBeamLinearKernelProjector2D.cpp b/tests/test_ParallelBeamLinearKernelProjector2D.cpp
deleted file mode 100644
index 4610aa5..0000000
--- a/tests/test_ParallelBeamLinearKernelProjector2D.cpp
+++ /dev/null
@@ -1,170 +0,0 @@
-/*
------------------------------------------------------------------------
-Copyright: 2010-2018, imec Vision Lab, University of Antwerp
-           2014-2018, CWI, Amsterdam
-
-Contact: astra@astra-toolbox.com
-Website: http://www.astra-toolbox.com/
-
-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/>.
-
------------------------------------------------------------------------
-*/
-
-
-#define BOOST_TEST_DYN_LINK
-#include <boost/test/unit_test.hpp>
-#include <boost/test/auto_unit_test.hpp>
-#include <boost/test/floating_point_comparison.hpp>
-
-#include "astra/ParallelBeamLineKernelProjector2D.h"
-#include "astra/ParallelBeamLinearKernelProjector2D.h"
-#include "astra/ParallelBeamStripKernelProjector2D.h"
-#include "astra/ParallelProjectionGeometry2D.h"
-#include "astra/VolumeGeometry2D.h"
-#include "astra/ProjectionGeometry2D.h"
-
-#include <ctime>
-
-using astra::float32;
-
-struct TestParallelBeamLinearKernelProjector2D {
-        TestParallelBeamLinearKernelProjector2D()
-	{
-		astra::float32 angles[] = { 2.6f };
-		BOOST_REQUIRE( projGeom.initialize(1, 3, 1.0f, angles) );
-		BOOST_REQUIRE( volGeom.initialize(3, 2) );
-		BOOST_REQUIRE( proj.initialize(&projGeom, &volGeom) );
-	}
-        ~TestParallelBeamLinearKernelProjector2D()
-	{
-
-	}
-
-	astra::CParallelBeamLinearKernelProjector2D proj;
-	astra::CParallelProjectionGeometry2D projGeom;
-	astra::CVolumeGeometry2D volGeom;
-};
-
-BOOST_FIXTURE_TEST_CASE( testParallelBeamLinearKernelProjector2D_General, TestParallelBeamLinearKernelProjector2D )
-{
-
-}
-
-
-// Compute linear kernel for a single volume pixel/detector pixel combination
-float32 compute_linear_kernel(const astra::CProjectionGeometry2D& projgeom, const astra::CVolumeGeometry2D& volgeom,
-                  int iX, int iY, int iDet, float32 fAngle)
-{
-	// projection of center of volume pixel on detector array
-	float32 fDetProj = (iX - (volgeom.getGridColCount()-1.0f)/2.0f ) * volgeom.getPixelLengthX() * cos(fAngle) - (iY - (volgeom.getGridRowCount()-1.0f)/2.0f ) * volgeom.getPixelLengthY() * sin(fAngle);
-
-	// start of detector pixel on detector array
-	float32 fDetStart = projgeom.indexToDetectorOffset(iDet) - 0.5f;
-
-//	printf("(%d,%d,%d): %f in (%f,%f)\n", iX,iY,iDet,fDetProj, fDetStart, fDetStart+1.0f);
-
-	// projection of center of next volume pixel on detector array
-	float32 fDetStep;
-	// length of projection ray through volume pixel
-	float32 fWeight;
-
-	if (fabs(cos(fAngle)) > fabs(sin(fAngle))) {
-		fDetStep = volgeom.getPixelLengthY() * fabs(cos(fAngle));
-		fWeight = projgeom.getDetectorWidth() * volgeom.getPixelLengthX() * 1.0f / fabs(cos(fAngle));
-	} else {
-		fDetStep = volgeom.getPixelLengthX() * fabs(sin(fAngle));
-		fWeight = projgeom.getDetectorWidth() * volgeom.getPixelLengthY() * 1.0f / fabs(sin(fAngle));
-	}
-
-//	printf("step: %f\n   weight: %f\n", fDetStep, fWeight);
-
-	// center of detector pixel on detector array
-	float32 fDetCenter = fDetStart + 0.5f;
-
-	// unweighted contribution of this volume pixel:
-	// linear interpolation between
-	//  fDetCenter - fDetStep    |---> 0
-	//  fDetCenter               |---> 1
-	//  fDetCenter + fDetStep    |---> 0
-	float32 fBase;
-	if (fDetCenter <= fDetProj) {
-		fBase = (fDetCenter - (fDetProj - fDetStep))/fDetStep;
-	} else {
-		fBase = ((fDetProj + fDetStep) - fDetCenter)/fDetStep;
-	}
-//	printf("base: %f\n", fBase);
-	if (fBase < 0) fBase = 0;
-	return fBase * fWeight;
-}
-
-BOOST_AUTO_TEST_CASE( testParallelBeamLinearKernelProjector2D_Rectangles )
-{
-	astra::CParallelBeamLinearKernelProjector2D proj;
-	astra::CParallelProjectionGeometry2D projGeom;
-	astra::CVolumeGeometry2D volGeom;
-
-	const unsigned int iRandomTestCount = 100;
-
-	unsigned int iSeed = time(0);
-	srand(iSeed);
-
-	for (unsigned int iTest = 0; iTest < iRandomTestCount; ++iTest) {
-		int iDetectorCount = 1 + (rand() % 100);
-		int iRows = 1 + (rand() % 100);
-		int iCols = 1 + (rand() % 100);
-		
-		
-		astra::float32 angles[] = { rand() * 2.0f*astra::PI / RAND_MAX };
-		projGeom.initialize(1, iDetectorCount, 0.8f, angles);
-		volGeom.initialize(iCols, iRows);
-		proj.initialize(&projGeom, &volGeom);
-
-		int iMax = proj.getProjectionWeightsCount(0);
-		BOOST_REQUIRE(iMax > 0);
-
-		astra::SPixelWeight* pPix = new astra::SPixelWeight[iMax];
-		BOOST_REQUIRE(pPix);
-
-		astra::float32 fWeight = 0;
-		for (int iDet = 0; iDet < projGeom.getDetectorCount(); ++iDet) {
-			int iCount;
-			proj.computeSingleRayWeights(0, iDet, pPix, iMax, iCount); 
-			BOOST_REQUIRE(iCount <= iMax);
-
-			astra::float32 fW = 0;
-			for (int i = 0; i < iCount; ++i) {
-				float32 fTest = compute_linear_kernel(
-				            projGeom,
-				            volGeom,
-				            pPix[i].m_iIndex % volGeom.getGridColCount(),
-				            pPix[i].m_iIndex / volGeom.getGridColCount(),
-				            iDet,
-				            projGeom.getProjectionAngle(0));
-				BOOST_CHECK_SMALL( pPix[i].m_fWeight - fTest, 0.0004f);
-				fW += pPix[i].m_fWeight;
-			}
-
-			fWeight += fW;
-
-		}
-
-		delete[] pPix;
-	}
-}
-
-