From 569515f3e20ef3b3c2c4a777f38f45dc67e6f9b6 Mon Sep 17 00:00:00 2001 From: Wim van Aarle 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(pDataObject); - // plhs[0] = createProjectionGeometryStruct(pDataObject2->getGeometry()); - mexErrMsgTxt("Not implemented yet. \n"); + CFloat32ProjectionData3DMemory* pDataObject2 = dynamic_cast(pDataObject); + plhs[0] = createProjectionGeometryStruct(pDataObject2->getGeometry()); } else if (pDataObject->getType() == CFloat32Data3D::VOLUME) { CFloat32VolumeData3DMemory* pDataObject2 = dynamic_cast(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 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(_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(_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(_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 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