added get_geometry for 3d projection data objects

6 files changed, 141 insertions, 20 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);
 
diff --git a/src/ConeVecProjectionGeometry3D.cpp b/src/ConeVecProjectionGeometry3D.cpp
index cf8f76d..fd549e0 100644
--- a/src/ConeVecProjectionGeometry3D.cpp
+++ b/src/ConeVecProjectionGeometry3D.cpp
@@ -198,13 +198,13 @@ bool CConeVecProjectionGeometry3D::isEqual(const CProjectionGeometry3D * _pGeom2
 // is of type
 bool CConeVecProjectionGeometry3D::isOfType(const std::string& _sType) const
 {
-	 return (_sType == "cone3d_vec");
+	 return (_sType == "cone3d_vec" || _sType == "cone_vec");
 }
 
 //----------------------------------------------------------------------------------------
 void CConeVecProjectionGeometry3D::toXML(XMLNode* _sNode) const
 {
-	_sNode->addAttribute("type","cone3d_vec");
+	_sNode->addAttribute("type","cone_vec");
 	_sNode->addChildNode("DetectorRowCount", m_iDetectorRowCount);
 	_sNode->addChildNode("DetectorColCount", m_iDetectorColCount);
 	// TODO:
diff --git a/src/ParallelProjectionGeometry3D.cpp b/src/ParallelProjectionGeometry3D.cpp
index 8a13f0f..2400ff9 100644
--- a/src/ParallelProjectionGeometry3D.cpp
+++ b/src/ParallelProjectionGeometry3D.cpp
@@ -154,7 +154,7 @@ bool CParallelProjectionGeometry3D::isEqual(const CProjectionGeometry3D * _pGeom
 // is of type
 bool CParallelProjectionGeometry3D::isOfType(const std::string& _sType) const
 {
-	 return (_sType == "parallel");
+	 return (_sType == "parallel" || _sType == "parallel3d");
 }
 
 //----------------------------------------------------------------------------------------
diff --git a/test_geometry.m b/test_geometry.m
index f846d52..8952660 100644
--- a/test_geometry.m
+++ b/test_geometry.m
@@ -15,14 +15,23 @@ proj_geom_parallel = astra_create_proj_geom('parallel', 1, 64, 1:180);
 proj_geom_fanflat = astra_create_proj_geom('fanflat', 1, 64, 1:180, 0, 2000);
 proj_geom_fanflat_vec = astra_geom_2vec(proj_geom_fanflat);
 
+proj_geom_parallel3d = astra_create_proj_geom('parallel3d', 1, 1, 64, 64, 1:180);
+proj_geom_cone = astra_create_proj_geom('cone', 1, 1, 64, 64, 1:180, 0, 2000);
+proj_geom_parallel3d_vec = astra_geom_2vec(proj_geom_parallel3d);
+proj_geom_cone_vec = astra_geom_2vec(proj_geom_cone);
+
 %% create data objects
 vol2d_id = astra_mex_data2d('create', '-vol', vol_geom2d, 0);
 vol3d_id = astra_mex_data3d('create', '-vol', vol_geom3d, 0);
 
 proj_parallel_id = astra_mex_data2d('create', '-sino', proj_geom_parallel, 0);
 proj_fanflat_id = astra_mex_data2d('create', '-sino', proj_geom_fanflat, 0);
-proj_fanflatvec_id = astra_mex_data2d('create', '-sino', proj_geom_fanflat_vec, 0);
+proj_fanflat_vec_id = astra_mex_data2d('create', '-sino', proj_geom_fanflat_vec, 0);
 
+proj_parallel3d_id = astra_mex_data3d('create', '-sino', proj_geom_parallel3d, 0);
+proj_cone_id = astra_mex_data3d('create', '-sino', proj_geom_cone, 0);
+proj_parallel3d_vec_id = astra_mex_data3d('create', '-sino', proj_geom_parallel3d_vec, 0);
+proj_cone_vec_id = astra_mex_data3d('create', '-sino', proj_geom_cone_vec, 0);
 
 %% get geometries
 vol_geom2d_new = astra_mex_data2d('get_geometry', vol2d_id);
@@ -30,13 +39,12 @@ vol_geom3d_new = astra_mex_data3d('get_geometry', vol3d_id);
 
 proj_geom_parallel_new = astra_mex_data2d('get_geometry', proj_parallel_id);
 proj_geom_fanflat_new = astra_mex_data2d('get_geometry', proj_fanflat_id);
-proj_geom_fanflat_vec_new = astra_mex_data2d('get_geometry', proj_fanflatvec_id);
-
-proj_geom_fanflat_vec
-proj_geom_fanflat_vec_new
+proj_geom_fanflat_vec_new = astra_mex_data2d('get_geometry', proj_fanflat_vec_id);
 
-proj_geom_fanflat_vec.Vectors(110,:)
-proj_geom_fanflat_vec_new.Vectors(110,:)
+proj_geom_parallel3d_new = astra_mex_data3d('get_geometry', proj_parallel3d_id);
+proj_geom_cone_new = astra_mex_data3d('get_geometry', proj_cone_id);
+proj_geom_parallel3d_vec_new = astra_mex_data3d('get_geometry', proj_parallel3d_vec_id);
+proj_geom_cone_vec_new = astra_mex_data3d('get_geometry', proj_cone_vec_id);
 
 %%
 astra_clear;
\ No newline at end of file