From b2fc6c70434674d74551c3a6c01ffb3233499312 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Mon, 1 Jul 2013 22:34:11 +0000 Subject: Update version to 1.3 --- matlab/mex/astra_mex.cpp | 121 +++ matlab/mex/astra_mex_algorithm_c.cpp | 348 +++++++++ matlab/mex/astra_mex_algorithm_vc08.vcproj | 593 +++++++++++++++ matlab/mex/astra_mex_c.cpp | 127 ++++ matlab/mex/astra_mex_data2d_c.cpp | 667 +++++++++++++++++ matlab/mex/astra_mex_data2d_vc08.vcproj | 591 +++++++++++++++ matlab/mex/astra_mex_data3d_c.cpp | 1036 ++++++++++++++++++++++++++ matlab/mex/astra_mex_data3d_vc08.vcproj | 588 +++++++++++++++ matlab/mex/astra_mex_matrix_c.cpp | 437 +++++++++++ matlab/mex/astra_mex_matrix_vc08.vcproj | 591 +++++++++++++++ matlab/mex/astra_mex_projector3d_c.cpp | 433 +++++++++++ matlab/mex/astra_mex_projector3d_vc08.vcproj | 588 +++++++++++++++ matlab/mex/astra_mex_projector_c.cpp | 510 +++++++++++++ matlab/mex/astra_mex_projector_vc08.vcproj | 591 +++++++++++++++ matlab/mex/astra_mex_vc08.vcproj | 591 +++++++++++++++ matlab/mex/mex.def | 1 + matlab/mex/mexHelpFunctions.cpp | 642 ++++++++++++++++ matlab/mex/mexHelpFunctions.h | 76 ++ 18 files changed, 8531 insertions(+) create mode 100644 matlab/mex/astra_mex.cpp create mode 100644 matlab/mex/astra_mex_algorithm_c.cpp create mode 100644 matlab/mex/astra_mex_algorithm_vc08.vcproj create mode 100644 matlab/mex/astra_mex_c.cpp create mode 100644 matlab/mex/astra_mex_data2d_c.cpp create mode 100644 matlab/mex/astra_mex_data2d_vc08.vcproj create mode 100644 matlab/mex/astra_mex_data3d_c.cpp create mode 100644 matlab/mex/astra_mex_data3d_vc08.vcproj create mode 100644 matlab/mex/astra_mex_matrix_c.cpp create mode 100644 matlab/mex/astra_mex_matrix_vc08.vcproj create mode 100644 matlab/mex/astra_mex_projector3d_c.cpp create mode 100644 matlab/mex/astra_mex_projector3d_vc08.vcproj create mode 100644 matlab/mex/astra_mex_projector_c.cpp create mode 100644 matlab/mex/astra_mex_projector_vc08.vcproj create mode 100644 matlab/mex/astra_mex_vc08.vcproj create mode 100644 matlab/mex/mex.def create mode 100644 matlab/mex/mexHelpFunctions.cpp create mode 100644 matlab/mex/mexHelpFunctions.h (limited to 'matlab/mex') diff --git a/matlab/mex/astra_mex.cpp b/matlab/mex/astra_mex.cpp new file mode 100644 index 0000000..4b77f76 --- /dev/null +++ b/matlab/mex/astra_mex.cpp @@ -0,0 +1,121 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +$Id$ +*/ + +#include +#include "mexHelpFunctions.h" + +#include "astra/Globals.h" + +using namespace std; +using namespace astra; + + +//----------------------------------------------------------------------------------------- +/** astra_mex('credits'); + * + * Print Credits + */ +void astra_mex_credits(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + cout << "All Scale Tomographic Reconstruction Antwerp Toolbox (ASTRA-Toolbox) was developed at the University of Antwerp by" << endl; + cout << " * Joost Batenburg, PhD" << endl; + cout << " * Gert Merckx" << endl; + cout << " * Willem Jan Palenstijn" << endl; + cout << " * Tom Roelandts" << endl; + cout << " * Prof. Dr. Jan Sijbers" << endl; + cout << " * Wim van Aarle" << endl; + cout << " * Sander van der Maar" << endl; + cout << " * Gert Van Gompel, PhD" << endl; +} + +//----------------------------------------------------------------------------------------- +/** use_cuda = astra_mex('use_cuda'); + * + * Is CUDA enabled? + */ +void astra_mex_use_cuda(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(astra::cudaEnabled() ? 1 : 0); + } +} + +//----------------------------------------------------------------------------------------- +/** version_number = astra_mex('version'); + * + * Fetch the version number of the toolbox. + */ +void astra_mex_version(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(astra::getVersion()); + } else { + cout << "astra toolbox version " << astra::getVersionString() << endl; + } +} + +//----------------------------------------------------------------------------------------- + +static void printHelp() +{ + mexPrintf("Please specify a mode of operation.\n"); + mexPrintf(" Valid modes: version, use_cuda, credits\n"); +} + +//----------------------------------------------------------------------------------------- +/** + * ... = astra_mex(type,...); + */ +void mexFunction(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + + // INPUT0: Mode + string sMode = ""; + if (1 <= nrhs) { + sMode = mex_util_get_string(prhs[0]); + } else { + printHelp(); + return; + } + + // SWITCH (MODE) + if (sMode == std::string("version")) { + astra_mex_version(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("use_cuda")) { + astra_mex_use_cuda(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("credits")) { + astra_mex_credits(nlhs, plhs, nrhs, prhs); + } else { + printHelp(); + } + + return; +} + + diff --git a/matlab/mex/astra_mex_algorithm_c.cpp b/matlab/mex/astra_mex_algorithm_c.cpp new file mode 100644 index 0000000..7476ba4 --- /dev/null +++ b/matlab/mex/astra_mex_algorithm_c.cpp @@ -0,0 +1,348 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +$Id$ +*/ + +/** \file astra_mex_algorithm_c.cpp + * + * \brief Creates and manages algorithms (reconstruction,projection,...). + */ +#include +#include "mexHelpFunctions.h" + +#define USE_MATLAB_UNDOCUMENTED + +#ifdef USE_MATLAB_UNDOCUMENTED +extern "C" { bool utIsInterruptPending(); } + +#ifdef __linux__ +#define USE_PTHREADS_CTRLC +#include +#else +#include +#endif + +#endif + + + +#include "astra/Globals.h" + +#include "astra/AstraObjectManager.h" +#include "astra/AstraObjectFactory.h" + +#include "astra/XMLNode.h" +#include "astra/XMLDocument.h" + +using namespace std; +using namespace astra; +//----------------------------------------------------------------------------------------- +/** id = astra_mex_algorithm('create', cfg); + * + * Create and configure a new algorithm object. + * cfg: MATLAB struct containing the configuration parameters, see doxygen documentation for details. + * id: identifier of the algorithm object as it is now stored in the astra-library. + */ +void astra_mex_algorithm_create(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + + if (!mxIsStruct(prhs[1])) { + mexErrMsgTxt("Argument 1 not a valid MATLAB struct. \n"); + } + + // turn MATLAB struct to an XML-based Config object + XMLDocument* xml = struct2XML("Algorithm", prhs[1]); + Config cfg; + cfg.self = xml->getRootNode(); + + CAlgorithm* pAlg = CAlgorithmFactory::getSingleton().create(cfg.self->getAttribute("type")); + if (!pAlg) { + delete xml; + mexErrMsgTxt("Unknown algorithm. \n"); + return; + } + + // create algorithm + if (!pAlg->initialize(cfg)) { + delete xml; + delete pAlg; + mexErrMsgTxt("Algorithm not initialized. \n"); + return; + } + + delete xml; + + // store algorithm + int iIndex = CAlgorithmManager::getSingleton().store(pAlg); + + // step4: set output + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(iIndex); + } + +} + +#ifdef USE_MATLAB_UNDOCUMENTED + +#ifndef USE_PTHREADS_CTRLC + +// boost version +void waitForInterrupt_boost(CAlgorithm* _pAlg) +{ + boost::posix_time::milliseconds rel(2000); + + while (!utIsInterruptPending()) { + + // This is an interruption point. If the main thread calls + // interrupt(), this thread will terminate here. + boost::this_thread::sleep(rel); + } + + //mexPrintf("Aborting. Please wait.\n"); + + // One last quick check to see if the algorithm already finished + boost::this_thread::interruption_point(); + + _pAlg->signalAbort(); +} + +#else + +// pthreads version +void *waitForInterrupt_pthreads(void *threadid) +{ + CAlgorithm* _pAlg = (CAlgorithm*)threadid; + + while (!utIsInterruptPending()) { + usleep(50000); + pthread_testcancel(); + } + + //mexPrintf("Aborting. Please wait.\n"); + + // One last quick check to see if the algorithm already finished + pthread_testcancel(); + + _pAlg->signalAbort(); + + return 0; +} + +#endif +#endif + +//----------------------------------------------------------------------------------------- +/** astra_mex_algorithm('run', id); or astra_mex_algorithm('iterate', id, iterations); + * + * Run or do iterations on a certain algorithm. + * id: identifier of the algorithm object as stored in the astra-library. + * iterations: if the algorithm is iterative, this specifies the number of iterations to perform. + */ +void astra_mex_algorithm_run(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: get input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iAid = (int)(mxGetScalar(prhs[1])); + int iIterations = 0; + if (3 <= nrhs) { + iIterations = (int)(mxGetScalar(prhs[2])); + } + + // step2: get algorithm object + CAlgorithm* pAlg = CAlgorithmManager::getSingleton().get(iAid); + if (!pAlg) { + mexErrMsgTxt("Invalid algorithm ID.\n"); + return; + } + if (!pAlg->isInitialized()) { + mexErrMsgTxt("Algorithm not initialized. \n"); + return; + } + + // step3: perform actions +#ifndef USE_MATLAB_UNDOCUMENTED + + pAlg->run(iIterations); + +#elif defined(USE_PTHREADS_CTRLC) + + // Start a new thread to watch if the user pressed Ctrl-C + pthread_t thread; + pthread_create(&thread, 0, waitForInterrupt_pthreads, (void*)pAlg); + + pAlg->run(iIterations); + + // kill the watcher thread in case it's still running + pthread_cancel(thread); + pthread_join(thread, 0); + +#else + + // Start a new thread to watch if the user pressed Ctrl-C + boost::thread interruptThread(waitForInterrupt_boost, pAlg); + + pAlg->run(iIterations); + + // kill the watcher thread in case it's still running + interruptThread.interrupt(); + interruptThread.join(); + +#endif +} +//----------------------------------------------------------------------------------------- +/** astra_mex_algorithm('get_res_norm', id); + * + * Get the L2-norm of the residual sinogram. Not all algorithms + * support this operation. + */ +void astra_mex_algorithm_get_res_norm(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: get input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iAid = (int)(mxGetScalar(prhs[1])); + + // step2: get algorithm object + CAlgorithm* pAlg = CAlgorithmManager::getSingleton().get(iAid); + if (!pAlg) { + mexErrMsgTxt("Invalid algorithm ID.\n"); + return; + } + if (!pAlg->isInitialized()) { + mexErrMsgTxt("Algorithm not initialized. \n"); + return; + } + + CReconstructionAlgorithm2D* pAlg2D = dynamic_cast(pAlg); + CReconstructionAlgorithm3D* pAlg3D = dynamic_cast(pAlg); + + float res = 0.0f; + bool ok; + if (pAlg2D) + ok = pAlg2D->getResidualNorm(res); + else if (pAlg3D) + ok = pAlg3D->getResidualNorm(res); + else + ok = false; + + if (!ok) { + mexErrMsgTxt("Operation not supported.\n"); + return; + } + + plhs[0] = mxCreateDoubleScalar(res); +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_algorithm('delete', id1, id2, ...); + * + * Delete one or more algorithm objects currently stored in the astra-library. + * id1, id2, ... : identifiers of the algorithm objects as stored in the astra-library. + */ +void astra_mex_algorithm_delete(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: get algorithm ID + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + + for (int i = 1; i < nrhs; i++) { + int iAid = (int)(mxGetScalar(prhs[i])); + CAlgorithmManager::getSingleton().remove(iAid); + } +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_algorithm('clear'); + * + * Delete all algorithm objects currently stored in the astra-library. + */ +void astra_mex_algorithm_clear(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + CAlgorithmManager::getSingleton().clear(); +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_algorithm('info'); + * + * Print information about all the algorithm objects currently stored in the astra-library. + */ +void astra_mex_algorithm_info(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + mexPrintf("%s", astra::CAlgorithmManager::getSingleton().info().c_str()); +} + +//----------------------------------------------------------------------------------------- +static void printHelp() +{ + mexPrintf("Please specify a mode of operation.\n"); + mexPrintf("Valid modes: create, info, delete, clear, run/iterate, get_res_norm\n"); +} + +//----------------------------------------------------------------------------------------- +/** + * ... = astra_mex_algorithm(mode, ...); + */ +void mexFunction(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + // INPUT: Mode + string sMode = ""; + if (1 <= nrhs) { + sMode = mex_util_get_string(prhs[0]); + } else { + printHelp(); + return; + } + + // SWITCH (MODE) + if (sMode == "create") { + astra_mex_algorithm_create(nlhs, plhs, nrhs, prhs); + } else if (sMode == "info") { + astra_mex_algorithm_info(nlhs, plhs, nrhs, prhs); + } else if (sMode == "delete") { + astra_mex_algorithm_delete(nlhs, plhs, nrhs, prhs); + } else if (sMode == "clear") { + astra_mex_algorithm_clear(nlhs, plhs, nrhs, prhs); + } else if (sMode == "run" || sMode == "iterate") { + astra_mex_algorithm_run(nlhs, plhs, nrhs, prhs); + } else if (sMode == "get_res_norm") { + astra_mex_algorithm_get_res_norm(nlhs, plhs, nrhs, prhs); + } else { + printHelp(); + } + return; +} diff --git a/matlab/mex/astra_mex_algorithm_vc08.vcproj b/matlab/mex/astra_mex_algorithm_vc08.vcproj new file mode 100644 index 0000000..baa4c44 --- /dev/null +++ b/matlab/mex/astra_mex_algorithm_vc08.vcproj @@ -0,0 +1,593 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/matlab/mex/astra_mex_c.cpp b/matlab/mex/astra_mex_c.cpp new file mode 100644 index 0000000..0068664 --- /dev/null +++ b/matlab/mex/astra_mex_c.cpp @@ -0,0 +1,127 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +$Id$ +*/ + +/** \file astra_mex_c.cpp + * + * \brief Contains some basic "about" functions. + */ + +#include +#include "mexHelpFunctions.h" + +#include "astra/Globals.h" + +using namespace std; +using namespace astra; + + +//----------------------------------------------------------------------------------------- +/** astra_mex('credits'); + * + * Print Credits + */ +void astra_mex_credits(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + mexPrintf("All Scale Tomographic Reconstruction Antwerp Toolbox (ASTRA-Toolbox) was developed at the University of Antwerp by\n"); + mexPrintf(" * Prof. dr. Joost Batenburg\n"); + mexPrintf(" * Andrei Dabravolski\n"); + mexPrintf(" * Gert Merckx\n"); + mexPrintf(" * Willem Jan Palenstijn\n"); + mexPrintf(" * Tom Roelandts\n"); + mexPrintf(" * Prof. dr. Jan Sijbers\n"); + mexPrintf(" * dr. Wim van Aarle\n"); + mexPrintf(" * Sander van der Maar\n"); + mexPrintf(" * dr. Gert Van Gompel\n"); +} + +//----------------------------------------------------------------------------------------- +/** use_cuda = astra_mex('use_cuda'); + * + * Is CUDA enabled? + */ +void astra_mex_use_cuda(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(astra::cudaEnabled() ? 1 : 0); + } +} + +//----------------------------------------------------------------------------------------- +/** version_number = astra_mex('version'); + * + * Fetch the version number of the toolbox. + */ +void astra_mex_version(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(astra::getVersion()); + } else { + mexPrintf("astra toolbox version %s\n", astra::getVersionString()); + } +} + +//----------------------------------------------------------------------------------------- + +static void printHelp() +{ + mexPrintf("Please specify a mode of operation.\n"); + mexPrintf(" Valid modes: version, use_cuda, credits\n"); +} + +//----------------------------------------------------------------------------------------- +/** + * ... = astra_mex(type,...); + */ +void mexFunction(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + + // INPUT0: Mode + string sMode = ""; + if (1 <= nrhs) { + sMode = mex_util_get_string(prhs[0]); + } else { + printHelp(); + return; + } + + // SWITCH (MODE) + if (sMode == std::string("version")) { + astra_mex_version(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("use_cuda")) { + astra_mex_use_cuda(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("credits")) { + astra_mex_credits(nlhs, plhs, nrhs, prhs); + } else { + printHelp(); + } + + return; +} + + diff --git a/matlab/mex/astra_mex_data2d_c.cpp b/matlab/mex/astra_mex_data2d_c.cpp new file mode 100644 index 0000000..99fb38e --- /dev/null +++ b/matlab/mex/astra_mex_data2d_c.cpp @@ -0,0 +1,667 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +$Id$ +*/ + +/** \file astra_mex_data2d_c.cpp + * + * \brief Creates, manages and manipulates 2D volume and projection data objects. + */ +#include +#include "mexHelpFunctions.h" + +#include + +#include "astra/Globals.h" + +#include "astra/AstraObjectManager.h" + +#include "astra/Float32ProjectionData2D.h" +#include "astra/Float32VolumeData2D.h" +#include "astra/SparseMatrixProjectionGeometry2D.h" +#include "astra/FanFlatProjectionGeometry2D.h" +#include "astra/FanFlatVecProjectionGeometry2D.h" + +using namespace std; +using namespace astra; + +//----------------------------------------------------------------------------------------- +/** astra_mex_data2d('delete', id1, id2, ...); + * + * Delete one or more data objects currently stored in the astra-library. + * id1, id2, ... : identifiers of the 2d data objects as stored in the astra-library. + */ +void astra_mex_data2d_delete(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: read input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + + // step2: delete all specified data objects + for (int i = 1; i < nrhs; i++) { + int iDataID = (int)(mxGetScalar(prhs[i])); + CData2DManager::getSingleton().remove(iDataID); + } +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_data2d('clear'); + * + * Delete all data objects currently stored in the astra-library. + */ +void astra_mex_data2d_clear(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + CData2DManager::getSingleton().clear(); +} + +//----------------------------------------------------------------------------------------- +/** id = astra_mex_data2d('create', datatype, geometry, data); + * + * Create a new data 2d object in the astra-library. + * type: '-vol' for volume data, '-sino' for projection data + * geom: MATLAB struct with the geometry for the data + * data: Optional. Can be either a MATLAB matrix containing the data. In that case the dimensions + * should match that of the geometry of the object. It can also be a single value, in which case + * the entire data will be set to that value. If this isn't specified all values are set to 0. + * id: identifier of the 2d data object as it is now stored in the astra-library. + */ +void astra_mex_data2d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArray* prhs[]) +{ + // step1: get datatype + if (nrhs < 3) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + + string sDataType = mex_util_get_string(prhs[1]); + CFloat32Data2D* pDataObject2D = NULL; + + if (nrhs >= 4 && !(mex_is_scalar(prhs[3])|| mxIsDouble(prhs[3]) || mxIsLogical(prhs[3]) || mxIsSingle(prhs[3]) )) { + mexErrMsgTxt("Data must be single, double or logical."); + return; + } + + // SWITCH DataType + if (sDataType == "-vol") { + // Read geometry + if (!mxIsStruct(prhs[2])) { + mexErrMsgTxt("Argument 3 is not a valid MATLAB struct.\n"); + } + XMLDocument* xml = struct2XML(string("VolumeGeometry"), prhs[2]); + if (!xml) + return; + Config cfg; + cfg.self = xml->getRootNode(); + CVolumeGeometry2D* pGeometry = new CVolumeGeometry2D(); + if (!pGeometry->initialize(cfg)) { + mexErrMsgTxt("Geometry class not initialized. \n"); + delete xml; + delete pGeometry; + return; + } + // If data is specified, check dimensions + if (nrhs >= 4 && !mex_is_scalar(prhs[3])) { + if (pGeometry->getGridColCount() != mxGetN(prhs[3]) || pGeometry->getGridRowCount() != mxGetM(prhs[3])) { + mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); + delete xml; + delete pGeometry; + return; + } + } + // Initialize data object + pDataObject2D = new CFloat32VolumeData2D(pGeometry); + delete pGeometry; + delete xml; + } + else if (sDataType == "-sino") { + // Read geometry + if (!mxIsStruct(prhs[2])) { + mexErrMsgTxt("Argument 3 is not a valid MATLAB struct.\n"); + } + XMLDocument* xml = struct2XML("ProjectionGeometry", prhs[2]); + if (!xml) + return; + Config cfg; + cfg.self = xml->getRootNode(); + // FIXME: Change how the base class is created. (This is duplicated + // in 'change_geometry' and Projector2D.cpp.) + std::string type = cfg.self->getAttribute("type"); + CProjectionGeometry2D* pGeometry; + if (type == "sparse_matrix") { + pGeometry = new CSparseMatrixProjectionGeometry2D(); + } else if (type == "fanflat") { + //CFanFlatProjectionGeometry2D* pFanFlatProjectionGeometry = new CFanFlatProjectionGeometry2D(); + //pFanFlatProjectionGeometry->initialize(Config(node)); + //m_pProjectionGeometry = pFanFlatProjectionGeometry; + pGeometry = new CFanFlatProjectionGeometry2D(); + } else if (type == "fanflat_vec") { + pGeometry = new CFanFlatVecProjectionGeometry2D(); + } else { + pGeometry = new CParallelProjectionGeometry2D(); + } + if (!pGeometry->initialize(cfg)) { + mexErrMsgTxt("Geometry class not initialized. \n"); + delete pGeometry; + delete xml; + return; + } + // If data is specified, check dimensions + if (nrhs >= 4 && !mex_is_scalar(prhs[3])) { + if (pGeometry->getDetectorCount() != mxGetN(prhs[3]) || pGeometry->getProjectionAngleCount() != mxGetM(prhs[3])) { + mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); + delete pGeometry; + delete xml; + return; + } + } + // Initialize data object + pDataObject2D = new CFloat32ProjectionData2D(pGeometry); + delete pGeometry; + delete xml; + } + else { + mexErrMsgTxt("Invalid datatype. Please specify '-vol' or '-sino'. \n"); + return; + } + + // Check initialization + if (!pDataObject2D->isInitialized()) { + mexErrMsgTxt("Couldn't initialize data object.\n"); + delete pDataObject2D; + return; + } + + // Store data + if (nrhs == 3) { + for (int i = 0; i < pDataObject2D->getSize(); ++i) { + pDataObject2D->getData()[i] = 0.0f; + } + } + + // Store data + if (nrhs >= 4) { + // fill with scalar value + if (mex_is_scalar(prhs[3])) { + float32 fValue = (float32)mxGetScalar(prhs[3]); + for (int i = 0; i < pDataObject2D->getSize(); ++i) { + pDataObject2D->getData()[i] = fValue; + } + } + // fill with array value + else { + const mwSize* dims = mxGetDimensions(prhs[3]); + // Check Data dimensions + if (pDataObject2D->getWidth() != mxGetN(prhs[3]) || pDataObject2D->getHeight() != mxGetM(prhs[3])) { + mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); + return; + } + + // logical data + if (mxIsLogical(prhs[3])) { + bool* pbMatlabData = mxGetLogicals(prhs[3]); + int i = 0; + int col, row; + for (col = 0; col < dims[1]; ++col) { + for (row = 0; row < dims[0]; ++row) { + pDataObject2D->getData2D()[row][col] = (float32)pbMatlabData[i]; + ++i; + } + } + // double data + } else if (mxIsDouble(prhs[3])) { + double* pdMatlabData = mxGetPr(prhs[3]); + int i = 0; + int col, row; + for (col = 0; col < dims[1]; ++col) { + for (row = 0; row < dims[0]; ++row) { + pDataObject2D->getData2D()[row][col] = pdMatlabData[i]; + ++i; + } + } + // single data + } else if (mxIsSingle(prhs[3])) { + const float* pfMatlabData = (const float *)mxGetData(prhs[3]); + int i = 0; + int col, row; + for (col = 0; col < dims[1]; ++col) { + for (row = 0; row < dims[0]; ++row) { + pDataObject2D->getData2D()[row][col] = pfMatlabData[i]; + ++i; + } + } + } else { + ASTRA_ASSERT(false); + } + } + } + + // step4: store data object + int iIndex = CData2DManager::getSingleton().store(pDataObject2D); + + // step5: return data id + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(iIndex); + } + +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_data2d('store', id, data); + * + * Store data in an existing astra 2d dataobject with a MATLAB matrix or with a scalar value. + * id: identifier of the 2d data object as stored in the astra-library. + * data: can be either a MATLAB matrix containing the data. In that case the dimensions should match that of the geometry of the object. It can also be a single value, in which case the entire data will be set to that value. + */ +void astra_mex_data2d_store(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: input + if (nrhs < 3) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + if (!mxIsDouble(prhs[1])) { + mexErrMsgTxt("Identifier should be a scalar value. \n"); + return; + } + int iDataID = (int)(mxGetScalar(prhs[1])); + + if (!(mex_is_scalar(prhs[2]) || mxIsDouble(prhs[2]) || mxIsLogical(prhs[2]) || mxIsSingle(prhs[2]))) { + mexErrMsgTxt("Data must be single, double or logical."); + return; + } + + // step2: get data object + CFloat32Data2D* pDataObject = astra::CData2DManager::getSingleton().get(iDataID); + if (!pDataObject || !pDataObject->isInitialized()) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; + } + + // step3: insert data + // fill with scalar value + if (mex_is_scalar(prhs[2])) { + float32 fValue = (float32)mxGetScalar(prhs[2]); + for (int i = 0; i < pDataObject->getSize(); ++i) { + pDataObject->getData()[i] = fValue; + } + } else { + // Check Data dimensions + if (pDataObject->getWidth() != mxGetN(prhs[2]) || pDataObject->getHeight() != mxGetM(prhs[2])) { + mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); + return; + } + const mwSize* dims = mxGetDimensions(prhs[2]); + + // logical data + if (mxIsLogical(prhs[2])) { + bool* pbMatlabData = mxGetLogicals(prhs[2]); + int i = 0; + int col, row; + for (col = 0; col < dims[1]; ++col) { + for (row = 0; row < dims[0]; ++row) { + pDataObject->getData2D()[row][col] = (float32)pbMatlabData[i]; + ++i; + } + } + // double data + } else if (mxIsDouble(prhs[2])) { + double* pdMatlabData = mxGetPr(prhs[2]); + int i = 0; + int col, row; + for (col = 0; col < dims[1]; ++col) { + for (row = 0; row < dims[0]; ++row) { + pDataObject->getData2D()[row][col] = pdMatlabData[i]; + ++i; + } + } + // single data + } else if (mxIsSingle(prhs[2])) { + const float* pfMatlabData = (const float *)mxGetData(prhs[2]); + int i = 0; + int col, row; + for (col = 0; col < dims[1]; ++col) { + for (row = 0; row < dims[0]; ++row) { + pDataObject->getData2D()[row][col] = pfMatlabData[i]; + ++i; + } + } + } else { + ASTRA_ASSERT(false); + } + } +} + +//----------------------------------------------------------------------------------------- +/** geom = astra_mex_data2d('get_geometry', id); + * + * Fetch the geometry of a 2d data object stored in the astra-library. + * id: identifier of the 2d data object as stored in the astra-library. + * geom: MATLAB-struct containing information about the used geometry. + */ +void astra_mex_data2d_get_geometry(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + if (!mxIsDouble(prhs[1])) { + mexErrMsgTxt("Identifier should be a scalar value. \n"); + return; + } + int iDataID = (int)(mxGetScalar(prhs[1])); + + // step2: get data object + CFloat32Data2D* pDataObject = astra::CData2DManager::getSingleton().get(iDataID); + if (!pDataObject || !pDataObject->isInitialized()) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; + } + + // create output + if (1 <= nlhs) { + if (pDataObject->getType() == CFloat32Data2D::PROJECTION) { + CFloat32ProjectionData2D* pDataObject2 = dynamic_cast(pDataObject); + plhs[0] = createProjectionGeometryStruct(pDataObject2->getGeometry()); + } + else if (pDataObject->getType() == CFloat32Data2D::VOLUME) { + CFloat32VolumeData2D* pDataObject2 = dynamic_cast(pDataObject); + plhs[0] = createVolumeGeometryStruct(pDataObject2->getGeometry()); + } + } +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_data2d('change_geometry', id, geom); + * + * Change the associated geometry of a 2d data object (volume or sinogram) + * id: identifier of the 2d data object as stored in the astra-library. + * geom: the new geometry struct, as created by astra_create_vol/proj_geom + */ +void astra_mex_data2d_change_geometry(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: check input + if (nrhs < 3) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + if (!mxIsDouble(prhs[1])) { + mexErrMsgTxt("Identifier should be a scalar value. \n"); + return; + } + + // step2: get data object + int iDataID = (int)(mxGetScalar(prhs[1])); + CFloat32Data2D* pDataObject = astra::CData2DManager::getSingleton().get(iDataID); + if (!pDataObject || !pDataObject->isInitialized()) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; + } + + CFloat32ProjectionData2D* pSinogram = dynamic_cast(pDataObject); + + if (pSinogram) { + // Projection data + + // Read geometry + if (!mxIsStruct(prhs[2])) { + mexErrMsgTxt("Argument 3 is not a valid MATLAB struct.\n"); + } + XMLDocument* xml = struct2XML("ProjectionGeometry", prhs[2]); + Config cfg; + cfg.self = xml->getRootNode(); + // FIXME: Change how the base class is created. (This is duplicated + // in 'create' and Projector2D.cpp.) + std::string type = cfg.self->getAttribute("type"); + CProjectionGeometry2D* pGeometry; + if (type == "sparse_matrix") { + pGeometry = new CSparseMatrixProjectionGeometry2D(); + } else if (type == "fanflat") { + //CFanFlatProjectionGeometry2D* pFanFlatProjectionGeometry = new CFanFlatProjectionGeometry2D(); + //pFanFlatProjectionGeometry->initialize(Config(node)); + //m_pProjectionGeometry = pFanFlatProjectionGeometry; + pGeometry = new CFanFlatProjectionGeometry2D(); + } else if (type == "fanflat_vec") { + pGeometry = new CFanFlatVecProjectionGeometry2D(); + } else { + pGeometry = new CParallelProjectionGeometry2D(); + } + if (!pGeometry->initialize(cfg)) { + mexErrMsgTxt("Geometry class not initialized. \n"); + delete pGeometry; + delete xml; + return; + } + // If data is specified, check dimensions + if (pGeometry->getDetectorCount() != pSinogram->getDetectorCount() || pGeometry->getProjectionAngleCount() != pSinogram->getAngleCount()) { + mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); + delete pGeometry; + delete xml; + return; + } + + // If ok, change geometry + pSinogram->changeGeometry(pGeometry); + delete pGeometry; + delete xml; + + return; + } + + CFloat32VolumeData2D* pVolume = dynamic_cast(pDataObject); + + if (pVolume) { + // Volume data + + // Read geometry + if (!mxIsStruct(prhs[2])) { + mexErrMsgTxt("Argument 3 is not a valid MATLAB struct.\n"); + } + XMLDocument* xml = struct2XML(string("VolumeGeometry"), prhs[2]); + Config cfg; + cfg.self = xml->getRootNode(); + CVolumeGeometry2D* pGeometry = new CVolumeGeometry2D(); + if (!pGeometry->initialize(cfg)) { + mexErrMsgTxt("Geometry class not initialized. \n"); + delete xml; + delete pGeometry; + return; + } + // If data is specified, check dimensions + if (pGeometry->getGridColCount() != pVolume->getWidth() || pGeometry->getGridRowCount() != pVolume->getHeight()) { + mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); + delete xml; + delete pGeometry; + return; + } + + // If ok, change geometry + pVolume->changeGeometry(pGeometry); + delete xml; + delete pGeometry; + + } + + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; +} + +//----------------------------------------------------------------------------------------- +/** data = astra_mex_data2d('get', id); + * + * Fetch data from the astra-library to a MATLAB matrix. + * id: identifier of the 2d data object as stored in the astra-library. + * data: MATLAB data + */ +void astra_mex_data2d_get(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: check input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + if (!mxIsDouble(prhs[1])) { + mexErrMsgTxt("Identifier should be a scalar value. \n"); + return; + } + + // step2: get data object + int iDataID = (int)(mxGetScalar(prhs[1])); + CFloat32Data2D* pDataObject = astra::CData2DManager::getSingleton().get(iDataID); + if (!pDataObject || !pDataObject->isInitialized()) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; + } + + // create output + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleMatrix(pDataObject->getHeight(), // # rows + pDataObject->getWidth(), // # cols + mxREAL); // datatype 64-bits + double* out = mxGetPr(plhs[0]); + int i = 0; + int row, col; + for (col = 0; col < pDataObject->getWidth(); ++col) { + for (row = 0; row < pDataObject->getHeight(); ++row) { + out[i] = pDataObject->getData2D()[row][col]; + ++i; + } + } + } + +} + +//----------------------------------------------------------------------------------------- +/** data = astra_mex_data2d('get_single', id); + * + * Fetch data from the astra-library to a MATLAB matrix. + * id: identifier of the 2d data object as stored in the astra-library. + * data: MATLAB data + */ +void astra_mex_data2d_get_single(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: check input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + if (!mxIsDouble(prhs[1])) { + mexErrMsgTxt("Identifier should be a scalar value. \n"); + return; + } + + // step2: get data object + int iDataID = (int)(mxGetScalar(prhs[1])); + CFloat32Data2D* pDataObject = astra::CData2DManager::getSingleton().get(iDataID); + if (!pDataObject || !pDataObject->isInitialized()) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; + } + + // create output + if (1 <= nlhs) { + mwSize dims[2]; + dims[0] = pDataObject->getHeight(); + dims[1] = pDataObject->getWidth(); + plhs[0] = mxCreateNumericArray(2, dims, mxSINGLE_CLASS, mxREAL); + float* out = (float *)mxGetData(plhs[0]); + int i = 0; + int row, col; + for (col = 0; col < pDataObject->getWidth(); ++col) { + for (row = 0; row < pDataObject->getHeight(); ++row) { + out[i] = pDataObject->getData2D()[row][col]; + ++i; + } + } + } + +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_data2d('info'); + * + * Print information about all the 2d data objects currently stored in the astra-library. + */ +void astra_mex_data2d_info(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + mexPrintf("%s", astra::CData2DManager::getSingleton().info().c_str()); +} + +//----------------------------------------------------------------------------------------- + +static void printHelp() +{ + mexPrintf("Please specify a mode of operation.\n"); + mexPrintf("Valid modes: get, get_single, delete, clear, set/store, create, get_geometry, change_geometry, info\n"); +} + +//----------------------------------------------------------------------------------------- +/** + * ... = astra_mex_data2d(type,...); + */ +void mexFunction(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + + // INPUT0: Mode + string sMode = ""; + if (1 <= nrhs) { + sMode = mex_util_get_string(prhs[0]); + } else { + printHelp(); + return; + } + + // SWITCH (MODE) + if (sMode == std::string("get")) { + astra_mex_data2d_get(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("get_single")) { + astra_mex_data2d_get_single(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("delete")) { + astra_mex_data2d_delete(nlhs, plhs, nrhs, prhs); + } else if (sMode == "clear") { + astra_mex_data2d_clear(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("store") || + sMode == std::string("set")) { + astra_mex_data2d_store(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("create")) { + astra_mex_data2d_create(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("get_geometry")) { + astra_mex_data2d_get_geometry(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("change_geometry")) { + astra_mex_data2d_change_geometry(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("info")) { + astra_mex_data2d_info(nlhs, plhs, nrhs, prhs); + } else { + printHelp(); + } + + return; +} + + diff --git a/matlab/mex/astra_mex_data2d_vc08.vcproj b/matlab/mex/astra_mex_data2d_vc08.vcproj new file mode 100644 index 0000000..8f1fc13 --- /dev/null +++ b/matlab/mex/astra_mex_data2d_vc08.vcproj @@ -0,0 +1,591 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/matlab/mex/astra_mex_data3d_c.cpp b/matlab/mex/astra_mex_data3d_c.cpp new file mode 100644 index 0000000..1af8844 --- /dev/null +++ b/matlab/mex/astra_mex_data3d_c.cpp @@ -0,0 +1,1036 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +$Id$ +*/ + +/** \file astra_mex_data3d_c.cpp + * + * \brief Creates, manages and manipulates 3D volume and projection data objects. + */ +#include +#include "mexHelpFunctions.h" + +#include + +#include "astra/Globals.h" + +#include "astra/AstraObjectManager.h" + +#include "astra/Float32ProjectionData2D.h" +#include "astra/Float32VolumeData2D.h" +#include "astra/Float32ProjectionData3D.h" +#include "astra/Float32ProjectionData3DMemory.h" +#include "astra/Float32VolumeData3D.h" +#include "astra/Float32VolumeData3DMemory.h" +#include "astra/ParallelProjectionGeometry3D.h" +#include "astra/ParallelVecProjectionGeometry3D.h" +#include "astra/ConeProjectionGeometry3D.h" +#include "astra/ConeVecProjectionGeometry3D.h" + +using namespace std; +using namespace astra; + + + +//----------------------------------------------------------------------------------------- +/** + * id = astra_mex_io_data('create', datatype, geometry, data); + * datatype: ['-vol','-sino','-sinocone'] + */ +void astra_mex_data3d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArray* prhs[]) +{ + // step1: get datatype + if (nrhs < 3) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + + string sDataType = mex_util_get_string(prhs[1]); + CFloat32Data3DMemory* pDataObject3D = NULL; + + if (nrhs >= 4 && !(mex_is_scalar(prhs[3]) || mxIsDouble(prhs[3]) || mxIsSingle(prhs[3]))) { + mexErrMsgTxt("Data must be single or double."); + return; + } + + mwSize dims[3]; + + // SWITCH DataType + if (sDataType == "-vol") { + + // Read geometry + if (!mxIsStruct(prhs[2])) { + mexErrMsgTxt("Argument 3 is not a valid MATLAB struct.\n"); + } + Config cfg; + XMLDocument* xml = struct2XML("VolumeGeometry", prhs[2]); + if (!xml) + return; + cfg.self = xml->getRootNode(); + CVolumeGeometry3D* pGeometry = new CVolumeGeometry3D(); + if (!pGeometry->initialize(cfg)) { + mexErrMsgTxt("Geometry class not initialized. \n"); + delete pGeometry; + delete xml; + return; + } + delete xml; + + // If data is specified, check dimensions + if (nrhs >= 4 && !mex_is_scalar(prhs[3])) { + get3DMatrixDims(prhs[3], dims); + if (pGeometry->getGridColCount() != dims[0] || pGeometry->getGridRowCount() != dims[1] || pGeometry->getGridSliceCount() != dims[2]) { + mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); + delete pGeometry; + return; + } + } + + // Initialize data object + pDataObject3D = new CFloat32VolumeData3DMemory(pGeometry); + delete pGeometry; + } + + else if (sDataType == "-sino" || sDataType == "-proj3d") { + + // Read geometry + if (!mxIsStruct(prhs[2])) { + mexErrMsgTxt("Argument 3 is not a valid MATLAB struct.\n"); + } + XMLDocument* xml = struct2XML("ProjectionGeometry", prhs[2]); + if (!xml) + return; + Config cfg; + cfg.self = xml->getRootNode(); + + // FIXME: Change how the base class is created. (This is duplicated + // in Projector2D.cpp.) + std::string type = cfg.self->getAttribute("type"); + CProjectionGeometry3D* pGeometry = 0; + if (type == "parallel3d") { + pGeometry = new CParallelProjectionGeometry3D(); + } else if (type == "parallel3d_vec") { + pGeometry = new CParallelVecProjectionGeometry3D(); + } else if (type == "cone") { + pGeometry = new CConeProjectionGeometry3D(); + } else if (type == "cone_vec") { + pGeometry = new CConeVecProjectionGeometry3D(); + } else { + mexErrMsgTxt("Invalid geometry type.\n"); + return; + } + + if (!pGeometry->initialize(cfg)) { + mexErrMsgTxt("Geometry class not initialized. \n"); + delete pGeometry; + delete xml; + return; + } + delete xml; + + // If data is specified, check dimensions + if (nrhs >= 4 && !mex_is_scalar(prhs[3])) { + get3DMatrixDims(prhs[3], dims); + if (pGeometry->getDetectorColCount() != dims[0] || pGeometry->getProjectionCount() != dims[1] || pGeometry->getDetectorRowCount() != dims[2]) { + mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); + delete pGeometry; + return; + } + } + + // Initialize data object + pDataObject3D = new CFloat32ProjectionData3DMemory(pGeometry); + } + + else if (sDataType == "-sinocone") { + // Read geometry + if (!mxIsStruct(prhs[2])) { + mexErrMsgTxt("Argument 3 is not a valid MATLAB struct.\n"); + } + XMLDocument* xml = struct2XML("ProjectionGeometry", prhs[2]); + if (!xml) + return; + Config cfg; + cfg.self = xml->getRootNode(); + CConeProjectionGeometry3D* pGeometry = new CConeProjectionGeometry3D(); + if (!pGeometry->initialize(cfg)) { + mexErrMsgTxt("Geometry class not initialized. \n"); + delete xml; + delete pGeometry; + return; + } + delete xml; + // If data is specified, check dimensions + if (nrhs >= 4 && !mex_is_scalar(prhs[3])) { + get3DMatrixDims(prhs[3], dims); + if (pGeometry->getDetectorRowCount() != dims[2] || pGeometry->getProjectionCount() != dims[1] || pGeometry->getDetectorColCount() != dims[0]) { + mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); + delete pGeometry; + return; + } + } + // Initialize data object + pDataObject3D = new CFloat32ProjectionData3DMemory(pGeometry); + delete pGeometry; + } + else { + mexErrMsgTxt("Invalid datatype. Please specify '-vol' or '-proj3d'. \n"); + return; + } + + // Check initialization + if (!pDataObject3D->isInitialized()) { + mexErrMsgTxt("Couldn't initialize data object.\n"); + delete pDataObject3D; + return; + } + + // Store data + + // fill with scalar value + if (nrhs < 4 || mex_is_scalar(prhs[3])) { + float32 fValue = 0.0f; + if (nrhs >= 4) + fValue = (float32)mxGetScalar(prhs[3]); + for (int i = 0; i < pDataObject3D->getSize(); ++i) { + pDataObject3D->getData()[i] = fValue; + } + } + // fill with array value + else if (mxIsDouble(prhs[3])) { + double* pdMatlabData = mxGetPr(prhs[3]); + int i = 0; + int col, row, slice; + for (slice = 0; slice < dims[2]; ++slice) { + for (row = 0; row < dims[1]; ++row) { + for (col = 0; col < dims[0]; ++col) { + // TODO: Benchmark and remove triple indexing? + pDataObject3D->getData3D()[slice][row][col] = pdMatlabData[i]; + ++i; + } + } + } + } + else if (mxIsSingle(prhs[3])) { + const float* pfMatlabData = (const float*)mxGetData(prhs[3]); + int i = 0; + int col, row, slice; + for (slice = 0; slice < dims[2]; ++slice) { + for (row = 0; row < dims[1]; ++row) { + for (col = 0; col < dims[0]; ++col) { + // TODO: Benchmark and remove triple indexing? + pDataObject3D->getData3D()[slice][row][col] = pfMatlabData[i]; + ++i; + } + } + } + } + pDataObject3D->updateStatistics(); + + // step4: store data object + int iIndex = CData3DManager::getSingleton().store(pDataObject3D); + + // step5: return data id + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(iIndex); + } + +} + +//----------------------------------------------------------------------------------------- +/** + * [id] = astra_mex_io_data('create_cache', config); + */ +void astra_mex_data3d_create_cache(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ +// if (nrhs < 2) { +// mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); +// return; +// } +// +// if (!mxIsStruct(prhs[1])) { +// mexErrMsgTxt("Argument 1 not a valid MATLAB struct. \n"); +// } +// +// // turn MATLAB struct to an XML-based Config object +// XMLDocument* xml = struct2XML("Data3D", prhs[1]); +// Config cfg; +// cfg.self = xml->getRootNode(); +// +// // create dataobject +// string sType = cfg.self->getAttribute("type"); +// int iIndex; +// if (sType == "ProjectionCached") { +// CFloat32ProjectionData3DCached* pData = new CFloat32ProjectionData3DCached(cfg); +// iIndex = CData3DManager::getSingleton().store(pData); +// } +//// else if (sType == "VolumeCached") { +//// CFloat32VolumeData3DCached* pData = new CFloat32VolumeData3DCached(cfg); +//// pData->initialize(cfg); +//// iIndex = CData3DManager::getSingleton().store(pData); +//// } +// +// // step4: set output +// if (1 <= nlhs) { +// plhs[0] = mxCreateDoubleScalar(iIndex); +// } + +} + + +//----------------------------------------------------------------------------------------- +/** + * data = astra_mex_data3d('get', id); + * + * Fetch data from the astra-library to a MATLAB matrix. + * id: identifier of the 3d data object as stored in the astra-library. + * data: MATLAB data + + */ +void astra_mex_data3d_get(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iDataID = (int)(mxGetScalar(prhs[1])); + + // step2: get data object + CFloat32Data3DMemory* pDataObject = dynamic_cast(astra::CData3DManager::getSingleton().get(iDataID)); + if (!pDataObject || !pDataObject->isInitialized()) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; + } + + // create output + if (1 <= nlhs) { + mwSize dims[3]; + dims[0] = pDataObject->getWidth(); + dims[1] = pDataObject->getHeight(); + dims[2] = pDataObject->getDepth(); + + plhs[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL); + double* out = mxGetPr(plhs[0]); + + int i = 0; + for (int slice = 0; slice < pDataObject->getDepth(); slice++) { + for (int row = 0; row < pDataObject->getHeight(); row++) { + for (int col = 0; col < pDataObject->getWidth(); col++) { + // TODO: Benchmark and remove triple indexing? + out[i] = pDataObject->getData3D()[slice][row][col]; + ++i; + } + } + } + } + +} + +//----------------------------------------------------------------------------------------- +/** + * data = astra_mex_data3d('get_single', id); + * + * Fetch data from the astra-library to a MATLAB matrix. + * id: identifier of the 3d data object as stored in the astra-library. + * data: MATLAB data + + */ +void astra_mex_data3d_get_single(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iDataID = (int)(mxGetScalar(prhs[1])); + + // step2: get data object + CFloat32Data3DMemory* pDataObject = dynamic_cast(astra::CData3DManager::getSingleton().get(iDataID)); + if (!pDataObject || !pDataObject->isInitialized()) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; + } + + // create output + if (1 <= nlhs) { + mwSize dims[3]; + dims[0] = pDataObject->getWidth(); + dims[1] = pDataObject->getHeight(); + dims[2] = pDataObject->getDepth(); + + plhs[0] = mxCreateNumericArray(3, dims, mxSINGLE_CLASS, mxREAL); + float* out = (float *)mxGetData(plhs[0]); + + int i = 0; + for (int slice = 0; slice < pDataObject->getDepth(); slice++) { + for (int row = 0; row < pDataObject->getHeight(); row++) { + for (int col = 0; col < pDataObject->getWidth(); col++) { + // TODO: Benchmark and remove triple indexing? + out[i] = pDataObject->getData3D()[slice][row][col]; + ++i; + } + } + } + } + +} + + +//----------------------------------------------------------------------------------------- +/** + * astra_mex_data3d('store', id, data); + * + * Store MATLAB matrix data in the astra-library. + * id: identifier of the 3d data object as stored in the astra-library. + * data: MATLAB data + + */ +void astra_mex_data3d_store(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: input + if (nrhs < 3) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iDataID = (int)(mxGetScalar(prhs[1])); + + // step2: get data object + CFloat32Data3DMemory* pDataObject = dynamic_cast(astra::CData3DManager::getSingleton().get(iDataID)); + if (!pDataObject || !pDataObject->isInitialized()) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; + } + + if (!(mex_is_scalar(prhs[2]) || mxIsDouble(prhs[2]) || mxIsSingle(prhs[2]))) { + mexErrMsgTxt("Data must be single or double."); + return; + } + + // fill with scalar value + if (mex_is_scalar(prhs[2])) { + float32 fValue = (float32)mxGetScalar(prhs[2]); + for (int i = 0; i < pDataObject->getSize(); ++i) { + pDataObject->getData()[i] = fValue; + } + } + // fill with array value + else if (mxIsDouble(prhs[2])) { + mwSize dims[3]; + get3DMatrixDims(prhs[2], dims); + if (dims[0] != pDataObject->getWidth() || dims[1] != pDataObject->getHeight() || dims[2] != pDataObject->getDepth()) { + mexErrMsgTxt("Data object dimensions don't match.\n"); + return; + + } + double* pdMatlabData = mxGetPr(prhs[2]); + int i = 0; + int col, row, slice; + for (slice = 0; slice < dims[2]; ++slice) { + for (row = 0; row < dims[1]; ++row) { + for (col = 0; col < dims[0]; ++col) { + // TODO: Benchmark and remove triple indexing? + pDataObject->getData3D()[slice][row][col] = pdMatlabData[i]; + ++i; + } + } + } + } + else if (mxIsSingle(prhs[2])) { + mwSize dims[3]; + get3DMatrixDims(prhs[2], dims); + if (dims[0] != pDataObject->getWidth() || dims[1] != pDataObject->getHeight() || dims[2] != pDataObject->getDepth()) { + mexErrMsgTxt("Data object dimensions don't match.\n"); + return; + + } + const float* pfMatlabData = (const float *)mxGetData(prhs[2]); + int i = 0; + int col, row, slice; + for (slice = 0; slice < dims[2]; ++slice) { + for (row = 0; row < dims[1]; ++row) { + for (col = 0; col < dims[0]; ++col) { + // TODO: Benchmark and remove triple indexing? + pDataObject->getData3D()[slice][row][col] = pfMatlabData[i]; + ++i; + } + } + } + } + pDataObject->updateStatistics(); +} + + +//----------------------------------------------------------------------------------------- +/** + * [id] = astra_mex_io_data('fetch_slice', id, slicenr); + */ +void astra_mex_data3d_fetch_slice_z(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ +// // step1: get input +// if (nrhs < 3) { +// mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); +// return; +// } +// int iDid = (int)(mxGetScalar(prhs[1])); +// int iSliceNr = (int)(mxGetScalar(prhs[2])); +// +// // Get data object +// CFloat32Data3D* pData = CData3DManager::getSingleton().get(iDid); +// if (!pData) { +// mexErrMsgTxt("DataObject not valid. \n"); +// return; +// } +// +// CFloat32Data2D* res = NULL; +// // Projection Data +// if (pData->getType() == CFloat32Data3D::PROJECTION) { +// CFloat32ProjectionData3D* pData2 = dynamic_cast(pData); +//// res = pData2->fetchSlice(iSliceNr); +// } +// // Volume Data +// else if (pData->getType() == CFloat32Data3D::VOLUME) { +// CFloat32VolumeData3D* pData2 = dynamic_cast(pData); +//// res = pData2->fetchSliceZ(iSliceNr); +// } +// // Error +// else { +// mexErrMsgTxt("DataObject not valid. \n"); +// return; +// } +// +// // store data +// int iIndex = CData2DManager::getSingleton().store(res); +// +// // step4: set output +// if (1 <= nlhs) { +// plhs[0] = mxCreateDoubleScalar(iIndex); +// } +} + +//----------------------------------------------------------------------------------------- +/** + * astra_mex_io_data('returnSlice', id, slicenr); + */ +void astra_mex_data3d_return_slice_z(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ +// // step1: get input +// if (nrhs < 3) { +// mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); +// return; +// } +// int iDid = (int)(mxGetScalar(prhs[1])); +// int iSliceNr = (int)(mxGetScalar(prhs[2])); +// +// // Get data object +// CFloat32Data3D* pData = CData3DManager::getSingleton().get(iDid); +// if (!pData) { +// mexErrMsgTxt("DataObject not valid. \n"); +// return; +// } +// +// // Projection Data +// if (pData->getType() == CFloat32Data3D::PROJECTION) { +// CFloat32ProjectionData3D* pData2 = dynamic_cast(pData); +//// TODO: think about returning slices +//// pData2->returnSlice(iSliceNr); +// } +// // Volume Data +// else if (pData->getType() == CFloat32Data3D::VOLUME) { +// CFloat32VolumeData3D* pData2 = dynamic_cast(pData); +//// TODO: think about returning slices +//// pData2->returnSliceZ(iSliceNr); +// } +// // Error +// else { +// mexErrMsgTxt("DataObject not valid. \n"); +// return; +// } +} + +//----------------------------------------------------------------------------------------- +/** + * [id] = astra_mex_io_data('fetch_projection', id, slicenr); + */ +void astra_mex_data3d_fetch_projection(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + //// step1: get input + //if (nrhs < 3) { + // mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + // return; + //} + //int iDid = (int)(mxGetScalar(prhs[1])); + //int iProjectionNr = (int)(mxGetScalar(prhs[2])); + + //// Get data object + //CFloat32Data3D* pData = CData3DManager::getSingleton().get(iDid); + //if (!pData) { + // mexErrMsgTxt("DataObject not valid. \n"); + // return; + //} + + //CFloat32Data2D* res = NULL; + //// Projection Data + //if (pData->getType() == CFloat32Data3D::PROJECTION) { + // CFloat32ProjectionData3D* pData2 = dynamic_cast(pData); + // res = pData2->fetchProjection(iProjectionNr); + //} + //// Error + //else { + // mexErrMsgTxt("DataObject not valid. \n"); + // return; + //} + // + //// store data + //int iIndex = CData2DManager::getSingleton().store(res); + + //// step4: set output + //if (1 <= nlhs) { + // plhs[0] = mxCreateDoubleScalar(iIndex); + //} +} + +//----------------------------------------------------------------------------------------- +/** + * astra_mex_io_data('return_projection', id, slicenr); + */ +void astra_mex_data3d_return_projection(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + //// step1: get input + //if (nrhs < 3) { + // mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + // return; + //} + //int iDid = (int)(mxGetScalar(prhs[1])); + //int iProjectionNr = (int)(mxGetScalar(prhs[2])); + + //// Get data object + //CFloat32Data3D* pData = CData3DManager::getSingleton().get(iDid); + //if (!pData) { + // mexErrMsgTxt("DataObject not valid. \n"); + // return; + //} + + //// Projection Data + //if (pData->getType() == CFloat32Data3D::PROJECTION) { + // CFloat32ProjectionData3D* pData2 = dynamic_cast(pData); + //// pData2->returnProjection(iProjectionNr); + //} + //// Error + //else { + // mexErrMsgTxt("DataObject not valid. \n"); + // return; + //} +} + +//----------------------------------------------------------------------------------------- +/** + * [id] = astra_mex_io_data('fetch_projection', id, slicenr); + */ +void astra_mex_data3d_fetch_slice_x(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + //// step1: get input + //if (nrhs < 3) { + // mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + // return; + //} + //int iDid = (int)(mxGetScalar(prhs[1])); + //int iSliceNr = (int)(mxGetScalar(prhs[2])); + + //// Get data object + //CFloat32Data3D* pData = CData3DManager::getSingleton().get(iDid); + //if (!pData) { + // mexErrMsgTxt("DataObject not valid. \n"); + // return; + //} + + //CFloat32Data2D* res = NULL; + //// Projection Data + //if (pData->getType() == CFloat32Data3D::VOLUME) { + // CFloat32VolumeData3D* pData2 = dynamic_cast(pData); + // res = pData2->fetchSliceX(iSliceNr); + //} + //// Error + //else { + // mexErrMsgTxt("DataObject not valid. \n"); + // return; + //} + // + //// store data + //int iIndex = CData2DManager::getSingleton().store(res); + + //// step4: set output + //if (1 <= nlhs) { + // plhs[0] = mxCreateDoubleScalar(iIndex); + //} +} + +//----------------------------------------------------------------------------------------- +/** + * astra_mex_io_data('return_slice_x', id, slicenr); + */ +void astra_mex_data3d_return_slice_x(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ +// // step1: get input +// if (nrhs < 3) { +// mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); +// return; +// } +// int iDid = (int)(mxGetScalar(prhs[1])); +// int iSliceNr = (int)(mxGetScalar(prhs[2])); +// +// // Get data object +// CFloat32Data3D* pData = CData3DManager::getSingleton().get(iDid); +// if (!pData) { +// mexErrMsgTxt("DataObject not valid. \n"); +// return; +// } +// +// // Projection Data +// if (pData->getType() == CFloat32Data3D::VOLUME) { +// CFloat32VolumeData3D* pData2 = dynamic_cast(pData); +//// TODO: think about returning slices +//// pData2->returnSliceX(iSliceNr); +// } +// // Error +// else { +// mexErrMsgTxt("DataObject not valid. \n"); +// return; +// } +} + + +//----------------------------------------------------------------------------------------- +/** + * [id] = astra_mex_io_data('fetch_slice_y', id, slicenr); + */ +void astra_mex_data3d_fetch_slice_y(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + //// step1: get input + //if (nrhs < 3) { + // mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + // return; + //} + //int iDid = (int)(mxGetScalar(prhs[1])); + //int iSliceNr = (int)(mxGetScalar(prhs[2])); + + //// Get data object + //CFloat32Data3D* pData = CData3DManager::getSingleton().get(iDid); + //if (!pData) { + // mexErrMsgTxt("DataObject not valid. \n"); + // return; + //} + + //CFloat32Data2D* res = NULL; + //// Projection Data + //if (pData->getType() == CFloat32Data3D::VOLUME) { + // CFloat32VolumeData3D* pData2 = dynamic_cast(pData); + // res = pData2->fetchSliceY(iSliceNr); + //} + //// Error + //else { + // mexErrMsgTxt("DataObject not valid. \n"); + // return; + //} + // + //// store data + //int iIndex = CData2DManager::getSingleton().store(res); + + //// step4: set output + //if (1 <= nlhs) { + // plhs[0] = mxCreateDoubleScalar(iIndex); + //} +} + +//----------------------------------------------------------------------------------------- +/** + * astra_mex_io_data('return_slice_y', id, slicenr); + */ +void astra_mex_data3d_return_slice_y(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ +// // step1: get input +// if (nrhs < 3) { +// mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); +// return; +// } +// int iDid = (int)(mxGetScalar(prhs[1])); +// int iSliceNr = (int)(mxGetScalar(prhs[2])); +// +// // Get data object +// CFloat32Data3D* pData = CData3DManager::getSingleton().get(iDid); +// if (!pData) { +// mexErrMsgTxt("DataObject not valid. \n"); +// return; +// } +// +// // Projection Data +// if (pData->getType() == CFloat32Data3D::VOLUME) { +// CFloat32VolumeData3D* pData2 = dynamic_cast(pData); +//// TODO: think about returning slices +//// pData2->returnSliceY(iSliceNr); +// } +// // Error +// else { +// mexErrMsgTxt("DataObject not valid. \n"); +// return; +// } +} + +//----------------------------------------------------------------------------------------- +/** + * [dim_x dim_y dim_z] = astra_mex_io_data('dimensions', id); + */ +void astra_mex_data3d_dimensions(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: get input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iDid = (int)(mxGetScalar(prhs[1])); + + // step2: get data object + CFloat32Data3D* pData; + if (!(pData = CData3DManager::getSingleton().get(iDid))) { + mexErrMsgTxt("DataObject not valid. \n"); + return; + } + + // step3: output + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(pData->getWidth()); + } + if (2 <= nlhs) { + plhs[1] = mxCreateDoubleScalar(pData->getHeight()); + } + if (3 <= nlhs) { + plhs[2] = mxCreateDoubleScalar(pData->getDepth()); + } +} + +//----------------------------------------------------------------------------------------- +/** + * [geom] = astra_mex_data3d('geometry', id); + */ +void astra_mex_data3d_geometry(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + //// Get input + //if (nrhs < 2) { + // mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + // return; + //} + //int iDid = (int)(mxGetScalar(prhs[1])); + + //// Get data object + //CFloat32Data3D* pData = CData3DManager::getSingleton().get(iDid); + //if (!pData) { + // mexErrMsgTxt("DataObject not valid. \n"); + // return; + //} + + //// Projection Data + //if (pData->getType() == CFloat32Data3D::PROJECTION) { + // CFloat32ProjectionData3D* pData2 = dynamic_cast(pData); + // CProjectionGeometry3D* pProjGeom = pData2->getGeometry(); + // XMLDocument* config = pProjGeom->toXML(); + + // if (1 <= nlhs) { + // plhs[0] = XML2struct(config); + // } + //} + //// Volume Data + //else if (pData->getType() == CFloat32Data3D::VOLUME) { + //// CFloat32VolumeData3D* pData2 = dynamic_cast(pData); + //// CVolumeGeometry2D* pVolGeom = pData2->getGeometry2D(iSliceNr); + //// if (1 <= nlhs) { + //// plhs[0] = createVolumeGeometryStruct(pVolGeom); + //// } + //} + //// Error + //else { + // mexErrMsgTxt("Type not valid. \n"); + // return; + //} +} + +//----------------------------------------------------------------------------------------- +/** + * [geom_xml] = astra_mex_data3d('geometry_xml', id); + */ +void astra_mex_data3d_geometry_xml(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + //// Get input + //if (nrhs < 2) { + // mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + // return; + //} + //int iDid = (int)(mxGetScalar(prhs[1])); + + //// Get data object + //CFloat32Data3D* pData = CData3DManager::getSingleton().get(iDid); + //if (!pData) { + // mexErrMsgTxt("DataObject not valid. \n"); + // return; + //} + + //// Projection Data + //if (pData->getType() == CFloat32Data3D::PROJECTION) { + // CFloat32ProjectionData3D* pData2 = dynamic_cast(pData); + // CProjectionGeometry3D* pProjGeom = pData2->getGeometry(); + // XMLDocument* config = pProjGeom->toXML(); + + // if (1 <= nlhs) { + // plhs[0] = mxCreateString(config->getRootNode()->toString().c_str()); + // } + //} + //// Volume Data + //else if (pData->getType() == CFloat32Data3D::VOLUME) { + //// CFloat32VolumeData3D* pData2 = dynamic_cast(pData); + //// CVolumeGeometry2D* pVolGeom = pData2->getGeometry2D(iSliceNr); + //// if (1 <= nlhs) { + //// plhs[0] = createVolumeGeometryStruct(pVolGeom); + //// } + //} + //// Error + //else { + // mexErrMsgTxt("Type not valid. \n"); + // return; + //} +} +//----------------------------------------------------------------------------------------- +/** + * astra_mex_data3d('delete', did1, did2, ...); + */ +void astra_mex_data3d_delete(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: read input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + + for (int i = 1; i < nrhs; i++) { + int iDataID = (int)(mxGetScalar(prhs[i])); + CData3DManager::getSingleton().remove(iDataID); + } +} + +//----------------------------------------------------------------------------------------- +/** + * astra_mex_data3d('clear'); + */ +void astra_mex_data3d_clear(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + CData3DManager::getSingleton().clear(); +} + +//----------------------------------------------------------------------------------------- +/** + * astra_mex_data3d('info'); + */ +void astra_mex_data3d_info(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + mexPrintf("%s", astra::CData3DManager::getSingleton().info().c_str()); +} + +//----------------------------------------------------------------------------------------- + +static void printHelp() +{ + mexPrintf("Please specify a mode of operation.\n"); + mexPrintf("Valid modes: create, create_cache, get, get_single, delete, clear, info\n"); + mexPrintf(" fetch_projection, return_projection, fetch_slice[_z],\n"); + mexPrintf(" return_slice[_z], fetch_slice_x, return slice_x\n"); + mexPrintf(" fetch_slice_y, return slice_y, dimensions, geometry\n"); + mexPrintf(" geometry_xml\n"); +} + + +//----------------------------------------------------------------------------------------- +/** + * ... = astra_mex_io_data(mode,...); + */ +void mexFunction(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + + // INPUT: Mode + string sMode = ""; + if (1 <= nrhs) { + sMode = mex_util_get_string(prhs[0]); + } else { + printHelp(); + return; + } + + // 3D data + if (sMode == std::string("create")) { + astra_mex_data3d_create(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("create_cache")) { + astra_mex_data3d_create_cache(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("get")) { + astra_mex_data3d_get(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("get_single")) { + astra_mex_data3d_get_single(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("store") || + sMode == std::string("set")) { + astra_mex_data3d_store(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("delete")) { + astra_mex_data3d_delete(nlhs, plhs, nrhs, prhs); + } else if (sMode == "clear") { + astra_mex_data3d_clear(nlhs, plhs, nrhs, prhs); + } else if (sMode == "info") { + astra_mex_data3d_info(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("fetch_projection")) { + astra_mex_data3d_fetch_projection(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("return_projection")) { + astra_mex_data3d_return_projection(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("fetch_slice") || sMode == std::string("fetch_slice_z")) { + astra_mex_data3d_fetch_slice_z(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("return_slice") || sMode == std::string("return_slice_z")) { + astra_mex_data3d_return_slice_z(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("fetch_slice_x")) { + astra_mex_data3d_fetch_slice_x(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("return_slice_x")) { + astra_mex_data3d_return_slice_x(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("fetch_slice_y")) { + astra_mex_data3d_fetch_slice_y(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("return_slice_y")) { + astra_mex_data3d_return_slice_y(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("dimensions")) { + astra_mex_data3d_dimensions(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("geometry")) { + astra_mex_data3d_geometry(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("geometry_xml")) { + astra_mex_data3d_geometry_xml(nlhs, plhs, nrhs, prhs); + } else { + printHelp(); + } + + return; +} + + diff --git a/matlab/mex/astra_mex_data3d_vc08.vcproj b/matlab/mex/astra_mex_data3d_vc08.vcproj new file mode 100644 index 0000000..2e69c16 --- /dev/null +++ b/matlab/mex/astra_mex_data3d_vc08.vcproj @@ -0,0 +1,588 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/matlab/mex/astra_mex_matrix_c.cpp b/matlab/mex/astra_mex_matrix_c.cpp new file mode 100644 index 0000000..accaab5 --- /dev/null +++ b/matlab/mex/astra_mex_matrix_c.cpp @@ -0,0 +1,437 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +$Id$ +*/ + +/** \file astra_mex_matrix_c.cpp + * + * \brief Create sparse (projection) matrices in the ASTRA workspace + */ +#include +#include "mexHelpFunctions.h" + +#include + +#include "astra/Globals.h" + +#include "astra/AstraObjectManager.h" + +#include "astra/SparseMatrix.h" + +using namespace std; +using namespace astra; + +//----------------------------------------------------------------------------------------- +/** astra_mex_matrix('delete', id1, id2, ...); + * + * Delete one or more data objects currently stored in the astra-library. + * id1, id2, ... : identifiers of the 2d data objects as stored in the astra-library. + */ +void astra_mex_matrix_delete(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: read input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + + // step2: delete all specified data objects + for (int i = 1; i < nrhs; i++) { + int iDataID = (int)(mxGetScalar(prhs[i])); + CMatrixManager::getSingleton().remove(iDataID); + } +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_matrix('clear'); + * + * Delete all data objects currently stored in the astra-library. + */ +void astra_mex_matrix_clear(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + CMatrixManager::getSingleton().clear(); +} + + + +static bool matlab_to_astra(const mxArray* _rhs, CSparseMatrix* _pMatrix) +{ + // Check input + if (!mxIsSparse (_rhs)) { + mexErrMsgTxt("Argument is not a valid MATLAB sparse matrix.\n"); + return false; + } + if (!_pMatrix->isInitialized()) { + mexErrMsgTxt("Couldn't initialize data object.\n"); + return false; + } + + unsigned int iHeight = mxGetM(_rhs); + unsigned int iWidth = mxGetN(_rhs); + unsigned long lSize = mxGetNzmax(_rhs); + + if (_pMatrix->m_lSize < lSize || _pMatrix->m_iHeight < iHeight) { + // TODO: support resizing? + mexErrMsgTxt("Matrix too large to store in this object.\n"); + return false; + } + + // Transpose matrix, as matlab stores a matrix column-by-column + // but we want it row-by-row. + // 1. Compute sizes of rows. We store these in _pMatrix->m_plRowStarts. + // 2. Fill data structure + // Complexity: O( #rows + #entries ) + + for (unsigned int i = 0; i <= iHeight; ++i) + _pMatrix->m_plRowStarts[i] = 0; + + mwIndex *colStarts = mxGetJc(_rhs); + mwIndex *rowIndices = mxGetIr(_rhs); + double *floatValues = 0; + bool *boolValues = 0; + bool bLogical = mxIsLogical(_rhs); + if (bLogical) + boolValues = mxGetLogicals(_rhs); + else + floatValues = mxGetPr(_rhs); + + for (mwIndex i = 0; i < colStarts[iWidth]; ++i) { + unsigned int iRow = rowIndices[i]; + assert(iRow < iHeight); + _pMatrix->m_plRowStarts[iRow+1]++; + } + + // Now _pMatrix->m_plRowStarts[i+1] is the number of entries in row i + + for (unsigned int i = 1; i <= iHeight; ++i) + _pMatrix->m_plRowStarts[i] += _pMatrix->m_plRowStarts[i-1]; + + // Now _pMatrix->m_plRowStarts[i+1] is the number of entries in rows <= i, + // so the intended start of row i+1 + + int iCol = 0; + for (mwIndex i = 0; i < colStarts[iWidth]; ++i) { + while (i >= colStarts[iCol+1]) + iCol++; + + unsigned int iRow = rowIndices[i]; + assert(iRow < iHeight); + float32 fVal; + if (bLogical) + fVal = (float32)boolValues[i]; + else + fVal = (float32)floatValues[i]; + + unsigned long lIndex = _pMatrix->m_plRowStarts[iRow]++; + _pMatrix->m_pfValues[lIndex] = fVal; + _pMatrix->m_piColIndices[lIndex] = iCol; + } + + // Now _pMatrix->m_plRowStarts[i] is the start of row i+1 + + for (unsigned int i = iHeight; i > 0; --i) + _pMatrix->m_plRowStarts[i] = _pMatrix->m_plRowStarts[i-1]; + _pMatrix->m_plRowStarts[0] = 0; + +#if 0 + // Debugging: dump matrix + for (unsigned int i = 0; i < iHeight; ++i) { + printf("Row %d: %ld-%ld\n", i, _pMatrix->m_plRowStarts[i], _pMatrix->m_plRowStarts[i+1]); + for (unsigned long j = _pMatrix->m_plRowStarts[i]; j < _pMatrix->m_plRowStarts[i+1]; ++j) { + printf("(%d,%d) = %f\n", i, _pMatrix->m_piColIndices[j], _pMatrix->m_pfValues[j]); + } + } +#endif + + return true; +} + +static bool astra_to_matlab(const CSparseMatrix* _pMatrix, mxArray*& _lhs) +{ + if (!_pMatrix->isInitialized()) { + mexErrMsgTxt("Uninitialized data object.\n"); + return false; + } + + unsigned int iHeight = _pMatrix->m_iHeight; + unsigned int iWidth = _pMatrix->m_iWidth; + unsigned long lSize = _pMatrix->m_lSize; + + _lhs = mxCreateSparse(iHeight, iWidth, lSize, mxREAL); + if (!mxIsSparse (_lhs)) { + mexErrMsgTxt("Couldn't initialize matlab sparse matrix.\n"); + return false; + } + + mwIndex *colStarts = mxGetJc(_lhs); + mwIndex *rowIndices = mxGetIr(_lhs); + double *floatValues = mxGetPr(_lhs); + + for (unsigned int i = 0; i <= iWidth; ++i) + colStarts[i] = 0; + + for (unsigned int i = 0; i < _pMatrix->m_plRowStarts[iHeight]; ++i) { + unsigned int iCol = _pMatrix->m_piColIndices[i]; + assert(iCol < iWidth); + colStarts[iCol+1]++; + } + // Now _pMatrix->m_plRowStarts[i+1] is the number of entries in row i + + for (unsigned int i = 1; i <= iWidth; ++i) + colStarts[i] += colStarts[i-1]; + // Now _pMatrix->m_plRowStarts[i+1] is the number of entries in rows <= i, + // so the intended start of row i+1 + + unsigned int iRow = 0; + for (unsigned int i = 0; i < _pMatrix->m_plRowStarts[iHeight]; ++i) { + while (i >= _pMatrix->m_plRowStarts[iRow+1]) + iRow++; + + unsigned int iCol = _pMatrix->m_piColIndices[i]; + assert(iCol < iWidth); + double fVal = _pMatrix->m_pfValues[i]; + unsigned long lIndex = colStarts[iCol]++; + floatValues[lIndex] = fVal; + rowIndices[lIndex] = iRow; + } + // Now _pMatrix->m_plRowStarts[i] is the start of row i+1 + + for (unsigned int i = iWidth; i > 0; --i) + colStarts[i] = colStarts[i-1]; + colStarts[0] = 0; + + return true; +} + +//----------------------------------------------------------------------------------------- +/** id = astra_mex_matrix('create', data); + * + * Create a new matrix object in the astra-library. + * data: a sparse MATLAB matrix containing the data. + * id: identifier of the matrix object as it is now stored in the astra-library. + */ +void astra_mex_matrix_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArray* prhs[]) +{ + // step1: get datatype + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + + if (!mxIsSparse (prhs[1])) { + mexErrMsgTxt("Argument is not a valid MATLAB sparse matrix.\n"); + return; + } + + unsigned int iHeight = mxGetM(prhs[1]); + unsigned int iWidth = mxGetN(prhs[1]); + unsigned long lSize = mxGetNzmax(prhs[1]); + + CSparseMatrix* pMatrix = new CSparseMatrix(iHeight, iWidth, lSize); + + // Check initialization + if (!pMatrix->isInitialized()) { + mexErrMsgTxt("Couldn't initialize data object.\n"); + delete pMatrix; + return; + } + + bool bResult = matlab_to_astra(prhs[1], pMatrix); + + if (!bResult) { + mexErrMsgTxt("Failed to create data object.\n"); + delete pMatrix; + return; + } + + // store data object + int iIndex = CMatrixManager::getSingleton().store(pMatrix); + + // return data id + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(iIndex); + } +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_matrix('store', id, data); + * + * Store a sparse MATLAB matrix in an existing astra matrix dataobject. + * id: identifier of the 2d data object as stored in the astra-library. + * data: a sparse MATLAB matrix. + */ +void astra_mex_matrix_store(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: input + if (nrhs < 3) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + if (!mxIsDouble(prhs[1])) { + mexErrMsgTxt("Identifier should be a scalar value. \n"); + return; + } + int iDataID = (int)(mxGetScalar(prhs[1])); + + // step2: get data object + CSparseMatrix* pMatrix = astra::CMatrixManager::getSingleton().get(iDataID); + if (!pMatrix || !pMatrix->isInitialized()) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; + } + + bool bResult = matlab_to_astra(prhs[2], pMatrix); + if (!bResult) { + mexErrMsgTxt("Failed to store matrix.\n"); + } +} + +//----------------------------------------------------------------------------------------- +/** geom = astra_mex_matrix('get_size', id); + * + * Fetch the dimensions and size of a matrix stored in the astra-library. + * id: identifier of the 2d data object as stored in the astra-library. + * geom: a 1x2 matrix containing [rows, columns] + */ +void astra_mex_matrix_get_size(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + if (!mxIsDouble(prhs[1])) { + mexErrMsgTxt("Identifier should be a scalar value. \n"); + return; + } + int iDataID = (int)(mxGetScalar(prhs[1])); + + // step2: get data object + CSparseMatrix* pMatrix = astra::CMatrixManager::getSingleton().get(iDataID); + if (!pMatrix || !pMatrix->isInitialized()) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; + } + + // create output + // TODO +} + +//----------------------------------------------------------------------------------------- +/** data = astra_mex_matrix('get', id); + * + * Fetch data from the astra-library to a MATLAB matrix. + * id: identifier of the matrix data object as stored in the astra-library. + * data: MATLAB + */ +void astra_mex_matrix_get(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: check input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + if (!mxIsDouble(prhs[1])) { + mexErrMsgTxt("Identifier should be a scalar value. \n"); + return; + } + int iDataID = (int)(mxGetScalar(prhs[1])); + + // step2: get data object + CSparseMatrix* pMatrix = astra::CMatrixManager::getSingleton().get(iDataID); + if (!pMatrix || !pMatrix->isInitialized()) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; + } + + // create output + if (1 <= nlhs) { + bool bResult = astra_to_matlab(pMatrix, plhs[0]); + if (!bResult) { + mexErrMsgTxt("Failed to get matrix.\n"); + } + } + +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_matrix('info'); + * + * Print information about all the matrix objects currently stored in the astra-library. + */ +void astra_mex_matrix_info(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + mexPrintf("%s", astra::CMatrixManager::getSingleton().info().c_str()); +} + +//----------------------------------------------------------------------------------------- + +static void printHelp() +{ + mexPrintf("Please specify a mode of operation.\n"); + mexPrintf("Valid modes: get, delete, clear, store, create, get_size, info\n"); +} + +//----------------------------------------------------------------------------------------- +/** + * ... = astra_mex_matrix(type,...); + */ +void mexFunction(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + + // INPUT0: Mode + string sMode = ""; + if (1 <= nrhs) { + sMode = mex_util_get_string(prhs[0]); + } else { + printHelp(); + return; + } + + // SWITCH (MODE) + if (sMode == std::string("get")) { + astra_mex_matrix_get(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("delete")) { + astra_mex_matrix_delete(nlhs, plhs, nrhs, prhs); + } else if (sMode == "clear") { + astra_mex_matrix_clear(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("store")) { + astra_mex_matrix_store(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("create")) { + astra_mex_matrix_create(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("get_size")) { + astra_mex_matrix_get_size(nlhs, plhs, nrhs, prhs); + } else if (sMode == std::string("info")) { + astra_mex_matrix_info(nlhs, plhs, nrhs, prhs); + } else { + printHelp(); + } + + return; +} + + diff --git a/matlab/mex/astra_mex_matrix_vc08.vcproj b/matlab/mex/astra_mex_matrix_vc08.vcproj new file mode 100644 index 0000000..47509f6 --- /dev/null +++ b/matlab/mex/astra_mex_matrix_vc08.vcproj @@ -0,0 +1,591 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/matlab/mex/astra_mex_projector3d_c.cpp b/matlab/mex/astra_mex_projector3d_c.cpp new file mode 100644 index 0000000..1385863 --- /dev/null +++ b/matlab/mex/astra_mex_projector3d_c.cpp @@ -0,0 +1,433 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +$Id$ +*/ + +/** \file astra_mex_projector3d_c.cpp + * + * \brief Create and manage 3d projectors in the ASTRA workspace + */ + +#include +#include "mexHelpFunctions.h" + +#include "astra/Globals.h" + +#include "astra/Projector3D.h" +#include "astra/AstraObjectManager.h" +#include "astra/AstraObjectFactory.h" + +#include "astra/ProjectionGeometry3D.h" +#include "astra/VolumeGeometry3D.h" + +#include +#include + +using namespace std; +using namespace astra; + +//----------------------------------------------------------------------------------------- +/** +* [pid] = astra_mex_projector('create', cfgstruct); +*/ +void astra_mex_projector3d_create(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + + if (!mxIsStruct(prhs[1])) { + mexErrMsgTxt("Argument 1 not a valid MATLAB struct. \n"); + } + + // turn MATLAB struct to an XML-based Config object + XMLDocument* xml = struct2XML("Projector3D", prhs[1]); + Config cfg; + cfg.self = xml->getRootNode(); + + // create algorithm + CProjector3D* pProj = CProjector3DFactory::getSingleton().create(cfg); + if (pProj == NULL) { + mexErrMsgTxt("Error creating Projector3D. \n"); + return; + } + + // store projector + int iIndex = CProjector3DManager::getSingleton().store(pProj); + + // step4: set output + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(iIndex); + } + + } + +//----------------------------------------------------------------------------------------- +/** +* astra_mex_projector3d('destroy', pid1, pid2, ...); +*/ +void astra_mex_projector3d_destroy(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: read input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + + for (int i = 1; i < nrhs; i++) { + int iPid = (int)(mxGetScalar(prhs[i])); + CProjector3DManager::getSingleton().remove(iPid); + } +} + +//----------------------------------------------------------------------------------------- +/** + * astra_mex_projector3d('clear'); + */ +void astra_mex_projector3d_clear(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + CProjector3DManager::getSingleton().clear(); +} + + +//----------------------------------------------------------------------------------------- +/** +* [proj_geom] = astra_mex_projector3d('get_projection_geometry', pid); +*/ +void astra_mex_projector3d_get_projection_geometry(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: read input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iPid = (int)(mxGetScalar(prhs[1])); + + // step2: get projector + CProjector3D* pProjector; + if (!(pProjector = CProjector3DManager::getSingleton().get(iPid))) { + mexErrMsgTxt("Projector not found.\n"); + return; + } + + // step3: get projection_geometry and turn it into a MATLAB struct + //if (1 <= nlhs) { + // plhs[0] = createProjectionGeometryStruct(pProjector->getProjectionGeometry()); + //} +} + +//----------------------------------------------------------------------------------------- +/** +* [recon_geom] = astra_mex_projector3d('get_volume_geometry', pid); +*/ +void astra_mex_projector3d_get_reconstruction_geometry(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: read input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iPid = (int)(mxGetScalar(prhs[1])); + + // step2: get projector + CProjector3D* pProjector; + if (!(pProjector = CProjector3DManager::getSingleton().get(iPid))) { + mexErrMsgTxt("Projector not found.\n"); + return; + } + + // step3: get projection_geometry and turn it into a MATLAB struct + //if (1 <= nlhs) { + // plhs[0] = createVolumeGeometryStruct(pProjector->getVolumeGeometry()); + //} +} + +//----------------------------------------------------------------------------------------- +/** +* [weights] = astra_mex_projector3d('weights_single_ray', pid, projection_index, detector_index); +*/ +void astra_mex_projector_weights_single_ray(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + //// step1: get input + //if (nrhs < 4) { + // mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + // return; + //} + //int iPid = (int)(mxGetScalar(prhs[1])); + //int iProjectionIndex = (int)(mxGetScalar(prhs[2])); + //int iDetectorIndex = (int)(mxGetScalar(prhs[3])); + + //// step2: get projector + //CProjector3D* pProjector; + //if (!(pProjector = CProjector3DManager::getSingleton().get(iPid))) { + // mexErrMsgTxt("Projector not found.\n"); + // return; + //} + // + //// step3: create output vars + //int iStoredPixelCount; + //int iMaxPixelCount = pProjector->getProjectionWeightsCount(iProjectionIndex); + //SWeightedPixel* pPixelsWeights = new SWeightedPixel3D[iMaxPixelCount]; + // + //// step4: perform operation + //pProjector->computeSingleRayWeights(iProjectionIndex, + // iDetectorIndex, + // pPixelsWeights, + // iMaxPixelCount, + // iStoredPixelCount); + + //// step5: return output + //if (1 <= nlhs) { + // mwSize dims[2]; + // dims[0] = iStoredPixelCount; + // dims[1] = 2; + + // plhs[0] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL); + // double* out = mxGetPr(plhs[0]); + + // for (int col = 0; col < iStoredPixelCount; col++) { + // out[col] = pPixelsWeights[col].m_iIndex; + // out[iStoredPixelCount+col] = pPixelsWeights[col].m_fWeight; + // //cout << pPixelsWeights[col].m_iIndex << " " << pPixelsWeights[col].m_fWeight <getProjectionWeightsCount(iProjectionIndex)]; +// int* piRayStoredPixelCount = new int[pProjector->getProjectionGeometry()->getDetectorCount()]; +// +// // step4: perform operation +// pProjector->computeProjectionRayWeights(iProjectionIndex, pPixelsWheights, piRayStoredPixelCount); +// +// // step5: return output +// if (1 <= nlhs) { +// // get basic values +// int iMatrixSize = pProjector->getVolumeGeometry()->getWindowLengthX() * +// pProjector->getVolumeGeometry()->getWindowLengthY(); +// int iDetectorCount = pProjector->getProjectionGeometry()->getDetectorCount(); +// int iTotalStoredPixelCount = 0; +// for (int i = 0; i < iDetectorCount; i++) { +// iTotalStoredPixelCount += piRayStoredPixelCount[i]; +// } +// +// // create matlab sparse matrix +// plhs[0] = mxCreateSparse(iMatrixSize, // number of rows (#pixels) +// iDetectorCount, // number of columns (#detectors) +// iTotalStoredPixelCount, // number of non-zero elements +// mxREAL); // element type +// double* values = mxGetPr(plhs[0]); +// mwIndex* rows = mxGetIr(plhs[0]); +// mwIndex* cols = mxGetJc(plhs[0]); +// +// int currentBase = 0; +// int currentIndex = 0; +// for (int i = 0; i < iDetectorCount; i++) { +// for (int j = 0; j < piRayStoredPixelCount[i]; j++) { +// values[currentIndex + j] = pPixelsWheights[currentBase + j].m_fWeight; +// rows[currentIndex + j] = pPixelsWheights[currentBase + j].m_iIndex; +// } +// +// currentBase += pProjector->getProjectionWeightsCount(iProjectionIndex) / pProjector->getProjectionGeometry()->getDetectorCount(); +// currentIndex += piRayStoredPixelCount[i]; +// } +// cols[0] = piRayStoredPixelCount[0]; +// for (int j = 1; j < iDetectorCount; j++) { +// cols[j] = cols[j-1] + piRayStoredPixelCount[j]; +// } +// cols[iDetectorCount] = iTotalStoredPixelCount; +// } +// +//} +// +////----------------------------------------------------------------------------------------- +///** +//* output = astra_mex_projector('splat', pid, x, y); +//*/ +//void astra_mex_projector_splat(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +//{ +// // step1: get input +// if (nrhs < 4) { +// mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); +// return; +// } +// int iPid = (int)(mxGetScalar(prhs[1])); +// int iX = (int)(mxGetScalar(prhs[2])); +// int iY = (int)(mxGetScalar(prhs[3])); +// +// // step2: get projector +// CProjector2D* pProjector; +// if (!(pProjector = CProjectorManager::getSingleton().get(iPid))) { +// mexErrMsgTxt("Projector not found.\n"); +// return; +// } +// +// // step3: perform action +// vector detinfo = pProjector->projectPoint(iX, iY); +// +// // step4: output +// if (nlhs <= 1) { +// plhs[0] = mxCreateDoubleMatrix(detinfo.size(), // # rows +// 2, // # cols +// mxREAL); // datatype 32-bits +// double* out = mxGetPr(plhs[0]); +// +// // fill up output +// int i = 0; +// for (int x = 0; x < detinfo.size() ; x++) { +// out[i] = detinfo[x].m_iAngleIndex; +// i++; +// } +// for (int x = 0; x < detinfo.size() ; x++) { +// out[i] = detinfo[x].m_iDetectorIndex; +// i++; +// } +// } +// +// +//} + +//----------------------------------------------------------------------------------------- +/** result = astra_mex_projector3d('is_cuda', id); + * + * Return is the specified projector is a cuda projector. + * id: identifier of the projector object as stored in the astra-library. + */ +void astra_mex_projector3d_is_cuda(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: get input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iPid = (int)(mxGetScalar(prhs[1])); + + // step2: get projector + CProjector3D* pProjector = CProjector3DManager::getSingleton().get(iPid); + if (!pProjector || !pProjector->isInitialized()) { + mexErrMsgTxt("Projector not initialized.\n"); + return; + } + +#ifdef ASTRA_CUDA + CCudaProjector3D* pCP = dynamic_cast(pProjector); + plhs[0] = mxCreateLogicalScalar(pCP ? 1 : 0); +#else + plhs[0] = mxCreateLogicalScalar(0); +#endif +} + + +//----------------------------------------------------------------------------------------- +/** + * astra_mex_projector3d('help'); + */ +void astra_mex_projector3d_help(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + cout << "astra_mex_projector3d help:" < + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/matlab/mex/astra_mex_projector_c.cpp b/matlab/mex/astra_mex_projector_c.cpp new file mode 100644 index 0000000..5cbe502 --- /dev/null +++ b/matlab/mex/astra_mex_projector_c.cpp @@ -0,0 +1,510 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +$Id$ +*/ + +/** \file astra_mex_projector_c.cpp + * + * \brief Create and manage 2d projectors in the ASTRA workspace + */ +#include "astra/Globals.h" + +#include +#include "mexHelpFunctions.h" + +#include "astra/AstraObjectManager.h" +#include "astra/Projector2D.h" +#include "astra/AstraObjectFactory.h" + +#include "astra/Float32VolumeData2D.h" + +#include "astra/ProjectionGeometry2D.h" +#include "astra/ParallelProjectionGeometry2D.h" +#include "astra/VolumeGeometry2D.h" + + +#include +#include + +using namespace std; +using namespace astra; + +//----------------------------------------------------------------------------------------- +/** id = astra_mex_projector('create', cfg); + * + * Create and configure a new projector object. + * cfg: MATLAB struct containing the configuration parameters, see doxygen documentation for details. + * id: identifier of the projector object as it is now stored in the astra-library. + */ +void astra_mex_projector_create(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + int iIndex = 0; + + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + + if (!mxIsStruct(prhs[1])) { + mexErrMsgTxt("Argument 1 not a valid MATLAB struct. \n"); + } + + + // turn MATLAB struct to an XML-based Config object + XMLDocument* xml = struct2XML("Projector2D", prhs[1]); + Config cfg; + cfg.self = xml->getRootNode(); + + // create algorithm + CProjector2D* pProj = CProjector2DFactory::getSingleton().create(cfg); + if (pProj == NULL) { + delete xml; + mexErrMsgTxt("Error creating projector. \n"); + return; + } + + delete xml; + + // store projector + iIndex = CProjector2DManager::getSingleton().store(pProj); + + // step4: set output + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(iIndex); + } +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_projector('delete', id1, id2, ...); + * + * Delete one or more projector objects currently stored in the astra-library. + * id1, id2, ... : identifiers of the projector objects as stored in the astra-library. + */ +void astra_mex_projector_delete(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: read input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + + for (int i = 1; i < nrhs; i++) { + int iPid = (int)(mxGetScalar(prhs[i])); + CProjector2DManager::getSingleton().remove(iPid); + } +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_projector('clear'); + * + * Delete all projector objects currently stored in the astra-library. + */ +void astra_mex_projector_clear(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + CProjector2DManager::getSingleton().clear(); +} + +//----------------------------------------------------------------------------------------- +/** astra_mex_projector('info'); + * + * Print information about all the projector objects currently stored in the astra-library. + */ +void astra_mex_projector_info(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + mexPrintf("%s", astra::CProjector2DManager::getSingleton().info().c_str()); +} + +//----------------------------------------------------------------------------------------- +/** proj_geom = astra_mex_projector('projection_geometry', id); + * + * Fetch the projection geometry of a certain projector. + * id: identifier of the projector object as stored in the astra-library. + * proj_geom: MATLAB struct containing all information about the projection geometry +*/ +void astra_mex_projector_projection_geometry(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: read input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iPid = (int)(mxGetScalar(prhs[1])); + + // step2: get projector + CProjector2D* pProjector = CProjector2DManager::getSingleton().get(iPid); + if (!pProjector || !pProjector->isInitialized()) { + mexErrMsgTxt("Projector not initialized.\n"); + return; + } + + // step3: get projection_geometry and turn it into a MATLAB struct + if (1 <= nlhs) { + plhs[0] = createProjectionGeometryStruct(pProjector->getProjectionGeometry()); + } +} + +//----------------------------------------------------------------------------------------- +/** vol_geom = astra_mex_projector('volume_geometry', id); + * + * Fetch the volume geometry of a certain projector. + * id: identifier of the projector object as stored in the astra-library. + * vol_geom: MATLAB struct containing all information about the volume geometry + */ +void astra_mex_projector_volume_geometry(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: read input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iPid = (int)(mxGetScalar(prhs[1])); + + // step2: get projector + CProjector2D* pProjector = CProjector2DManager::getSingleton().get(iPid); + if (!pProjector || !pProjector->isInitialized()) { + mexErrMsgTxt("Projector not initialized.\n"); + return; + } + + // step3: get projection_geometry and turn it into a MATLAB struct + if (1 <= nlhs) { + plhs[0] = createVolumeGeometryStruct(pProjector->getVolumeGeometry()); + } +} + +//----------------------------------------------------------------------------------------- +/** weights = astra_mex_projector('weights_single_ray', id, projection_index, detector_index); + * + * Calculate the nonzero weights of a certain projection ray. + * id: identifier of the projector object as stored in the astra-library. + * projection_index: index of the projection angle + * detector_index: index of the detector + * weights: list of computed weights [pixel_index, weight] + */ +void astra_mex_projector_weights_single_ray(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: get input + if (nrhs < 4) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iPid = (int)(mxGetScalar(prhs[1])); + int iProjectionIndex = (int)(mxGetScalar(prhs[2])); + int iDetectorIndex = (int)(mxGetScalar(prhs[3])); + + // step2: get projector + CProjector2D* pProjector = CProjector2DManager::getSingleton().get(iPid); + if (!pProjector || !pProjector->isInitialized()) { + mexErrMsgTxt("Projector not initialized.\n"); + return; + } + + // step3: create output vars + int iStoredPixelCount; + int iMaxPixelCount = pProjector->getProjectionWeightsCount(iProjectionIndex); + SPixelWeight* pPixelsWeights = new SPixelWeight[iMaxPixelCount]; + + // step4: perform operation + pProjector->computeSingleRayWeights(iProjectionIndex, + iDetectorIndex, + pPixelsWeights, + iMaxPixelCount, + iStoredPixelCount); + + // step5: return output + if (1 <= nlhs) { + mwSize dims[2]; + dims[0] = iStoredPixelCount; + dims[1] = 2; + + plhs[0] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL); + double* out = mxGetPr(plhs[0]); + + for (int col = 0; col < iStoredPixelCount; col++) { + out[col] = pPixelsWeights[col].m_iIndex; + out[iStoredPixelCount+col] = pPixelsWeights[col].m_fWeight; + } + } + + // garbage collection + delete[] pPixelsWeights; +} + +//----------------------------------------------------------------------------------------- +/** weights = astra_mex_projector('weights_projection', id, projection_index); + * + * Calculate the nonzero weights of all rays in a certain projection. + * id: identifier of the projector object as stored in the astra-library. + * projection_index: index of the projection angle + * weights: sparse matrix containing the rows of the projection matric belonging to the requested projection angle. + */ +void astra_mex_projector_weights_projection(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: get input + if (nrhs < 3) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iPid = (int)(mxGetScalar(prhs[1])); + int iProjectionIndex = (int)(mxGetScalar(prhs[2])); + + // step2: get projector + CProjector2D* pProjector = CProjector2DManager::getSingleton().get(iPid); + if (!pProjector || !pProjector->isInitialized()) { + mexErrMsgTxt("Projector not initialized.\n"); + return; + } + + // step3: create output vars + SPixelWeight* pPixelsWheights = new SPixelWeight[pProjector->getProjectionWeightsCount(iProjectionIndex)]; + int* piRayStoredPixelCount = new int[pProjector->getProjectionGeometry()->getDetectorCount()]; + + // step4: perform operation + pProjector->computeProjectionRayWeights(iProjectionIndex, pPixelsWheights, piRayStoredPixelCount); + + // step5: return output + if (1 <= nlhs) { + // get basic values + int iMatrixSize = pProjector->getVolumeGeometry()->getWindowLengthX() * + pProjector->getVolumeGeometry()->getWindowLengthY(); + int iDetectorCount = pProjector->getProjectionGeometry()->getDetectorCount(); + int iTotalStoredPixelCount = 0; + for (int i = 0; i < iDetectorCount; i++) { + iTotalStoredPixelCount += piRayStoredPixelCount[i]; + } + + // create matlab sparse matrix + plhs[0] = mxCreateSparse(iMatrixSize, // number of rows (#pixels) + iDetectorCount, // number of columns (#detectors) + iTotalStoredPixelCount, // number of non-zero elements + mxREAL); // element type + double* values = mxGetPr(plhs[0]); + mwIndex* rows = mxGetIr(plhs[0]); + mwIndex* cols = mxGetJc(plhs[0]); + + int currentBase = 0; + int currentIndex = 0; + for (int i = 0; i < iDetectorCount; i++) { + for (int j = 0; j < piRayStoredPixelCount[i]; j++) { + values[currentIndex + j] = pPixelsWheights[currentBase + j].m_fWeight; + rows[currentIndex + j] = pPixelsWheights[currentBase + j].m_iIndex; + } + + currentBase += pProjector->getProjectionWeightsCount(iProjectionIndex) / pProjector->getProjectionGeometry()->getDetectorCount(); + currentIndex += piRayStoredPixelCount[i]; + } + cols[0] = piRayStoredPixelCount[0]; + for (int j = 1; j < iDetectorCount; j++) { + cols[j] = cols[j-1] + piRayStoredPixelCount[j]; + } + cols[iDetectorCount] = iTotalStoredPixelCount; + } + + delete[] pPixelsWheights; + delete[] piRayStoredPixelCount; +} + +//----------------------------------------------------------------------------------------- +/** hit_detectors = astra_mex_projector('splat', id, col, row); + * + * Create a list of detector indices which have a nonzero contribution to the projection matrix for a pixel [row,col]. + * id: identifier of the projector object as stored in the astra-library. + * col: column of the pixel + * row: row of the pixel + * hit_detectors: list of detector indices [angle_index, det_index] that are hit + */ +void astra_mex_projector_splat(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: get input + if (nrhs < 4) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iPid = (int)(mxGetScalar(prhs[1])); + int iX = (int)(mxGetScalar(prhs[2])); + int iY = (int)(mxGetScalar(prhs[3])); + + // step2: get projector + CProjector2D* pProjector = CProjector2DManager::getSingleton().get(iPid); + if (!pProjector || !pProjector->isInitialized()) { + mexErrMsgTxt("Projector not initialized.\n"); + return; + } + + // step3: perform action + vector detinfo = pProjector->projectPoint(iX, iY); + + // step4: output + if (nlhs <= 1) { + plhs[0] = mxCreateDoubleMatrix(detinfo.size(), // # rows + 2, // # cols + mxREAL); // datatype 32-bits + double* out = mxGetPr(plhs[0]); + + // fill up output + int i = 0; + int x; + for (x = 0; x < detinfo.size(); ++x) { + out[i] = detinfo[x].m_iAngleIndex; + ++i; + } + for (x = 0; x < detinfo.size(); ++x) { + out[i] = detinfo[x].m_iDetectorIndex; + ++i; + } + } + + +} + +//----------------------------------------------------------------------------------------- +/** matrix_id = astra_mex_projector('matrix', id); + * + * Create an explicit projection matrix for this projector. + * It returns an ID usable with astra_mex_matrix(). + * id: identifier of the projector object as stored in the astra-library. + */ +void astra_mex_projector_matrix(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: get input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iPid = (int)(mxGetScalar(prhs[1])); + + // step2: get projector + CProjector2D* pProjector = CProjector2DManager::getSingleton().get(iPid); + if (!pProjector || !pProjector->isInitialized()) { + mexErrMsgTxt("Projector not initialized.\n"); + return; + } + + CSparseMatrix* pMatrix = pProjector->getMatrix(); + if (!pMatrix || !pMatrix->isInitialized()) { + mexErrMsgTxt("Couldn't initialize data object.\n"); + delete pMatrix; + return; + } + + // store data object + int iIndex = CMatrixManager::getSingleton().store(pMatrix); + + // return data id + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(iIndex); + } +} + +//----------------------------------------------------------------------------------------- +/** result = astra_mex_projector('is_cuda', id); + * + * Return is the specified projector is a cuda projector. + * id: identifier of the projector object as stored in the astra-library. + */ +void astra_mex_projector_is_cuda(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) +{ + // step1: get input + if (nrhs < 2) { + mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); + return; + } + int iPid = (int)(mxGetScalar(prhs[1])); + + // step2: get projector + CProjector2D* pProjector = CProjector2DManager::getSingleton().get(iPid); + if (!pProjector || !pProjector->isInitialized()) { + mexErrMsgTxt("Projector not initialized.\n"); + return; + } + +#ifdef ASTRA_CUDA + CCudaProjector2D* pCP = dynamic_cast(pProjector); + plhs[0] = mxCreateLogicalScalar(pCP ? 1 : 0); +#else + plhs[0] = mxCreateLogicalScalar(0); +#endif +} + + + +//----------------------------------------------------------------------------------------- + +static void printHelp() +{ + mexPrintf("Please specify a mode of operation.\n"); + mexPrintf("Valid modes: create, delete, clear, info, projection_geometry,\n"); + mexPrintf(" volume_geometry, weights_single_ray, weights_projection\n"); + mexPrintf(" splat, matrix, is_cuda\n"); +} + + +//----------------------------------------------------------------------------------------- +/** + * ... = astra_mex_projector(mode, ...); + */ +void mexFunction(int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + // INPUT: Mode + string sMode = ""; + if (1 <= nrhs) { + sMode = mex_util_get_string(prhs[0]); + } else { + printHelp(); + return; + } + + // SWITCH (MODE) + if (sMode == "create") { + astra_mex_projector_create(nlhs, plhs, nrhs, prhs); + } else if (sMode == "delete") { + astra_mex_projector_delete(nlhs, plhs, nrhs, prhs); + } else if (sMode == "clear") { + astra_mex_projector_clear(nlhs, plhs, nrhs, prhs); + } else if (sMode == "info") { + astra_mex_projector_info(nlhs, plhs, nrhs, prhs); + } else if (sMode == "projection_geometry") { + astra_mex_projector_projection_geometry(nlhs, plhs, nrhs, prhs); + } else if (sMode == "volume_geometry") { + astra_mex_projector_volume_geometry(nlhs, plhs, nrhs, prhs); + } else if (sMode == "weights_single_ray") { + astra_mex_projector_weights_single_ray(nlhs, plhs, nrhs, prhs); + } else if (sMode == "weights_projection") { + astra_mex_projector_weights_projection(nlhs, plhs, nrhs, prhs); + } else if (sMode == "splat") { + astra_mex_projector_splat(nlhs, plhs, nrhs, prhs); + } else if (sMode == "matrix") { + astra_mex_projector_matrix(nlhs, plhs, nrhs, prhs); + } else if (sMode == "is_cuda") { + astra_mex_projector_is_cuda(nlhs, plhs, nrhs, prhs); + } else { + printHelp(); + } + return; +} + + diff --git a/matlab/mex/astra_mex_projector_vc08.vcproj b/matlab/mex/astra_mex_projector_vc08.vcproj new file mode 100644 index 0000000..1380061 --- /dev/null +++ b/matlab/mex/astra_mex_projector_vc08.vcproj @@ -0,0 +1,591 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/matlab/mex/astra_mex_vc08.vcproj b/matlab/mex/astra_mex_vc08.vcproj new file mode 100644 index 0000000..58c1e0a --- /dev/null +++ b/matlab/mex/astra_mex_vc08.vcproj @@ -0,0 +1,591 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/matlab/mex/mex.def b/matlab/mex/mex.def new file mode 100644 index 0000000..c54c4e0 --- /dev/null +++ b/matlab/mex/mex.def @@ -0,0 +1 @@ +EXPORTS mexFunction \ No newline at end of file diff --git a/matlab/mex/mexHelpFunctions.cpp b/matlab/mex/mexHelpFunctions.cpp new file mode 100644 index 0000000..4105ee1 --- /dev/null +++ b/matlab/mex/mexHelpFunctions.cpp @@ -0,0 +1,642 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +$Id$ +*/ + +/** \file mexHelpFunctions.cpp + * + * \brief Contains some functions for interfacing matlab with c data structures + */ +#include "mexHelpFunctions.h" + +#include "astra/SparseMatrixProjectionGeometry2D.h" +#include "astra/FanFlatVecProjectionGeometry2D.h" +#include "astra/AstraObjectManager.h" + +using namespace std; +using namespace astra; + + +//----------------------------------------------------------------------------------------- +// get string from matlab +std::string mex_util_get_string(const mxArray* pInput) +{ + if (!mxIsChar(pInput)) { + return ""; + } + mwSize iLength = mxGetNumberOfElements(pInput) + 1; + char* buf = new char[iLength]; + mxGetString(pInput, buf, iLength); + std::string res = std::string(buf); + delete[] buf; + return res; +} + +//----------------------------------------------------------------------------------------- +// is option +bool isOption(std::list lOptions, std::string sOption) +{ + return std::find(lOptions.begin(), lOptions.end(), sOption) != lOptions.end(); +} + +//----------------------------------------------------------------------------------------- +// turn a matlab struct into a c++ map +std::map parseStruct(const mxArray* pInput) +{ + std::map res; + + // check type + if (!mxIsStruct(pInput)) { + mexErrMsgTxt("Input must be a struct."); + return res; + } + + // get field names + int nfields = mxGetNumberOfFields(pInput); + for (int i = 0; i < nfields; i++) { + std::string sFieldName = std::string(mxGetFieldNameByNumber(pInput, i)); + res[sFieldName] = mxGetFieldByNumber(pInput,0,i); + } + return res; +} + +//----------------------------------------------------------------------------------------- +// turn a c++ map into a matlab struct +mxArray* buildStruct(std::map mInput) +{ + mwSize dims[2] = {1, 1}; + mxArray* res = mxCreateStructArray(2,dims,0,0); + + for (std::map::iterator it = mInput.begin(); it != mInput.end(); it++) { + mxAddField(res, (*it).first.c_str()); + mxSetField(res, 0, (*it).first.c_str(), (*it).second); + } + return res; +} + +//----------------------------------------------------------------------------------------- +// parse projection geometry data +astra::CProjectionGeometry2D* parseProjectionGeometryStruct(const mxArray* prhs) +{ + // parse struct + std::map mStruct = parseStruct(prhs); + + // create projection geometry object + string type = mex_util_get_string(mStruct["type"]); + if (type == "parallel") { + + // detector_width + float32 fDetWidth = 1.0f; + mxArray* tmp = mStruct["detector_width"]; + if (tmp != NULL) { + fDetWidth = (float32)(mxGetScalar(tmp)); + } + + // detector_count + int iDetCount = 100; + tmp = mStruct["detector_count"]; + if (tmp != NULL) { + iDetCount = (int)(mxGetScalar(tmp)); + } + + // angles + float32* pfAngles; + int iAngleCount; + tmp = mStruct["projection_angles"]; + if (tmp != NULL) { + double* angleValues = mxGetPr(tmp); + iAngleCount = mxGetN(tmp) * mxGetM(tmp); + pfAngles = new float32[iAngleCount]; + for (int i = 0; i < iAngleCount; i++) { + pfAngles[i] = angleValues[i]; + } + } else { + mexErrMsgTxt("'angles' not specified, error."); + return NULL; + } + + // create projection geometry + return new astra::CParallelProjectionGeometry2D(iAngleCount, // number of projections + iDetCount, // number of detectors + fDetWidth, // width of the detectors + pfAngles); // angles array + } + + else if (type == "fanflat") { + + // detector_width + float32 fDetWidth = 1.0f; + mxArray* tmp = mStruct["detector_width"]; + if (tmp != NULL) { + fDetWidth = (float32)(mxGetScalar(tmp)); + } + + // detector_count + int iDetCount = 100; + tmp = mStruct["detector_count"]; + if (tmp != NULL) { + iDetCount = (int)(mxGetScalar(tmp)); + } + + // angles + float32* pfAngles; + int iAngleCount; + tmp = mStruct["projection_angles"]; + if (tmp != NULL) { + double* angleValues = mxGetPr(tmp); + iAngleCount = mxGetN(tmp) * mxGetM(tmp); + pfAngles = new float32[iAngleCount]; + for (int i = 0; i < iAngleCount; i++) { + pfAngles[i] = angleValues[i]; + } + } else { + mexErrMsgTxt("'angles' not specified, error."); + return NULL; + } + + // origin_source_dist + int iDistOriginSource = 100; + tmp = mStruct["origin_source_dist"]; + if (tmp != NULL) { + iDistOriginSource = (int)(mxGetScalar(tmp)); + } + + // origin_det_dist + int iDistOriginDet = 100; + tmp = mStruct["origin_det_dist"]; + if (tmp != NULL) { + iDistOriginDet = (int)(mxGetScalar(tmp)); + } + + // create projection geometry + return new astra::CFanFlatProjectionGeometry2D(iAngleCount, // number of projections + iDetCount, // number of detectors + fDetWidth, // width of the detectors + pfAngles, // angles array + iDistOriginSource, // distance origin source + iDistOriginDet); // distance origin detector + } + + else { + mexPrintf("Only parallel and fanflat projection geometry implemented."); + return NULL; + } +} + +//----------------------------------------------------------------------------------------- +// create projection geometry data +mxArray* createProjectionGeometryStruct(astra::CProjectionGeometry2D* _pProjGeom) +{ + // temporary map to store the data for the MATLAB struct + std::map mGeometryInfo; + + // detectorCount + mGeometryInfo["DetectorCount"] = mxCreateDoubleScalar(_pProjGeom->getDetectorCount()); + + if (!_pProjGeom->isOfType("fanflat_vec")) { + // detectorWidth + mGeometryInfo["DetectorWidth"] = mxCreateDoubleScalar(_pProjGeom->getDetectorWidth()); + + // pfProjectionAngles + mxArray* pAngles = mxCreateDoubleMatrix(1, _pProjGeom->getProjectionAngleCount(), mxREAL); + double* out = mxGetPr(pAngles); + for (int i = 0; i < _pProjGeom->getProjectionAngleCount(); i++) { + out[i] = _pProjGeom->getProjectionAngle(i); + } + mGeometryInfo["ProjectionAngles"] = pAngles; + } + else { + astra::CFanFlatVecProjectionGeometry2D* pVecGeom = dynamic_cast(_pProjGeom); + mxArray* pVectors = mxCreateDoubleMatrix(1, pVecGeom->getProjectionAngleCount()*6, mxREAL); + double* out = mxGetPr(pVectors); + int iDetCount = pVecGeom->getDetectorCount(); + for (int i = 0; i < pVecGeom->getProjectionAngleCount(); i++) { + const SFanProjection* p = &pVecGeom->getProjectionVectors()[i]; + out[6*i + 0] = p->fSrcX; + out[6*i + 1] = p->fSrcY; + out[6*i + 2] = p->fDetSX + 0.5f*iDetCount*p->fDetUX; + out[6*i + 3] = p->fDetSY + 0.5f*iDetCount*p->fDetUY; + out[6*i + 4] = p->fDetUX; + out[6*i + 5] = p->fDetUY; + } + mGeometryInfo["Vectors"] = pVectors; + } + + // parallel specific options + if (_pProjGeom->isOfType("parallel")) { + // type + mGeometryInfo["type"] = mxCreateString("parallel"); + } + // fanflat specific options + else if (_pProjGeom->isOfType("fanflat")) { + astra::CFanFlatProjectionGeometry2D* pFanFlatGeom = dynamic_cast(_pProjGeom); + // detectorCount + mGeometryInfo["DistanceOriginSource"] = mxCreateDoubleScalar(pFanFlatGeom->getOriginSourceDistance()); + // detectorWidth + mGeometryInfo["DistanceOriginDetector"] = mxCreateDoubleScalar(pFanFlatGeom->getOriginDetectorDistance()); + // type + mGeometryInfo["type"] = mxCreateString("fanflat"); + } + else if (_pProjGeom->isOfType("sparse_matrix")) { + astra::CSparseMatrixProjectionGeometry2D* pSparseMatrixGeom = dynamic_cast(_pProjGeom); + mGeometryInfo["type"] = mxCreateString("sparse_matrix"); + mGeometryInfo["MatrixID"] = mxCreateDoubleScalar(CMatrixManager::getSingleton().getIndex(pSparseMatrixGeom->getMatrix())); + } + else if(_pProjGeom->isOfType("fanflat_vec")) { + mGeometryInfo["type"] = mxCreateString("fanflat_vec"); + } + + // build and return the MATLAB struct + return buildStruct(mGeometryInfo); +} + +//----------------------------------------------------------------------------------------- +// parse reconstruction geometry data +astra::CVolumeGeometry2D* parseVolumeGeometryStruct(const mxArray* prhs) +{ + // parse struct + std::map mStruct = parseStruct(prhs); + + std::map mOptions = parseStruct(mStruct["option"]); + + // GridColCount + int iWindowColCount = 128; + mxArray* tmp = mStruct["GridColCount"]; + if (tmp != NULL) { + iWindowColCount = (int)(mxGetScalar(tmp)); + } + + // GridRowCount + int iWindowRowCount = 128; + tmp = mStruct["GridRowCount"]; + if (tmp != NULL) { + iWindowRowCount = (int)(mxGetScalar(tmp)); + } + + // WindowMinX + float32 fWindowMinX = - iWindowColCount / 2; + tmp = mOptions["WindowMinX"]; + if (tmp != NULL) { + fWindowMinX = (float32)(mxGetScalar(tmp)); + } + + // WindowMaxX + float32 fWindowMaxX = iWindowColCount / 2; + tmp = mOptions["WindowMaxX"]; + if (tmp != NULL) { + fWindowMaxX = (float32)(mxGetScalar(tmp)); + } + + // WindowMinY + float32 fWindowMinY = - iWindowRowCount / 2; + tmp = mOptions["WindowMinY"]; + if (tmp != NULL) { + fWindowMinY = (float32)(mxGetScalar(tmp)); + } + + // WindowMaxX + float32 fWindowMaxY = iWindowRowCount / 2; + tmp = mOptions["WindowMaxY"]; + if (tmp != NULL) { + fWindowMaxY = (float32)(mxGetScalar(tmp)); + } + + // create and return reconstruction geometry + return new astra::CVolumeGeometry2D(iWindowColCount, iWindowRowCount, + fWindowMinX, fWindowMinY, + fWindowMaxX, fWindowMaxY); +} + +//----------------------------------------------------------------------------------------- +// create reconstruction geometry data +mxArray* createVolumeGeometryStruct(astra::CVolumeGeometry2D* _pReconGeom) +{ + // temporary map to store the data for the MATLAB struct + std::map mGeometryInfo; + + // fill up map + mGeometryInfo["GridColCount"] = mxCreateDoubleScalar(_pReconGeom->getGridColCount()); + mGeometryInfo["GridRowCount"] = mxCreateDoubleScalar(_pReconGeom->getGridRowCount()); + + std::map mGeometryOptions; + mGeometryOptions["WindowMinX"] = mxCreateDoubleScalar(_pReconGeom->getWindowMinX()); + mGeometryOptions["WindowMaxX"] = mxCreateDoubleScalar(_pReconGeom->getWindowMaxX()); + mGeometryOptions["WindowMinY"] = mxCreateDoubleScalar(_pReconGeom->getWindowMinY()); + mGeometryOptions["WindowMaxY"] = mxCreateDoubleScalar(_pReconGeom->getWindowMaxY()); + + mGeometryInfo["option"] = buildStruct(mGeometryOptions); + + // build and return the MATLAB struct + return buildStruct(mGeometryInfo); +} + + +//----------------------------------------------------------------------------------------- +string matlab2string(const mxArray* pField) +{ + // is string? + if (mxIsChar(pField)) { + return mex_util_get_string(pField); + } + + // is scalar? + if (mxIsNumeric(pField) && mxGetM(pField)*mxGetN(pField) == 1) { + return boost::lexical_cast(mxGetScalar(pField)); + } + + return ""; +} + +//----------------------------------------------------------------------------------------- +// Options struct to xml node +bool readOptions(XMLNode* node, const mxArray* pOptionStruct) +{ + // loop all fields + int nfields = mxGetNumberOfFields(pOptionStruct); + for (int i = 0; i < nfields; i++) { + std::string sFieldName = std::string(mxGetFieldNameByNumber(pOptionStruct, i)); + const mxArray* pField = mxGetFieldByNumber(pOptionStruct, 0, i); + + if (node->hasOption(sFieldName)) { + mexErrMsgTxt("Duplicate option"); + return false; + } + + // string or scalar + if (mxIsChar(pField) || mex_is_scalar(pField)) { + string sValue = matlab2string(pField); + node->addOption(sFieldName, sValue); + } else + // numerical array + if (mxIsNumeric(pField) && mxGetM(pField)*mxGetN(pField) > 1) { + if (!mxIsDouble(pField)) { + mexErrMsgTxt("Numeric input must be double."); + return false; + } + + XMLNode* listbase = node->addChildNode("Option"); + listbase->addAttribute("key", sFieldName); + listbase->addAttribute("listsize", mxGetM(pField)*mxGetN(pField)); + double* pdValues = mxGetPr(pField); + int index = 0; + for (unsigned int row = 0; row < mxGetM(pField); row++) { + for (unsigned int col = 0; col < mxGetN(pField); col++) { + XMLNode* item = listbase->addChildNode("ListItem"); + item->addAttribute("index", index); + item->addAttribute("value", pdValues[col*mxGetM(pField)+row]); + index++; + delete item; + } + } + delete listbase; + } else { + mexErrMsgTxt("Unsupported option type"); + return false; + } + } + return true; +} + +//----------------------------------------------------------------------------------------- +// struct to xml node +bool readStruct(XMLNode* root, const mxArray* pStruct) +{ + // loop all fields + int nfields = mxGetNumberOfFields(pStruct); + for (int i = 0; i < nfields; i++) { + + // field and fieldname + std::string sFieldName = std::string(mxGetFieldNameByNumber(pStruct, i)); + const mxArray* pField = mxGetFieldByNumber(pStruct, 0, i); + + // string + if (mxIsChar(pField)) { + string sValue = matlab2string(pField); + if (sFieldName == "type") { + root->addAttribute("type", sValue); + } else { + delete root->addChildNode(sFieldName, sValue); + } + } + + // scalar + if (mex_is_scalar(pField)) { + string sValue = matlab2string(pField); + delete root->addChildNode(sFieldName, sValue); + } + + // numerical array + if (mxIsNumeric(pField) && mxGetM(pField)*mxGetN(pField) > 1) { + if (!mxIsDouble(pField)) { + mexErrMsgTxt("Numeric input must be double."); + return false; + } + XMLNode* listbase = root->addChildNode(sFieldName); + listbase->addAttribute("listsize", mxGetM(pField)*mxGetN(pField)); + double* pdValues = mxGetPr(pField); + int index = 0; + for (unsigned int row = 0; row < mxGetM(pField); row++) { + for (unsigned int col = 0; col < mxGetN(pField); col++) { + XMLNode* item = listbase->addChildNode("ListItem"); + item->addAttribute("index", index); + item->addAttribute("value", pdValues[col*mxGetM(pField)+row]); + index++; + delete item; + } + } + delete listbase; + } + + + // not castable to a single string + if (mxIsStruct(pField)) { + if (sFieldName == "options" || sFieldName == "option" || sFieldName == "Options" || sFieldName == "Option") { + bool ret = readOptions(root, pField); + if (!ret) + return false; + } else { + XMLNode* newNode = root->addChildNode(sFieldName); + bool ret = readStruct(newNode, pField); + delete newNode; + if (!ret) + return false; + } + } + + } + + return true; +} + +//----------------------------------------------------------------------------------------- +// turn a MATLAB struct into an XML Document +XMLDocument* struct2XML(string rootname, const mxArray* pStruct) +{ + if (!mxIsStruct(pStruct)) { + mexErrMsgTxt("Input must be a struct."); + return NULL; + } + + // create the document + XMLDocument* doc = XMLDocument::createDocument(rootname); + XMLNode* rootnode = doc->getRootNode(); + + // read the struct + bool ret = readStruct(rootnode, pStruct); + //doc->getRootNode()->print(); + delete rootnode; + + if (!ret) { + delete doc; + doc = 0; + } + + return doc; +} + + + + + +//----------------------------------------------------------------------------------------- +// turn an std vector object to an mxArray +mxArray* vectorToMxArray(std::vector mInput) +{ + mxArray* res = mxCreateDoubleMatrix(1, mInput.size(), mxREAL); + double* pdData = mxGetPr(res); + for (unsigned int i = 0; i < mInput.size(); i++) { + pdData[i] = mInput[i]; + } + return res; +} + +//----------------------------------------------------------------------------------------- +// turn a vector> object to an mxArray +mxArray* vector2DToMxArray(std::vector > mInput) +{ + unsigned int sizex = mInput.size(); + if (sizex == 0) return mxCreateString("empty"); + unsigned int sizey = mInput[0].size(); + + mxArray* res = mxCreateDoubleMatrix(sizex, sizey, mxREAL); + double* pdData = mxGetPr(res); + for (unsigned int i = 0; i < sizex; i++) { + for (unsigned int j = 0; j < sizey && j < mInput[i].size(); j++) { + pdData[j*sizex+i] = mInput[i][j]; + } + } + return res; +} + +//----------------------------------------------------------------------------------------- +// turn a boost::any object to an mxArray +mxArray* anyToMxArray(boost::any _any) +{ + if (_any.type() == typeid(std::string)) { + std::string str = boost::any_cast(_any); + return mxCreateString(str.c_str()); + } + if (_any.type() == typeid(int)) { + return mxCreateDoubleScalar(boost::any_cast(_any)); + } + if (_any.type() == typeid(float32)) { + return mxCreateDoubleScalar(boost::any_cast(_any)); + } + if (_any.type() == typeid(std::vector)) { + return vectorToMxArray(boost::any_cast >(_any)); + } + if (_any.type() == typeid(std::vector >)) { + return vector2DToMxArray(boost::any_cast > >(_any)); + } + return NULL; +} +//----------------------------------------------------------------------------------------- +// return true ig the argument is a scalar +bool mex_is_scalar(const mxArray* pInput) +{ + return (mxIsNumeric(pInput) && mxGetM(pInput)*mxGetN(pInput) == 1); +} + +//----------------------------------------------------------------------------------------- +mxArray* XML2struct(astra::XMLDocument* xml) +{ + XMLNode* node = xml->getRootNode(); + mxArray* str = XMLNode2struct(xml->getRootNode()); + delete node; + return str; +} + +//----------------------------------------------------------------------------------------- +mxArray* XMLNode2struct(astra::XMLNode* node) +{ + std::map mList; + + // type_attribute + if (node->hasAttribute("type")) { + mList["type"] = mxCreateString(node->getAttribute("type").c_str()); + } + + list nodes = node->getNodes(); + for (list::iterator it = nodes.begin(); it != nodes.end(); it++) { + XMLNode* subnode = (*it); + // list + if (subnode->hasAttribute("listsize")) { + cout << "lkmdsqldqsjkl" << endl; + cout << " " << node->getContentNumericalArray().size() << endl; + mList[subnode->getName()] = vectorToMxArray(node->getContentNumericalArray()); + } + // string + else { + mList[subnode->getName()] = mxCreateString(subnode->getContent().c_str()); + } + delete subnode; + } + + return buildStruct(mList); +} + +void get3DMatrixDims(const mxArray* x, mwSize *dims) +{ + const mwSize* mdims = mxGetDimensions(x); + mwSize dimCount = mxGetNumberOfDimensions(x); + if (dimCount == 1) { + dims[0] = mdims[0]; + dims[1] = 1; + dims[2] = 1; + } else if (dimCount == 2) { + dims[0] = mdims[0]; + dims[1] = mdims[1]; + dims[2] = 1; + } else if (dimCount == 3) { + dims[0] = mdims[0]; + dims[1] = mdims[1]; + dims[2] = mdims[2]; + } else { + dims[0] = 0; + dims[1] = 0; + dims[2] = 0; + } +} diff --git a/matlab/mex/mexHelpFunctions.h b/matlab/mex/mexHelpFunctions.h new file mode 100644 index 0000000..425b4ef --- /dev/null +++ b/matlab/mex/mexHelpFunctions.h @@ -0,0 +1,76 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see . + +----------------------------------------------------------------------- +$Id$ +*/ + +#ifndef _INC_ASTRA_MEX_HELPFUNCTIONS +#define _INC_ASTRA_MEX_HELPFUNCTIONS + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "astra/Globals.h" +#include "astra/Utilities.h" + +#include "astra/ParallelProjectionGeometry2D.h" +#include "astra/FanFlatProjectionGeometry2D.h" +#include "astra/VolumeGeometry2D.h" + +#include "astra/XMLDocument.h" +#include "astra/XMLNode.h" + +std::string mex_util_get_string(const mxArray* pInput); +bool isOption(std::list lOptions, std::string sOption); + +bool mex_is_scalar(const mxArray* pInput); + +std::map parseStruct(const mxArray* pInput); +mxArray* buildStruct(std::map mInput); +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* createVolumeGeometryStruct(astra::CVolumeGeometry2D* _pReconGeom); + +astra::XMLDocument* struct2XML(string rootname, const mxArray* pStruct); + +mxArray* XML2struct(astra::XMLDocument* xml); +mxArray* XMLNode2struct(astra::XMLNode* xml); + +void get3DMatrixDims(const mxArray* x, mwSize *dims); + +#endif -- cgit v1.2.3