From 569515f3e20ef3b3c2c4a777f38f45dc67e6f9b6 Mon Sep 17 00:00:00 2001
From: Wim van Aarle <wimvanaarle@gmail.com>
Date: Tue, 24 Feb 2015 16:46:39 +0100
Subject: added get_geometry for 3d projection data objects
---
matlab/mex/astra_mex_data3d_c.cpp | 5 +-
matlab/mex/mexHelpFunctions.cpp | 117 +++++++++++++++++++++++++++++++++--
matlab/mex/mexHelpFunctions.h | 11 +++-
src/ConeVecProjectionGeometry3D.cpp | 4 +-
src/ParallelProjectionGeometry3D.cpp | 2 +-
test_geometry.m | 22 ++++---
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;
}
@@ -279,6 +279,115 @@ mxArray* createProjectionGeometryStruct(astra::CProjectionGeometry2D* _pProjGeom
return buildStruct(mGeometryInfo);
}
+//-----------------------------------------------------------------------------------------
+// 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
--
cgit v1.2.3