summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
-rw-r--r--cuda/3d/astra3d.cu150
-rw-r--r--cuda/3d/astra3d.h16
-rw-r--r--include/astra/GeometryUtil3D.h53
3 files changed, 219 insertions, 0 deletions
diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu
index 0b9c70b..f672d6c 100644
--- a/cuda/3d/astra3d.cu
+++ b/cuda/3d/astra3d.cu
@@ -40,6 +40,12 @@ $Id$
#include "arith3d.h"
#include "astra3d.h"
+#include "astra/ParallelProjectionGeometry3D.h"
+#include "astra/ParallelVecProjectionGeometry3D.h"
+#include "astra/ConeProjectionGeometry3D.h"
+#include "astra/ConeVecProjectionGeometry3D.h"
+#include "astra/VolumeGeometry3D.h"
+
#include <iostream>
using namespace astraCUDA3d;
@@ -137,6 +143,150 @@ static SPar3DProjection* genPar3DProjections(unsigned int iProjAngles,
+
+// adjust pProjs to normalize volume geometry
+template<typename ProjectionT>
+static bool convertAstraGeometry_internal(const CVolumeGeometry3D* pVolGeom,
+ unsigned int iProjectionAngleCount,
+ ProjectionT*& pProjs,
+ float& fOutputScale)
+{
+ assert(pVolGeom);
+ assert(pProjs);
+
+ // TODO: Relative instead of absolute
+ const float EPS = 0.00001f;
+ if (abs(pVolGeom->getPixelLengthX() - pVolGeom->getPixelLengthY()) > EPS)
+ return false;
+ if (abs(pVolGeom->getPixelLengthX() - pVolGeom->getPixelLengthZ()) > EPS)
+ return false;
+
+
+ // Translate
+ float dx = -(pVolGeom->getWindowMinX() + pVolGeom->getWindowMaxX()) / 2;
+ float dy = -(pVolGeom->getWindowMinY() + pVolGeom->getWindowMaxY()) / 2;
+ float dz = -(pVolGeom->getWindowMinZ() + pVolGeom->getWindowMaxZ()) / 2;
+
+ float factor = 1.0f / pVolGeom->getPixelLengthX();
+
+ for (int i = 0; i < iProjectionAngleCount; ++i) {
+ // CHECKME: Order of scaling and translation
+ pProjs[i].translate(dx, dy, dz);
+ pProjs[i].scale(factor);
+ }
+
+ // CHECKME: Check factor
+ fOutputScale *= pVolGeom->getPixelLengthX();
+
+ return true;
+}
+
+
+
+bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
+ const CParallelProjectionGeometry3D* pProjGeom,
+ SPar3DProjection*& pProjs, float& fOutputScale)
+{
+ assert(pVolGeom);
+ assert(pProjGeom);
+ assert(pProjGeom->getProjectionAngles());
+
+ int nth = pProjGeom->getProjectionCount();
+
+ pProjs = genPar3DProjections(nth,
+ pProjGeom->getDetectorColCount(),
+ pProjGeom->getDetectorRowCount(),
+ pProjGeom->getDetectorSpacingX(),
+ pProjGeom->getDetectorSpacingY(),
+ pProjGeom->getProjectionAngles());
+
+ bool ok;
+
+ fOutputScale = 1.0f;
+
+ ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale);
+
+ return ok;
+}
+
+bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
+ const CParallelVecProjectionGeometry3D* pProjGeom,
+ SPar3DProjection*& pProjs, float& fOutputScale)
+{
+ assert(pVolGeom);
+ assert(pProjGeom);
+ assert(pProjGeom->getProjectionVectors());
+
+ int nth = pProjGeom->getProjectionCount();
+
+ pProjs = new SPar3DProjection[nth];
+ for (int i = 0; i < nth; ++i)
+ pProjs[i] = pProjGeom->getProjectionVectors()[i];
+
+ bool ok;
+
+ fOutputScale = 1.0f;
+
+ ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale);
+
+ return ok;
+}
+
+bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
+ const CConeProjectionGeometry3D* pProjGeom,
+ SConeProjection*& pProjs, float& fOutputScale)
+{
+ assert(pVolGeom);
+ assert(pProjGeom);
+ assert(pProjGeom->getProjectionAngles());
+
+ int nth = pProjGeom->getProjectionCount();
+
+ pProjs = genConeProjections(nth,
+ pProjGeom->getDetectorColCount(),
+ pProjGeom->getDetectorRowCount(),
+ pProjGeom->getOriginSourceDistance(),
+ pProjGeom->getOriginDetectorDistance(),
+ pProjGeom->getDetectorSpacingX(),
+ pProjGeom->getDetectorSpacingY(),
+ pProjGeom->getProjectionAngles());
+
+ bool ok;
+
+ fOutputScale = 1.0f;
+
+ ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale);
+
+ return ok;
+}
+
+bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
+ const CConeVecProjectionGeometry3D* pProjGeom,
+ SConeProjection*& pProjs, float& fOutputScale)
+{
+ assert(pVolGeom);
+ assert(pProjGeom);
+ assert(pProjGeom->getProjectionAngles());
+
+ int nth = pProjGeom->getProjectionCount();
+
+ pProjs = new SConeProjection[nth];
+ for (int i = 0; i < nth; ++i)
+ pProjs[i] = pProjGeom->getProjectionVectors()[i];
+
+ bool ok;
+
+ fOutputScale = 1.0f;
+
+ ok = convertAstraGeometry_internal(pVolGeom, nth, pProjs, fOutputScale);
+
+ return ok;
+}
+
+
+
+
+
class AstraSIRT3d_internal {
public:
SDimensions3D dims;
diff --git a/cuda/3d/astra3d.h b/cuda/3d/astra3d.h
index f91fe26..47e252e 100644
--- a/cuda/3d/astra3d.h
+++ b/cuda/3d/astra3d.h
@@ -466,6 +466,22 @@ _AstraExport bool astraCudaFDK(float* pfVolume, const float* pfProjections,
bool bShortScan,
int iGPUIndex, int iVoxelSuperSampling);
+_AstraExport bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
+ const CParallelProjectionGeometry3D* pProjGeom,
+ SPar3DProjection*& pProjs, float& fOutputScale);
+
+_AstraExport bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
+ const CParallelVecProjectionGeometry3D* pProjGeom,
+ SPar3DProjection*& pProjs, float& fOutputScale);
+
+_AstraExport bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
+ const CConeProjectionGeometry3D* pProjGeom,
+ SConeProjection*& pProjs, float& fOutputScale);
+
+_AstraExport bool convertAstraGeometry(const CVolumeGeometry3D* pVolGeom,
+ const CConeVecProjectionGeometry3D* pProjGeom,
+ SConeProjection*& pProjs, float& fOutputScale);
+
}
diff --git a/include/astra/GeometryUtil3D.h b/include/astra/GeometryUtil3D.h
index 698372e..6ceac63 100644
--- a/include/astra/GeometryUtil3D.h
+++ b/include/astra/GeometryUtil3D.h
@@ -43,6 +43,33 @@ struct SConeProjection {
// the V-edge of a detector pixel
double fDetVX, fDetVY, fDetVZ;
+
+
+
+
+ void translate(double dx, double dy, double dz) {
+ fSrcX += dx;
+ fSrcY += dy;
+ fSrcZ += dz;
+ fDetSX += dx;
+ fDetSY += dy;
+ fDetSZ += dz;
+
+ }
+ void scale(double factor) {
+ fSrcX *= factor;
+ fSrcY *= factor;
+ fSrcZ *= factor;
+ fDetSX *= factor;
+ fDetSY *= factor;
+ fDetSZ *= factor;
+ fDetUX *= factor;
+ fDetUY *= factor;
+ fDetUZ *= factor;
+ fDetVX *= factor;
+ fDetVY *= factor;
+ fDetVZ *= factor;
+ }
};
struct SPar3DProjection {
@@ -57,6 +84,29 @@ struct SPar3DProjection {
// the V-edge of a detector pixel
double fDetVX, fDetVY, fDetVZ;
+
+
+
+
+ void translate(double dx, double dy, double dz) {
+ fDetSX += dx;
+ fDetSY += dy;
+ fDetSZ += dz;
+ }
+ void scale(double factor) {
+ fRayX *= factor;
+ fRayY *= factor;
+ fRayZ *= factor;
+ fDetSX *= factor;
+ fDetSY *= factor;
+ fDetSZ *= factor;
+ fDetUX *= factor;
+ fDetUY *= factor;
+ fDetUZ *= factor;
+ fDetVX *= factor;
+ fDetVY *= factor;
+ fDetVZ *= factor;
+ }
};
void computeBP_UV_Coeffs(const SPar3DProjection& proj,
@@ -68,6 +118,9 @@ void computeBP_UV_Coeffs(const SConeProjection& proj,
double &fVX, double &fVY, double &fVZ, double &fVC,
double &fDX, double &fDY, double &fDZ, double &fDC);
+
+
+
}
#endif