From b3e8338a7fa4c7ed9a5954ca02fa3126aefff530 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Tue, 2 Dec 2014 14:20:46 +0100
Subject: Add ProjectionGeometry3D::projectPoint for par3d and cone.
Par3d_vec and cone_vec to follow.
---
src/ConeProjectionGeometry3D.cpp | 32 ++++++++++++++++++++++++++++++++
src/ConeVecProjectionGeometry3D.cpp | 9 +++++++++
src/CudaForwardProjectionAlgorithm3D.cpp | 19 +++++++++++++++++++
src/ParallelProjectionGeometry3D.cpp | 16 ++++++++++++++++
src/ParallelVecProjectionGeometry3D.cpp | 8 ++++++++
5 files changed, 84 insertions(+)
(limited to 'src')
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;
+}
//----------------------------------------------------------------------------------------
--
cgit v1.2.3