summaryrefslogtreecommitdiffstats
path: root/src/GeometryUtil3D.cpp
diff options
context:
space:
mode:
authorWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2015-11-18 11:26:15 +0100
committerWillem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>2015-12-04 12:02:40 +0100
commitb14fb531ad9ae3d565f2cf28f5506408ab10dbed (patch)
tree892d400e72f185b6c4a06aca464343737b135120 /src/GeometryUtil3D.cpp
parent3ea35516aceec4f5817871a00008b109777ebb13 (diff)
downloadastra-b14fb531ad9ae3d565f2cf28f5506408ab10dbed.tar.gz
astra-b14fb531ad9ae3d565f2cf28f5506408ab10dbed.tar.bz2
astra-b14fb531ad9ae3d565f2cf28f5506408ab10dbed.tar.xz
astra-b14fb531ad9ae3d565f2cf28f5506408ab10dbed.zip
Add CompositeGeometryManager
This handles FP and BP operations on multiple data objects at once, splitting them to fit in GPU memory where necessary.
Diffstat (limited to 'src/GeometryUtil3D.cpp')
-rw-r--r--src/GeometryUtil3D.cpp172
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);
+}
+
+
}