diff options
Diffstat (limited to 'src/GeometryUtil3D.cpp')
| -rw-r--r-- | src/GeometryUtil3D.cpp | 172 |
1 files changed, 172 insertions, 0 deletions
diff --git a/src/GeometryUtil3D.cpp b/src/GeometryUtil3D.cpp index 52dd5a9..c6bfd8b 100644 --- a/src/GeometryUtil3D.cpp +++ b/src/GeometryUtil3D.cpp @@ -28,8 +28,96 @@ $Id$ #include "astra/GeometryUtil3D.h" +#include <cmath> + namespace astra { + +SConeProjection* genConeProjections(unsigned int iProjAngles, + unsigned int iProjU, + unsigned int iProjV, + double fOriginSourceDistance, + double fOriginDetectorDistance, + double fDetUSize, + double fDetVSize, + const float *pfAngles) +{ + SConeProjection base; + base.fSrcX = 0.0f; + base.fSrcY = -fOriginSourceDistance; + base.fSrcZ = 0.0f; + + base.fDetSX = iProjU * fDetUSize * -0.5f; + base.fDetSY = fOriginDetectorDistance; + base.fDetSZ = iProjV * fDetVSize * -0.5f; + + base.fDetUX = fDetUSize; + base.fDetUY = 0.0f; + base.fDetUZ = 0.0f; + + base.fDetVX = 0.0f; + base.fDetVY = 0.0f; + base.fDetVZ = fDetVSize; + + SConeProjection* p = new SConeProjection[iProjAngles]; + +#define ROTATE0(name,i,alpha) do { p[i].f##name##X = base.f##name##X * cos(alpha) - base.f##name##Y * sin(alpha); p[i].f##name##Y = base.f##name##X * sin(alpha) + base.f##name##Y * cos(alpha); p[i].f##name##Z = base.f##name##Z; } while(0) + + for (unsigned int i = 0; i < iProjAngles; ++i) { + ROTATE0(Src, i, pfAngles[i]); + ROTATE0(DetS, i, pfAngles[i]); + ROTATE0(DetU, i, pfAngles[i]); + ROTATE0(DetV, i, pfAngles[i]); + } + +#undef ROTATE0 + + return p; +} + +SPar3DProjection* genPar3DProjections(unsigned int iProjAngles, + unsigned int iProjU, + unsigned int iProjV, + double fDetUSize, + double fDetVSize, + const float *pfAngles) +{ + SPar3DProjection base; + base.fRayX = 0.0f; + base.fRayY = 1.0f; + base.fRayZ = 0.0f; + + base.fDetSX = iProjU * fDetUSize * -0.5f; + base.fDetSY = 0.0f; + base.fDetSZ = iProjV * fDetVSize * -0.5f; + + base.fDetUX = fDetUSize; + base.fDetUY = 0.0f; + base.fDetUZ = 0.0f; + + base.fDetVX = 0.0f; + base.fDetVY = 0.0f; + base.fDetVZ = fDetVSize; + + SPar3DProjection* p = new SPar3DProjection[iProjAngles]; + +#define ROTATE0(name,i,alpha) do { p[i].f##name##X = base.f##name##X * cos(alpha) - base.f##name##Y * sin(alpha); p[i].f##name##Y = base.f##name##X * sin(alpha) + base.f##name##Y * cos(alpha); p[i].f##name##Z = base.f##name##Z; } while(0) + + for (unsigned int i = 0; i < iProjAngles; ++i) { + ROTATE0(Ray, i, pfAngles[i]); + ROTATE0(DetS, i, pfAngles[i]); + ROTATE0(DetU, i, pfAngles[i]); + ROTATE0(DetV, i, pfAngles[i]); + } + +#undef ROTATE0 + + return p; +} + + + + // (See declaration in header for (mathematical) description of these functions) @@ -72,4 +160,88 @@ void computeBP_UV_Coeffs(const SConeProjection& proj, double &fUX, double &fUY, } +// TODO: Handle cases of rays parallel to coordinate planes + +void backprojectPointX(const SPar3DProjection& proj, double fU, double fV, + double fX, double &fY, double &fZ) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - px) / proj.fRayX; + + fY = py + a * proj.fRayY; + fZ = pz + a * proj.fRayZ; +} + +void backprojectPointY(const SPar3DProjection& proj, double fU, double fV, + double fY, double &fX, double &fZ) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - py) / proj.fRayY; + + fX = px + a * proj.fRayX; + fZ = pz + a * proj.fRayZ; + +} + +void backprojectPointZ(const SPar3DProjection& proj, double fU, double fV, + double fZ, double &fX, double &fY) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - pz) / proj.fRayZ; + + fX = px + a * proj.fRayX; + fY = py + a * proj.fRayY; +} + + + +void backprojectPointX(const SConeProjection& proj, double fU, double fV, + double fX, double &fY, double &fZ) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fX - proj.fSrcX) / (px - proj.fSrcX); + + fY = proj.fSrcY + a * (py - proj.fSrcY); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); +} + +void backprojectPointY(const SConeProjection& proj, double fU, double fV, + double fY, double &fX, double &fZ) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fY - proj.fSrcY) / (py - proj.fSrcY); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fZ = proj.fSrcZ + a * (pz - proj.fSrcZ); +} + +void backprojectPointZ(const SConeProjection& proj, double fU, double fV, + double fZ, double &fX, double &fY) +{ + double px = proj.fDetSX + fU * proj.fDetUX + fV * proj.fDetVX; + double py = proj.fDetSY + fU * proj.fDetUY + fV * proj.fDetVY; + double pz = proj.fDetSZ + fU * proj.fDetUZ + fV * proj.fDetVZ; + + double a = (fZ - proj.fSrcZ) / (pz - proj.fSrcZ); + + fX = proj.fSrcX + a * (px - proj.fSrcX); + fY = proj.fSrcY + a * (py - proj.fSrcY); +} + + } |