Merge branch 'consistent_scaling'

14 files changed, 132 insertions, 64 deletions

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