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. --- include/astra/ConeProjectionGeometry3D.h | 5 ++++ include/astra/ConeVecProjectionGeometry3D.h | 4 ++++ include/astra/ParallelProjectionGeometry3D.h | 4 ++++ include/astra/ParallelVecProjectionGeometry3D.h | 4 +++- include/astra/ProjectionGeometry3D.h | 15 +++++++++++- src/ConeProjectionGeometry3D.cpp | 32 +++++++++++++++++++++++++ src/ConeVecProjectionGeometry3D.cpp | 9 +++++++ src/CudaForwardProjectionAlgorithm3D.cpp | 19 +++++++++++++++ src/ParallelProjectionGeometry3D.cpp | 16 +++++++++++++ src/ParallelVecProjectionGeometry3D.cpp | 8 +++++++ 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; +} //---------------------------------------------------------------------------------------- -- cgit v1.2.3