Adjust linear/cuda kernels to line integral scaling

4 files changed, 12 insertions, 23 deletions

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/par_fp.cu b/cuda/2d/par_fp.cu
index 0835301..e03381c 100644
--- a/cuda/2d/par_fp.cu
+++ b/cuda/2d/par_fp.cu
@@ -115,10 +115,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 +192,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;
diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.inl b/include/astra/ParallelBeamLinearKernelProjector2D.inl
index 10d4892..1acd422 100644
--- a/include/astra/ParallelBeamLinearKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLinearKernelProjector2D.inl
@@ -154,16 +154,14 @@ 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);
+			lengthPerRow = 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);
+			lengthPerCol = m_pVolumeGeometry->getPixelLengthY() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX);
 			deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
 		}
 
diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py
index 755bd27..5647053 100644
--- a/tests/python/test_line2d.py
+++ b/tests/python/test_line2d.py
@@ -486,17 +486,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]
@@ -504,9 +496,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]
@@ -514,7 +506,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)):