added get_geometry for 3d projection data objects

3 files changed, 123 insertions, 10 deletions

diff --git a/matlab/mex/astra_mex_data3d_c.cpp b/matlab/mex/astra_mex_data3d_c.cpp
index 47316f5..84bc0a4 100644
--- a/matlab/mex/astra_mex_data3d_c.cpp
+++ b/matlab/mex/astra_mex_data3d_c.cpp
@@ -295,9 +295,8 @@ void astra_mex_data3d_get_geometry(int nlhs, mxArray* plhs[], int nrhs, const mx
 	// create output
 	if (1 <= nlhs) {
 		if (pDataObject->getType() == CFloat32Data3D::PROJECTION) {
-			// CFloat32ProjectionData2D* pDataObject2 = dynamic_cast<CFloat32ProjectionData2D*>(pDataObject);
-			// plhs[0] = createProjectionGeometryStruct(pDataObject2->getGeometry());
-			mexErrMsgTxt("Not implemented yet. \n");
+			CFloat32ProjectionData3DMemory* pDataObject2 = dynamic_cast<CFloat32ProjectionData3DMemory*>(pDataObject);
+			plhs[0] = createProjectionGeometryStruct(pDataObject2->getGeometry());
 		}
 		else if (pDataObject->getType() == CFloat32Data3D::VOLUME) {
 			CFloat32VolumeData3DMemory* pDataObject2 = dynamic_cast<CFloat32VolumeData3DMemory*>(pDataObject);
diff --git a/matlab/mex/mexHelpFunctions.cpp b/matlab/mex/mexHelpFunctions.cpp
index 63d2003..654f5c6 100644
--- a/matlab/mex/mexHelpFunctions.cpp
+++ b/matlab/mex/mexHelpFunctions.cpp
@@ -207,7 +207,7 @@ astra::CProjectionGeometry2D* parseProjectionGeometryStruct(const mxArray* prhs)
 }
 
 //-----------------------------------------------------------------------------------------
-// create projection geometry data
+// create 2D projection geometry struct
 mxArray* createProjectionGeometryStruct(astra::CProjectionGeometry2D* _pProjGeom)
 {
 	// temporary map to store the data for the MATLAB struct
@@ -224,7 +224,7 @@ mxArray* createProjectionGeometryStruct(astra::CProjectionGeometry2D* _pProjGeom
 		for (int i = 0; i < _pProjGeom->getProjectionAngleCount(); i++) {
 			out[i] = _pProjGeom->getProjectionAngle(i);
 		}
-		mGeometryInfo["ProjectionAngles"] = pAngles;	
+		mGeometryInfo["ProjectionAngles"] = pAngles;
 	}
 
 	// fanflat
@@ -242,7 +242,7 @@ mxArray* createProjectionGeometryStruct(astra::CProjectionGeometry2D* _pProjGeom
 		for (int i = 0; i < pFanFlatGeom->getProjectionAngleCount(); i++) {
 			out[i] = pFanFlatGeom->getProjectionAngle(i);
 		}
-		mGeometryInfo["ProjectionAngles"] = pAngles;	
+		mGeometryInfo["ProjectionAngles"] = pAngles;
 	}
 
 	// fanflat_vec
@@ -263,7 +263,7 @@ mxArray* createProjectionGeometryStruct(astra::CProjectionGeometry2D* _pProjGeom
 			out[2*iAngleCount + i] = p->fDetSX + 0.5f*iDetCount*p->fDetUX;
 			out[3*iAngleCount + i] = p->fDetSY + 0.5f*iDetCount*p->fDetUY;
 			out[4*iAngleCount + i] = p->fDetUX;
-			out[5*iAngleCount + i] = p->fDetUY;			
+			out[5*iAngleCount + i] = p->fDetUY;
 		}
 		mGeometryInfo["Vectors"] = pVectors;
 	}
@@ -280,6 +280,115 @@ mxArray* createProjectionGeometryStruct(astra::CProjectionGeometry2D* _pProjGeom
 }
 
 //-----------------------------------------------------------------------------------------
+// create 3D projection geometry struct
+mxArray* createProjectionGeometryStruct(astra::CProjectionGeometry3D* _pProjGeom)
+{
+	// temporary map to store the data for the MATLAB struct
+	std::map<std::string, mxArray*> mGeometryInfo;
+
+	// parallel beam
+	if (_pProjGeom->isOfType("parallel3d")) {
+		mGeometryInfo["type"] = mxCreateString("parallel3d");
+		mGeometryInfo["DetectorRowCount"] = mxCreateDoubleScalar(_pProjGeom->getDetectorRowCount());
+		mGeometryInfo["DetectorColCount"] = mxCreateDoubleScalar(_pProjGeom->getDetectorColCount());
+		mGeometryInfo["DetectorSpacingX"] = mxCreateDoubleScalar(_pProjGeom->getDetectorSpacingX());
+		mGeometryInfo["DetectorSpacingY"] = mxCreateDoubleScalar(_pProjGeom->getDetectorSpacingY());
+
+		mxArray* pAngles = mxCreateDoubleMatrix(1, _pProjGeom->getProjectionCount(), mxREAL);
+		double* out = mxGetPr(pAngles);
+		for (int i = 0; i < _pProjGeom->getProjectionCount(); i++) {
+			out[i] = _pProjGeom->getProjectionAngle(i);
+		}
+		mGeometryInfo["ProjectionAngles"] = pAngles;
+	}
+
+	// parallel beam vector
+	if (_pProjGeom->isOfType("parallel3d_vec")) {
+		astra::CParallelVecProjectionGeometry3D* pVecGeom = dynamic_cast<astra::CParallelVecProjectionGeometry3D*>(_pProjGeom);
+
+		mGeometryInfo["type"] = mxCreateString("parallel3d_vec");
+		mGeometryInfo["DetectorRowCount"] = mxCreateDoubleScalar(pVecGeom->getDetectorRowCount());
+		mGeometryInfo["DetectorColCount"] = mxCreateDoubleScalar(pVecGeom->getDetectorColCount());
+
+		mxArray* pVectors = mxCreateDoubleMatrix(pVecGeom->getProjectionCount(), 12, mxREAL);
+		double* out = mxGetPr(pVectors);
+		int iDetRowCount = pVecGeom->getDetectorRowCount();
+		int iDetColCount = pVecGeom->getDetectorColCount();
+		int iAngleCount = pVecGeom->getProjectionCount();
+		for (int i = 0; i < pVecGeom->getProjectionCount(); i++) {
+			const SPar3DProjection* p = &pVecGeom->getProjectionVectors()[i];
+			out[ 0*iAngleCount + i] = p->fRayX;
+			out[ 1*iAngleCount + i] = p->fRayY;
+			out[ 2*iAngleCount + i] = p->fRayZ;
+			out[ 3*iAngleCount + i] = p->fDetSX + 0.5f*iDetRowCount*p->fDetUX + 0.5f*iDetColCount*p->fDetVX;
+			out[ 4*iAngleCount + i] = p->fDetSY + 0.5f*iDetRowCount*p->fDetUY + 0.5f*iDetColCount*p->fDetVY;
+			out[ 5*iAngleCount + i] = p->fDetSZ + 0.5f*iDetRowCount*p->fDetUZ + 0.5f*iDetColCount*p->fDetVZ;
+			out[ 6*iAngleCount + i] = p->fDetUX;
+			out[ 7*iAngleCount + i] = p->fDetUY;
+			out[ 8*iAngleCount + i] = p->fDetUZ;
+			out[ 9*iAngleCount + i] = p->fDetVX;
+			out[10*iAngleCount + i] = p->fDetVY;
+			out[11*iAngleCount + i] = p->fDetVZ;
+		}
+		mGeometryInfo["Vectors"] = pVectors;
+	}
+
+	// cone beam
+	else if (_pProjGeom->isOfType("cone")) {
+		astra::CConeProjectionGeometry3D* pConeGeom = dynamic_cast<astra::CConeProjectionGeometry3D*>(_pProjGeom);
+
+		mGeometryInfo["type"] = mxCreateString("cone");
+		mGeometryInfo["DetectorRowCount"] = mxCreateDoubleScalar(pConeGeom->getDetectorRowCount());
+		mGeometryInfo["DetectorColCount"] = mxCreateDoubleScalar(pConeGeom->getDetectorColCount());
+		mGeometryInfo["DetectorSpacingX"] = mxCreateDoubleScalar(pConeGeom->getDetectorSpacingX());
+		mGeometryInfo["DetectorSpacingY"] = mxCreateDoubleScalar(pConeGeom->getDetectorSpacingY());
+		mGeometryInfo["DistanceOriginSource"] = mxCreateDoubleScalar(pConeGeom->getOriginSourceDistance());
+		mGeometryInfo["DistanceOriginDetector"] = mxCreateDoubleScalar(pConeGeom->getOriginDetectorDistance());	
+
+		mxArray* pAngles = mxCreateDoubleMatrix(1, pConeGeom->getProjectionCount(), mxREAL);
+		double* out = mxGetPr(pAngles);
+		for (int i = 0; i < pConeGeom->getProjectionCount(); i++) {
+			out[i] = pConeGeom->getProjectionAngle(i);
+		}
+		mGeometryInfo["ProjectionAngles"] = pAngles;
+	}
+
+	// cone beam vector
+	else if (_pProjGeom->isOfType("cone_vec")) {
+		astra::CConeVecProjectionGeometry3D* pConeVecGeom = dynamic_cast<astra::CConeVecProjectionGeometry3D*>(_pProjGeom);
+
+		mGeometryInfo["type"] = mxCreateString("cone_vec");
+		mGeometryInfo["DetectorRowCount"] = mxCreateDoubleScalar(pConeVecGeom->getDetectorRowCount());
+		mGeometryInfo["DetectorColCount"] = mxCreateDoubleScalar(pConeVecGeom->getDetectorColCount());
+
+		mxArray* pVectors = mxCreateDoubleMatrix(pConeVecGeom->getProjectionCount(), 12, mxREAL);
+		double* out = mxGetPr(pVectors);
+		int iDetRowCount = pConeVecGeom->getDetectorRowCount();
+		int iDetColCount = pConeVecGeom->getDetectorColCount();
+		int iAngleCount = pConeVecGeom->getProjectionCount();
+		for (int i = 0; i < pConeVecGeom->getProjectionCount(); i++) {
+			const SConeProjection* p = &pConeVecGeom->getProjectionVectors()[i];
+			out[ 0*iAngleCount + i] = p->fSrcX;
+			out[ 1*iAngleCount + i] = p->fSrcY;
+			out[ 2*iAngleCount + i] = p->fSrcZ;
+			out[ 3*iAngleCount + i] = p->fDetSX + 0.5f*iDetRowCount*p->fDetUX + 0.5f*iDetColCount*p->fDetVX;
+			out[ 4*iAngleCount + i] = p->fDetSY + 0.5f*iDetRowCount*p->fDetUY + 0.5f*iDetColCount*p->fDetVY;
+			out[ 5*iAngleCount + i] = p->fDetSZ + 0.5f*iDetRowCount*p->fDetUZ + 0.5f*iDetColCount*p->fDetVZ;
+			out[ 6*iAngleCount + i] = p->fDetUX;
+			out[ 7*iAngleCount + i] = p->fDetUY;
+			out[ 8*iAngleCount + i] = p->fDetUZ;
+			out[ 9*iAngleCount + i] = p->fDetVX;
+			out[10*iAngleCount + i] = p->fDetVY;
+			out[11*iAngleCount + i] = p->fDetVZ;
+		}
+		mGeometryInfo["Vectors"] = pVectors;
+	}
+
+	// build and return the MATLAB struct
+	return buildStruct(mGeometryInfo);
+}
+
+//-----------------------------------------------------------------------------------------
 // parse reconstruction geometry data
 astra::CVolumeGeometry2D* parseVolumeGeometryStruct(const mxArray* prhs)
 {
diff --git a/matlab/mex/mexHelpFunctions.h b/matlab/mex/mexHelpFunctions.h
index ae8acac..8b65a04 100644
--- a/matlab/mex/mexHelpFunctions.h
+++ b/matlab/mex/mexHelpFunctions.h
@@ -45,8 +45,12 @@ $Id$
 
 #include "astra/ParallelProjectionGeometry2D.h"
 #include "astra/FanFlatProjectionGeometry2D.h"
-#include "astra/VolumeGeometry2D.h"
+#include "astra/ParallelProjectionGeometry3D.h"
+#include "astra/ParallelVecProjectionGeometry3D.h"
+#include "astra/ConeProjectionGeometry3D.h"
+#include "astra/ConeVecProjectionGeometry3D.h"
 
+#include "astra/VolumeGeometry2D.h"
 #include "astra/VolumeGeometry3D.h"
 
 
@@ -65,10 +69,11 @@ mxArray* vectorToMxArray(std::vector<astra::float32> mInput);
 mxArray* anyToMxArray(boost::any _any);
 
 astra::CProjectionGeometry2D* parseProjectionGeometryStruct(const mxArray*);
-mxArray* createProjectionGeometryStruct(astra::CProjectionGeometry2D*);
-
 astra::CVolumeGeometry2D* parseVolumeGeometryStruct(const mxArray*);
 
+mxArray* createProjectionGeometryStruct(astra::CProjectionGeometry2D*);
+mxArray* createProjectionGeometryStruct(astra::CProjectionGeometry3D*);
+
 mxArray* createVolumeGeometryStruct(astra::CVolumeGeometry2D* _pReconGeom);
 mxArray* createVolumeGeometryStruct(astra::CVolumeGeometry3D* _pReconGeom);