From b3e8338a7fa4c7ed9a5954ca02fa3126aefff530 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn 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