Unify some parallel_vec parameter computations

4 files changed, 28 insertions, 34 deletions

diff --git a/cuda/2d/par_fp.cu b/cuda/2d/par_fp.cu
index 10ce683..a0b04ee 100644
--- a/cuda/2d/par_fp.cu
+++ b/cuda/2d/par_fp.cu
@@ -251,36 +251,14 @@ __global__ void FPvertical_simple(float* D_projData, unsigned int projPitch, uns
 static void convertAndUploadAngles(const SParProjection *projs, unsigned int nth, unsigned int ndets)
 {
 	float *angles = new float[nth];
-	float *offset = new float[nth];
+	float *offsets = new float[nth];
 	float *detsizes = new float[nth];
 
-	for (int i = 0; i < nth; ++i) {
-
-#warning TODO: Replace by getParParamaters
-
-		// Take part of DetU orthogonal to Ray
-		double ux = projs[i].fDetUX;
-		double uy = projs[i].fDetUY;
-
-		double t = (ux * projs[i].fRayX + uy * projs[i].fRayY) / (projs[i].fRayX * projs[i].fRayX + projs[i].fRayY * projs[i].fRayY);
-
-		ux -= t * projs[i].fRayX;
-		uy -= t * projs[i].fRayY;
-
-		double angle = atan2(uy, ux);
-
-		angles[i] = angle;
-
-		double norm2 = uy * uy + ux * ux;
-
-		detsizes[i] = sqrt(norm2);
-
-		// CHECKME: SIGNS?
-		offset[i] = -0.5*ndets - (projs[i].fDetSY*uy + projs[i].fDetSX*ux) / norm2;
-	}
+	for (int i = 0; i < nth; ++i)
+		getParParameters(projs[i], ndets, angles[i], detsizes[i], offsets[i]);
 
 	cudaMemcpyToSymbol(gC_angle, angles, nth*sizeof(float), 0, cudaMemcpyHostToDevice); 
-	cudaMemcpyToSymbol(gC_angle_offset, offset, nth*sizeof(float), 0, cudaMemcpyHostToDevice); 
+	cudaMemcpyToSymbol(gC_angle_offset, offsets, nth*sizeof(float), 0, cudaMemcpyHostToDevice);
 	cudaMemcpyToSymbol(gC_angle_detsize, detsizes, nth*sizeof(float), 0, cudaMemcpyHostToDevice); 
 }
 
diff --git a/cuda/2d/sirt.cu b/cuda/2d/sirt.cu
index 2516c6c..b393d7f 100644
--- a/cuda/2d/sirt.cu
+++ b/cuda/2d/sirt.cu
@@ -184,9 +184,8 @@ bool SIRT::doSlabCorrections()
 	float bound = cosf(1.3963f);
 	float* t = (float*)D_sinoData;
 	for (int i = 0; i < dims.iProjAngles; ++i) {
-		// TODO: Checkme
-		// TODO: Replace by getParParameters
-		double angle = atan2(parProjs[i].fRayX, -parProjs[i].fRayY);
+		float angle, detsize, offset;
+		getParParameters(parProjs[i], dims.iProjDets, angle, detsize, offset);
 		float f = fabs(cosf(angle));
 
 		if (f < bound)
diff --git a/include/astra/GeometryUtil2D.h b/include/astra/GeometryUtil2D.h
index 4d79353..914e40d 100644
--- a/include/astra/GeometryUtil2D.h
+++ b/include/astra/GeometryUtil2D.h
@@ -96,8 +96,11 @@ SFanProjection* genFanProjections(unsigned int iProjAngles,
                                   double fDetSize,
                                   const float *pfAngles);
 
+bool getParParameters(const SParProjection &proj, unsigned int iProjDets, float &fAngle, float &fDetSize, float &fOffset);
+
 bool getFanParameters(const SFanProjection &proj, unsigned int iProjDets, float &fAngle, float &fOriginSource, float &fOriginDetector, float &fDetSize, float &fOffset);
 
+
 }
 
 #endif
diff --git a/src/GeometryUtil2D.cpp b/src/GeometryUtil2D.cpp
index 1bca2bc..2ee6185 100644
--- a/src/GeometryUtil2D.cpp
+++ b/src/GeometryUtil2D.cpp
@@ -104,13 +104,27 @@ SFanProjection* genFanProjections(unsigned int iProjAngles,
 
 // Convert a SParProjection back into its set of "standard" circular parallel
 // beam parameters. This is always possible.
-bool getParParameters(const SParProjection &proj /* , ... */)
+bool getParParameters(const SParProjection &proj, unsigned int iProjDets, float &fAngle, float &fDetSize, float &fOffset)
 {
-	// angle
-	// det size
-	// offset
+	// Take part of DetU orthogonal to Ray
+	double ux = proj.fDetUX;
+	double uy = proj.fDetUY;
+
+	double t = (ux * proj.fRayX + uy * proj.fRayY) / (proj.fRayX * proj.fRayX + proj.fRayY * proj.fRayY);
+
+	ux -= t * proj.fRayX;
+	uy -= t * proj.fRayY;
+
+	double angle = atan2(uy, ux);
+
+	fAngle = (float)angle;
+
+	double norm2 = uy * uy + ux * ux;
+
+	fDetSize = (float)sqrt(norm2);
 
-	// (see convertAndUploadAngles in par_fp.cu)
+	// CHECKME: SIGNS?
+	fOffset = (float)(-0.5*iProjDets - (proj.fDetSY*uy + proj.fDetSX*ux) / norm2);
 
 	return true;
 }