Add ProjectionGeometry3D::projectPoint for par3d and cone.

Par3d_vec and cone_vec to follow.

10 files changed, 114 insertions, 2 deletions

diff --git a/include/astra/ConeProjectionGeometry3D.h b/include/astra/ConeProjectionGeometry3D.h
index aba8797..ce17353 100644
--- a/include/astra/ConeProjectionGeometry3D.h
+++ b/include/astra/ConeProjectionGeometry3D.h
@@ -185,6 +185,11 @@ public:
 	 * Returns a vector giving the projection direction for a projection and detector index
 	 */
 	virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex) const;
+
+	virtual void projectPoint(float32 fX, float32 fY, float32 fZ,
+	                          int iAngleIndex,
+	                          float32 &fU, float32 &fV) const;
+
 };
 
 // Returns the distance from the origin of the coordinate system to the source.
diff --git a/include/astra/ConeVecProjectionGeometry3D.h b/include/astra/ConeVecProjectionGeometry3D.h
index 1765cdd..e58200a 100644
--- a/include/astra/ConeVecProjectionGeometry3D.h
+++ b/include/astra/ConeVecProjectionGeometry3D.h
@@ -147,6 +147,10 @@ public:
 	 virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex) const;
 
 	const SConeProjection* getProjectionVectors() const { return m_pProjectionAngles; }
+
+	virtual void projectPoint(float32 fX, float32 fY, float32 fZ,
+	                          int iAngleIndex,
+	                          float32 &fU, float32 &fV) const;
 };
 
 } // namespace astra
diff --git a/include/astra/ParallelProjectionGeometry3D.h b/include/astra/ParallelProjectionGeometry3D.h
index a981104..907eead 100644
--- a/include/astra/ParallelProjectionGeometry3D.h
+++ b/include/astra/ParallelProjectionGeometry3D.h
@@ -147,6 +147,10 @@ public:
 	  */
 	virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex) const;
 
+	virtual void projectPoint(float32 fX, float32 fY, float32 fZ,
+	                          int iAngleIndex,
+	                          float32 &fU, float32 &fV) const;
+
 	 /**
 	  * Creates (= allocates) a 2D projection geometry used when projecting one slice using a 2D projector
 	  *
diff --git a/include/astra/ParallelVecProjectionGeometry3D.h b/include/astra/ParallelVecProjectionGeometry3D.h
index 0b4a766..0288f3e 100644
--- a/include/astra/ParallelVecProjectionGeometry3D.h
+++ b/include/astra/ParallelVecProjectionGeometry3D.h
@@ -149,7 +149,9 @@ public:
 
 	const SPar3DProjection* getProjectionVectors() const { return m_pProjectionAngles; }
 
-
+	virtual void projectPoint(float32 fX, float32 fY, float32 fZ,
+	                          int iAngleIndex,
+	                          float32 &fU, float32 &fV) const;
 };
 
 } // namespace astra
diff --git a/include/astra/ProjectionGeometry3D.h b/include/astra/ProjectionGeometry3D.h
index 20ad8f3..9e7f787 100644
--- a/include/astra/ProjectionGeometry3D.h
+++ b/include/astra/ProjectionGeometry3D.h
@@ -298,10 +298,23 @@ public:
 	 *
 	 * @param _iIndex	the index of the detector.
 	 * @param _iAngleIndex	output: index of angle
-	 * @param _iDetectorIndex output: index of dectecor
+	 * @param _iDetectorIndex output: index of detector
 	 */
 	virtual void indexToAngleDetectorIndex(int _iIndex, int& _iAngleIndex, int& _iDetectorIndex) const;
 
+	/** Project a point onto the detector. The 3D point coordinates
+	 * are in units. The output fU,fV are the (unrounded) indices of the
+	 * detector column and row.
+	 * This may fall outside of the actual detector.
+	 *
+	 * @param fX,fY,fZ	coordinates of the point to project
+	 * @param iAngleIndex	the index of the angle to use
+	 * @param fU,fV		the projected point.
+	 */
+	virtual void projectPoint(float32 fX, float32 fY, float32 fZ,
+	                          int iAngleIndex,
+	                          float32 &fU, float32 &fV) const = 0;
+
 	/** Returns true if the type of geometry defined in this class is the one specified in _sType.
 	 *
 	 * @param _sType geometry type to compare to.
diff --git a/src/ConeProjectionGeometry3D.cpp b/src/ConeProjectionGeometry3D.cpp
index 129e675..4cdb9e4 100644
--- a/src/ConeProjectionGeometry3D.cpp
+++ b/src/ConeProjectionGeometry3D.cpp
@@ -225,4 +225,36 @@ CVector3D CConeProjectionGeometry3D::getProjectionDirection(int _iProjectionInde
 	return ret;
 }
 
+void CConeProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ,
+                                                 int iAngleIndex,
+                                                 float32 &fU, float32 &fV) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	float alpha = m_pfProjectionAngles[iAngleIndex];
+
+	// Project point onto optical axis
+
+	// Projector direction is (cos(alpha), sin(alpha))
+	// Vector source->origin is (-sin(alpha), cos(alpha))
+
+	// Distance from source, projected on optical axis
+	float fD = -sin(alpha) * fX + cos(alpha) * fY + m_fOriginSourceDistance;
+
+	// Scale fZ to detector plane
+	fV = detectorOffsetYToRowIndexFloat( (fZ * (m_fOriginSourceDistance + m_fOriginDetectorDistance)) / fD );
+
+
+	// Orthogonal distance in XY-plane to optical axis
+	float fS = cos(alpha) * fX + sin(alpha) * fY;
+
+	// Scale fS to detector plane
+	fU = detectorOffsetXToColIndexFloat( (fS * (m_fOriginSourceDistance + m_fOriginDetectorDistance)) / fD );
+
+	fprintf(stderr, "alpha: %f, D: %f, V: %f, S: %f, U: %f\n", alpha, fD, fV, fS, fU);
+
+}
+
+
 } // end namespace astra
diff --git a/src/ConeVecProjectionGeometry3D.cpp b/src/ConeVecProjectionGeometry3D.cpp
index 875a2c7..29d84b7 100644
--- a/src/ConeVecProjectionGeometry3D.cpp
+++ b/src/ConeVecProjectionGeometry3D.cpp
@@ -220,6 +220,15 @@ CVector3D CConeVecProjectionGeometry3D::getProjectionDirection(int _iProjectionI
 	return CVector3D(p.fDetSX + (u+0.5)*p.fDetUX + (v+0.5)*p.fDetVX - p.fSrcX, p.fDetSY + (u+0.5)*p.fDetUY + (v+0.5)*p.fDetVY - p.fSrcY, p.fDetSZ + (u+0.5)*p.fDetUZ + (v+0.5)*p.fDetVZ - p.fSrcZ);
 }
 
+void CConeVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ,
+                                                 int iAngleIndex,
+                                                 float32 &fU, float32 &fV) const
+{
+#warning implementme
+	fU = 0.0f/0.0f;
+	fV = 0.0f/0.0f;
+}
+
 
 //----------------------------------------------------------------------------------------
 
diff --git a/src/CudaForwardProjectionAlgorithm3D.cpp b/src/CudaForwardProjectionAlgorithm3D.cpp
index f64620f..76b3419 100644
--- a/src/CudaForwardProjectionAlgorithm3D.cpp
+++ b/src/CudaForwardProjectionAlgorithm3D.cpp
@@ -251,6 +251,25 @@ void CCudaForwardProjectionAlgorithm3D::run(int)
 		projKernel = projector->getProjectionKernel();
 	}
 
+#if 0
+	// Debugging code that gives the coordinates of the corners of the volume
+	// projected on the detector.
+	{
+		float fX[] = { volgeom.getWindowMinX(), volgeom.getWindowMaxX() };
+		float fY[] = { volgeom.getWindowMinY(), volgeom.getWindowMaxY() };
+		float fZ[] = { volgeom.getWindowMinZ(), volgeom.getWindowMaxZ() };
+
+		for (int a = 0; a < projgeom->getProjectionCount(); ++a)
+		for (int i = 0; i < 2; ++i)
+		for (int j = 0; j < 2; ++j)
+		for (int k = 0; k < 2; ++k) {
+			float fU, fV;
+			projgeom->projectPoint(fX[i], fY[j], fZ[k], a, fU, fV);
+			fprintf(stderr, "%3d %c1,%c1,%c1 -> %12f %12f\n", a, i ? ' ' : '-', j ? ' ' : '-', k ? ' ' : '-', fU, fV);
+		}
+	}
+#endif
+
 	if (conegeom) {
 		astraCudaConeFP(m_pVolume->getDataConst(), m_pProjections->getData(),
 		                volgeom.getGridColCount(),
diff --git a/src/ParallelProjectionGeometry3D.cpp b/src/ParallelProjectionGeometry3D.cpp
index 2027f95..b5cdfb6 100644
--- a/src/ParallelProjectionGeometry3D.cpp
+++ b/src/ParallelProjectionGeometry3D.cpp
@@ -179,6 +179,22 @@ CVector3D CParallelProjectionGeometry3D::getProjectionDirection(int _iProjection
 	return CVector3D(fDirX, fDirY, fDirZ);
 }
 
+void CParallelProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ,
+                                                 int iAngleIndex,
+                                                 float32 &fU, float32 &fV) const
+{
+	ASTRA_ASSERT(iAngleIndex >= 0);
+	ASTRA_ASSERT(iAngleIndex < m_iProjectionAngleCount);
+
+	// V (detector row)
+	fV = detectorOffsetYToRowIndexFloat(fZ);
+
+	// U (detector column)
+	float alpha = m_pfProjectionAngles[iAngleIndex];
+	// projector direction is (cos(alpha), sin(alpha))
+	fU = detectorOffsetXToColIndexFloat(cos(alpha) * fX + sin(alpha) * fY);
+}
+
 CParallelProjectionGeometry2D * CParallelProjectionGeometry3D::createProjectionGeometry2D() const
 {
 	const float32 * pfProjectionAngles = getProjectionAngles(); //new float32[getProjectionCount()];
diff --git a/src/ParallelVecProjectionGeometry3D.cpp b/src/ParallelVecProjectionGeometry3D.cpp
index c1265dd..7c2d2cd 100644
--- a/src/ParallelVecProjectionGeometry3D.cpp
+++ b/src/ParallelVecProjectionGeometry3D.cpp
@@ -218,6 +218,14 @@ CVector3D CParallelVecProjectionGeometry3D::getProjectionDirection(int _iProject
 	return CVector3D(p.fRayX, p.fRayY, p.fRayZ);
 }
 
+void CParallelVecProjectionGeometry3D::projectPoint(float32 fX, float32 fY, float32 fZ,
+                                                 int iAngleIndex,
+                                                 float32 &fU, float32 &fV) const
+{
+#warning implementme
+	fU = 0.0f/0.0f;
+	fV = 0.0f/0.0f;
+}
 
 //----------------------------------------------------------------------------------------