From 4f525533a2037ad9d0fc8849e38ec72c32991bb4 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 30 Oct 2014 17:41:58 +0100 Subject: Refactor astra_mex_data3d MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit (Modified) patch by Nicola Viganò --- build/linux/Makefile.in | 6 +- matlab/mex/astra_mex_data3d_c.cpp | 543 +++-------------------------- matlab/mex/mexCopyDataHelpFunctions.cpp | 257 ++++++++++++++ matlab/mex/mexCopyDataHelpFunctions.h | 53 +++ matlab/mex/mexDataManagerHelpFunctions.cpp | 371 ++++++++++++++++++++ matlab/mex/mexDataManagerHelpFunctions.h | 90 +++++ 6 files changed, 827 insertions(+), 493 deletions(-) create mode 100644 matlab/mex/mexCopyDataHelpFunctions.cpp create mode 100644 matlab/mex/mexCopyDataHelpFunctions.h create mode 100644 matlab/mex/mexDataManagerHelpFunctions.cpp create mode 100644 matlab/mex/mexDataManagerHelpFunctions.h diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in index d39c044..685e1e5 100644 --- a/build/linux/Makefile.in +++ b/build/linux/Makefile.in @@ -35,6 +35,8 @@ endif ifeq ($(matlab),yes) CPPFLAGS+=-I$(MATLAB_ROOT)/extern/include -DMATLAB_MEX_FILE +CXXFLAGS+=-fopenmp +LDFLAGS+=-fopenmp endif BOOST_CPPFLAGS= @@ -197,7 +199,9 @@ TEST_OBJECTS=\ tests/test_XMLDocument.o MATLAB_CXX_OBJECTS=\ - matlab/mex/mexHelpFunctions.o + matlab/mex/mexHelpFunctions.o \ + matlab/mex/mexCopyDataHelpFunctions.o \ + matlab/mex/mexDataManagerHelpFunctions.o MATLAB_MEX=\ matlab/mex/astra_mex_algorithm_c.$(MEXSUFFIX) \ diff --git a/matlab/mex/astra_mex_data3d_c.cpp b/matlab/mex/astra_mex_data3d_c.cpp index a40e814..48893e7 100644 --- a/matlab/mex/astra_mex_data3d_c.cpp +++ b/matlab/mex/astra_mex_data3d_c.cpp @@ -32,6 +32,8 @@ $Id$ */ #include #include "mexHelpFunctions.h" +#include "mexCopyDataHelpFunctions.h" +#include "mexDataManagerHelpFunctions.h" #include @@ -68,186 +70,38 @@ void astra_mex_data3d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArra return; } - string sDataType = mex_util_get_string(prhs[1]); - CFloat32Data3DMemory* pDataObject3D = NULL; + const mxArray * const geometry = prhs[2]; + const mxArray * const data = nrhs > 3 ? prhs[3] : NULL; - if (nrhs >= 4 && !(mex_is_scalar(prhs[3]) || mxIsDouble(prhs[3]) || mxIsSingle(prhs[3]))) { - mexErrMsgTxt("Data must be single or double."); + if (!checkStructs(geometry)) { + mexErrMsgTxt("Argument 3 is not a valid MATLAB struct.\n"); 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); - delete 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"); + if (data && !checkDataType(data)) { + mexErrMsgTxt("Data must be single or double."); return; } - // Check initialization - if (!pDataObject3D->isInitialized()) { - mexErrMsgTxt("Couldn't initialize data object.\n"); - delete pDataObject3D; + const string sDataType = mex_util_get_string(prhs[1]); + + // step2: Allocate data + CFloat32Data3DMemory* pDataObject3D = + allocateDataObject(sDataType, geometry, data); + if (!pDataObject3D) { + // Error message was already set by the function 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; - } - } - } + // step3: Initialize data + if (!data) { + mxArray * emptyArray = mxCreateDoubleMatrix(0, 0, mxREAL); + copyMexToCFloat32Array(emptyArray, pDataObject3D->getData(), + pDataObject3D->getSize()); + mxDestroyArray(emptyArray); + } else { + copyMexToCFloat32Array(data, pDataObject3D->getData(), + pDataObject3D->getSize()); } pDataObject3D->updateStatistics(); @@ -258,7 +112,6 @@ void astra_mex_data3d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArra if (1 <= nlhs) { plhs[0] = mxCreateDoubleScalar(iIndex); } - } //----------------------------------------------------------------------------------------- @@ -274,209 +127,38 @@ void astra_mex_data3d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArra */ #ifdef USE_MATLAB_UNDOCUMENTED -extern "C" { -mxArray *mxCreateSharedDataCopy(const mxArray *pr); -bool mxUnshareArray(const mxArray *pr, const bool noDeepCopy); -mxArray *mxUnreference(const mxArray *pr); -#if 0 -// Unsupported in Matlab R2014b -bool mxIsSharedArray(const mxArray *pr); -#endif -} - -class CFloat32CustomMemoryMatlab3D : public CFloat32CustomMemory { -public: - // offset allows linking the data object to a sub-volume (in the z direction) - // offset is measured in floats. - CFloat32CustomMemoryMatlab3D(const mxArray* _pArray, bool bUnshare, size_t iOffset) { - //fprintf(stderr, "Passed:\narray: %p\tdata: %p\n", (void*)_pArray, (void*)mxGetData(_pArray)); - // First unshare the input array, so that we may modify it. - if (bUnshare) { -#if 0 - // Unsupported in Matlab R2014b - if (mxIsSharedArray(_pArray)) { - fprintf(stderr, "Performance note: unsharing shared array in link\n"); - } -#endif - mxUnshareArray(_pArray, false); - //fprintf(stderr, "Unshared:\narray: %p\tdata: %p\n", (void*)_pArray, (void*)mxGetData(_pArray)); - } - // Then create a (persistent) copy so the data won't be deleted - // or changed. - m_pLink = mxCreateSharedDataCopy(_pArray); - //fprintf(stderr, "SharedDataCopy:\narray: %p\tdata: %p\n", (void*)m_pLink, (void*)mxGetData(m_pLink)); - mexMakeArrayPersistent(m_pLink); - m_fPtr = (float *)mxGetData(m_pLink); - m_fPtr += iOffset; - } - virtual ~CFloat32CustomMemoryMatlab3D() { - // destroy the shared array - //fprintf(stderr, "Destroy:\narray: %p\tdata: %p\n", (void*)m_pLink, (void*)mxGetData(m_pLink)); - mxDestroyArray(m_pLink); - } -private: - mxArray* m_pLink; -}; void astra_mex_data3d_link(int& nlhs, mxArray* plhs[], int& nrhs, const mxArray* prhs[]) { - // TODO: Reduce code duplication with _create! - // TODO: Also allow empty argument to let this function create its own mxArray + // TODO: Allow empty argument to let this function create its own mxArray // step1: get datatype if (nrhs < 4) { mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); return; } - if (!mxIsStruct(prhs[2])) { - mexErrMsgTxt("Argument 3 is not a valid MATLAB struct.\n"); - } - - const mxArray* _pArray = prhs[3]; - - string sDataType = mex_util_get_string(prhs[1]); - CFloat32Data3DMemory* pDataObject3D = NULL; + const mxArray * const geometry = prhs[2]; + const mxArray * const data = nrhs > 3 ? prhs[3] : NULL; + const mxArray * const unshare = nrhs > 4 ? prhs[4] : NULL; + const mxArray * const zIndex = nrhs > 5 ? prhs[5] : NULL; - if (mex_is_scalar(_pArray) || !mxIsSingle(_pArray)) { - mexErrMsgTxt("Data must be a single array."); + if (!checkStructs(geometry)) { + mexErrMsgTxt("Argument 3 is not a valid MATLAB struct.\n"); return; } - unsigned int iZ = 0; - bool bUnshare = true; - if (nrhs > 4) { - if (!mex_is_scalar(prhs[4])) { - mexErrMsgTxt("Argument 5 (read-only) must be scalar"); - return; - } - // unshare the array if we're not linking read-only - bUnshare = !(bool)mxGetScalar(prhs[4]); - } - if (nrhs > 5) { - if (!mex_is_scalar(prhs[5])) { - mexErrMsgTxt("Argument 6 (Z) must be scalar"); - return; - } - iZ = (unsigned int)mxGetScalar(prhs[5]); - } - - mwSize dims[3]; - - // SWITCH DataType - if (sDataType == "-vol") { - // Read geometry - 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) { - get3DMatrixDims(prhs[3], dims); - bool ok = pGeometry->getGridColCount() == dims[0] && pGeometry->getGridRowCount() == dims[1]; - if (nrhs == 4) { - ok &= pGeometry->getGridSliceCount() == dims[2]; - } else if (nrhs > 4) { - // sub-volume in Z direction - ok &= iZ + pGeometry->getGridSliceCount() <= dims[2]; - } - if (!ok) { - mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); - delete pGeometry; - return; - } - } - - size_t iOffset = iZ; - iOffset *= dims[0]; - iOffset *= dims[1]; - CFloat32CustomMemoryMatlab3D* pHandle = new CFloat32CustomMemoryMatlab3D(_pArray, bUnshare, iOffset); - - // Initialize data object - pDataObject3D = new CFloat32VolumeData3DMemory(pGeometry, pHandle); - delete pGeometry; - } - else if (sDataType == "-sino" || sDataType == "-proj3d" || 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(); - - // 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) { - get3DMatrixDims(prhs[3], dims); - bool ok = pGeometry->getDetectorColCount() == dims[0] && pGeometry->getProjectionCount() == dims[1]; - if (nrhs == 4) { - ok &= pGeometry->getDetectorRowCount() == dims[2]; - } else if (nrhs > 4) { - // sub-volume in Z direction - ok &= iZ + pGeometry->getDetectorRowCount() <= dims[2]; - } - if (!ok) { - mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); - delete pGeometry; - return; - } - } - - size_t iOffset = iZ; - iOffset *= dims[0]; - iOffset *= dims[1]; - CFloat32CustomMemoryMatlab3D* pHandle = new CFloat32CustomMemoryMatlab3D(_pArray, bUnshare, iOffset); - - // Initialize data object - pDataObject3D = new CFloat32ProjectionData3DMemory(pGeometry, pHandle); - delete pGeometry; - } - else { - mexErrMsgTxt("Invalid datatype. Please specify '-vol' or '-proj3d'. \n"); + if (data && !checkDataType(data)) { + mexErrMsgTxt("Data must be single or double."); return; } - // Check initialization - if (!pDataObject3D->isInitialized()) { - mexErrMsgTxt("Couldn't initialize data object.\n"); - delete pDataObject3D; + string sDataType = mex_util_get_string(prhs[1]); + + // step2: Allocate data + CFloat32Data3DMemory* pDataObject3D = + allocateDataObject(sDataType, geometry, data, unshare, zIndex); + if (!pDataObject3D) { + // Error message was already set by the function return; } @@ -489,10 +171,7 @@ void astra_mex_data3d_link(int& nlhs, mxArray* plhs[], int& nrhs, const mxArray* if (1 <= nlhs) { plhs[0] = mxCreateDoubleScalar(iIndex); } - } - - #endif @@ -547,43 +226,8 @@ void astra_mex_data3d_create_cache(int nlhs, mxArray* plhs[], int nrhs, const mx */ 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; - } - } - } - } - +{ + generic_astra_mex_data3d_get(nlhs, plhs, nrhs, prhs); } //----------------------------------------------------------------------------------------- @@ -596,43 +240,8 @@ void astra_mex_data3d_get(int nlhs, mxArray* plhs[], int nrhs, const mxArray* pr */ 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; - } - } - } - } - +{ + generic_astra_mex_data3d_get(nlhs, plhs, nrhs, prhs); } @@ -652,70 +261,20 @@ void astra_mex_data3d_store(int nlhs, mxArray* plhs[], int nrhs, const mxArray* 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"); + if (!checkDataType(prhs[2])) { + mexErrMsgTxt("Data must be single or double."); return; } - if (!(mex_is_scalar(prhs[2]) || mxIsDouble(prhs[2]) || mxIsSingle(prhs[2]))) { - mexErrMsgTxt("Data must be single or double."); + // step2: get data object + CFloat32Data3DMemory* pDataObject = NULL; + if (!checkID(mxGetScalar(prhs[1]), pDataObject)) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); 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; - } - } - } - } + copyMexToCFloat32Array(prhs[2], pDataObject->getData(), pDataObject->getSize()); pDataObject->updateStatistics(); } diff --git a/matlab/mex/mexCopyDataHelpFunctions.cpp b/matlab/mex/mexCopyDataHelpFunctions.cpp new file mode 100644 index 0000000..bbcad30 --- /dev/null +++ b/matlab/mex/mexCopyDataHelpFunctions.cpp @@ -0,0 +1,257 @@ +/* +----------------------------------------------------------------------- +Copyright 2013 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 "mexCopyDataHelpFunctions.h" + +#include "mexHelpFunctions.h" + +#define HAVE_OMP + +#define ROUND_DOWN(x, s) ((x) & ~((s)-1)) + +#ifdef HAVE_OMP +# include +#endif + +#if defined(__SSE2__) +# include + +# define STORE_32F_64F_CORE8(in, out, count, base) \ + {\ + const __m128 inV0 = *((const __m128 *) &in[count + 0 + base]);\ + const __m128 inV1 = *((const __m128 *) &in[count + 4 + base]);\ +\ + *((__m128d *)&out[count + 0 + base]) = _mm_cvtps_pd(inV0);\ + *((__m128d *)&out[count + 2 + base]) = _mm_cvtps_pd(_mm_movehl_ps(inV0, inV0));\ +\ + *((__m128d *)&out[count + 4 + base]) = _mm_cvtps_pd(inV1);\ + *((__m128d *)&out[count + 6 + base]) = _mm_cvtps_pd(_mm_movehl_ps(inV1, inV1));\ + } +# define STORE_64F_32F_CORE8(in, out, count, base) \ + {\ + const __m128d inV0 = *((const __m128d *) &in[count + 0 + base]);\ + const __m128d inV1 = *((const __m128d *) &in[count + 2 + base]);\ +\ + const __m128d inV2 = *((const __m128d *) &in[count + 4 + base]);\ + const __m128d inV3 = *((const __m128d *) &in[count + 6 + base]);\ +\ + *((__m128 *)&out[count + 0 + base]) = _mm_movelh_ps(_mm_cvtpd_ps(inV0), _mm_cvtpd_ps(inV1));\ + *((__m128 *)&out[count + 4 + base]) = _mm_movelh_ps(_mm_cvtpd_ps(inV2), _mm_cvtpd_ps(inV3));\ + } +# define STORE_32F_32F_CORE8(in, out, count, base) \ + {\ + *((__m128 *)&out[count + 0 + base]) = *((const __m128 *)&in[count + 0 + base]);\ + *((__m128 *)&out[count + 4 + base]) = *((const __m128 *)&in[count + 4 + base]);\ + } +# define STORE_CONST_32F_CORE8(val, out, count, base) \ + {\ + *((__m128 *)&out[count + 0 + base]) = val;\ + *((__m128 *)&out[count + 4 + base]) = val;\ + } +#else +# define STORE_32F_64F_CORE8(in, out, count, base) \ + {\ + out[count + 0 + base] = (double)in[count + 0 + base];\ + out[count + 1 + base] = (double)in[count + 1 + base];\ + out[count + 2 + base] = (double)in[count + 2 + base];\ + out[count + 3 + base] = (double)in[count + 3 + base];\ + out[count + 4 + base] = (double)in[count + 4 + base];\ + out[count + 5 + base] = (double)in[count + 5 + base];\ + out[count + 6 + base] = (double)in[count + 6 + base];\ + out[count + 7 + base] = (double)in[count + 7 + base];\ + } +# define STORE_64F_32F_CORE8(in, out, count, base) \ + {\ + out[count + 0 + base] = (float)in[count + 0 + base];\ + out[count + 1 + base] = (float)in[count + 1 + base];\ + out[count + 2 + base] = (float)in[count + 2 + base];\ + out[count + 3 + base] = (float)in[count + 3 + base];\ + out[count + 4 + base] = (float)in[count + 4 + base];\ + out[count + 5 + base] = (float)in[count + 5 + base];\ + out[count + 6 + base] = (float)in[count + 6 + base];\ + out[count + 7 + base] = (float)in[count + 7 + base];\ + } +# define STORE_32F_32F_CORE8(in, out, count, base) \ + {\ + out[count + 0 + base] = in[count + 0 + base];\ + out[count + 1 + base] = in[count + 1 + base];\ + out[count + 2 + base] = in[count + 2 + base];\ + out[count + 3 + base] = in[count + 3 + base];\ + out[count + 4 + base] = in[count + 4 + base];\ + out[count + 5 + base] = in[count + 5 + base];\ + out[count + 6 + base] = in[count + 6 + base];\ + out[count + 7 + base] = in[count + 7 + base];\ + } +#endif + +const char * warnDataTypeNotSupported = "Data type not supported: nothing was copied"; + +void +_copyMexToCFloat32Array(const mxArray * const inArray, astra::float32 * const out) +{ + const long long tot_size = mxGetNumberOfElements(inArray); + const long long block = 32; + const long long totRoundedPixels = ROUND_DOWN(tot_size, block); + + if (mxIsDouble(inArray)) { + const double * const pdMatlabData = mxGetPr(inArray); + +#pragma omp for nowait + for (long long count = 0; count < totRoundedPixels; count += block) { + STORE_64F_32F_CORE8(pdMatlabData, out, count, 0); + STORE_64F_32F_CORE8(pdMatlabData, out, count, 8); + STORE_64F_32F_CORE8(pdMatlabData, out, count, 16); + STORE_64F_32F_CORE8(pdMatlabData, out, count, 24); + } +#pragma omp for nowait + for (long long count = totRoundedPixels; count < tot_size; count++) { + out[count] = pdMatlabData[count]; + } + } + else if (mxIsSingle(inArray)) { + const float * const pfMatlabData = (const float *)mxGetData(inArray); + +#pragma omp for nowait + for (long long count = 0; count < totRoundedPixels; count += block) { + STORE_32F_32F_CORE8(pfMatlabData, out, count, 0); + STORE_32F_32F_CORE8(pfMatlabData, out, count, 8); + STORE_32F_32F_CORE8(pfMatlabData, out, count, 16); + STORE_32F_32F_CORE8(pfMatlabData, out, count, 24); + } +#pragma omp for nowait + for (long long count = totRoundedPixels; count < tot_size; count++) { + out[count] = pfMatlabData[count]; + } + } + else { +#pragma omp single nowait + mexWarnMsgIdAndTxt("ASTRA_MEX:wrong_datatype", warnDataTypeNotSupported); + } +} + +void +_copyCFloat32ArrayToMex(const astra::float32 * const in, mxArray * const outArray) +{ + const long long tot_size = mxGetNumberOfElements(outArray); + const long long block = 32; + const long long tot_rounded_size = ROUND_DOWN(tot_size, block); + + if (mxIsDouble(outArray)) { + double * const pdMatlabData = mxGetPr(outArray); + +#pragma omp for nowait + for (long long count = 0; count < tot_rounded_size; count += block) { + STORE_32F_64F_CORE8(in, pdMatlabData, count, 0); + STORE_32F_64F_CORE8(in, pdMatlabData, count, 8); + STORE_32F_64F_CORE8(in, pdMatlabData, count, 16); + STORE_32F_64F_CORE8(in, pdMatlabData, count, 24); + } +#pragma omp for nowait + for (long long count = tot_rounded_size; count < tot_size; count++) { + pdMatlabData[count] = in[count]; + } + } + else if (mxIsSingle(outArray)) { + float * const pfMatlabData = (float *) mxGetData(outArray); + +#pragma omp for nowait + for (long long count = 0; count < tot_rounded_size; count += block) { + STORE_32F_32F_CORE8(in, pfMatlabData, count, 0); + STORE_32F_32F_CORE8(in, pfMatlabData, count, 8); + STORE_32F_32F_CORE8(in, pfMatlabData, count, 16); + STORE_32F_32F_CORE8(in, pfMatlabData, count, 24); + } +#pragma omp for nowait + for (long long count = tot_rounded_size; count < tot_size; count++) { + pfMatlabData[count] = in[count]; + } + } + else { +#pragma omp single nowait + mexWarnMsgIdAndTxt("ASTRA_MEX:wrong_datatype", warnDataTypeNotSupported); + } +} + +void +_initializeCFloat32Array(const astra::float32 & val, astra::float32 * const out, + const size_t & tot_size) +{ +#ifdef __SSE2__ + const long long block = 32; + const long long tot_rounded_size = ROUND_DOWN(tot_size, block); + + const __m128 vecVal = _mm_set1_ps(val); + + { +#pragma omp for nowait + for (long long count = 0; count < tot_rounded_size; count += block) { + STORE_CONST_32F_CORE8(vecVal, out, count, 0); + STORE_CONST_32F_CORE8(vecVal, out, count, 8); + STORE_CONST_32F_CORE8(vecVal, out, count, 16); + STORE_CONST_32F_CORE8(vecVal, out, count, 24); + } +#else + const long long tot_rounded_size = 0; + { +#endif +#pragma omp for nowait + for (long long count = tot_rounded_size; count < (long long)tot_size; count++) { + out[count] = val; + } + } +} + +void +copyMexToCFloat32Array(const mxArray * const in, + astra::float32 * const out, const size_t &tot_size) +{ +#pragma omp parallel + { + // fill with scalar value + if (mex_is_scalar(in)) { + astra::float32 fValue = 0.f; + if (!mxIsEmpty(in)) { + fValue = (astra::float32)mxGetScalar(in); + } + _initializeCFloat32Array(fValue, out, tot_size); + } + // fill with array value + else { + _copyMexToCFloat32Array(in, out); + } + } +} + +void +copyCFloat32ArrayToMex(const float * const in, mxArray * const outArray) +{ +#pragma omp parallel + { + _copyCFloat32ArrayToMex(in, outArray); + } +} diff --git a/matlab/mex/mexCopyDataHelpFunctions.h b/matlab/mex/mexCopyDataHelpFunctions.h new file mode 100644 index 0000000..d0cf3c6 --- /dev/null +++ b/matlab/mex/mexCopyDataHelpFunctions.h @@ -0,0 +1,53 @@ +/* +----------------------------------------------------------------------- +Copyright 2013 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 MEXCOPYDATAHELPFUNCTIONS_H_ +#define MEXCOPYDATAHELPFUNCTIONS_H_ + +#include + +#include "astra/Globals.h" + +#include + +void copyMexToCFloat32Array(const mxArray * const, astra::float32 * const, + const size_t &); +void copyCFloat32ArrayToMex(const astra::float32 * const, mxArray * const); + +template +mxArray * createEquivMexArray(const AType * const pDataObj) +{ + mwSize dims[3]; + dims[0] = pDataObj->getWidth(); + dims[1] = pDataObj->getHeight(); + dims[2] = pDataObj->getDepth(); + + return mxCreateNumericArray(3, dims, MType, mxREAL); +} + +#endif /* MEXCOPYDATAHELPFUNCTIONS_H_ */ diff --git a/matlab/mex/mexDataManagerHelpFunctions.cpp b/matlab/mex/mexDataManagerHelpFunctions.cpp new file mode 100644 index 0000000..2985a9d --- /dev/null +++ b/matlab/mex/mexDataManagerHelpFunctions.cpp @@ -0,0 +1,371 @@ +/* +----------------------------------------------------------------------- +Copyright 2013 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 "mexDataManagerHelpFunctions.h" + +#include "mexHelpFunctions.h" + +#include "astra/ParallelProjectionGeometry3D.h" +#include "astra/ParallelVecProjectionGeometry3D.h" +#include "astra/ConeProjectionGeometry3D.h" +#include "astra/ConeVecProjectionGeometry3D.h" +#include "astra/Float32VolumeData3DMemory.h" +#include "astra/Float32ProjectionData3DMemory.h" + +#define USE_MATLAB_UNDOCUMENTED + +#ifdef USE_MATLAB_UNDOCUMENTED +extern "C" { +mxArray *mxCreateSharedDataCopy(const mxArray *pr); +bool mxUnshareArray(const mxArray *pr, const bool noDeepCopy); +mxArray *mxUnreference(const mxArray *pr); +#if 0 +// Unsupported in Matlab R2014b +bool mxIsSharedArray(const mxArray *pr); +#endif +} + +class CFloat32CustomMemoryMatlab3D : public astra::CFloat32CustomMemory { +public: + // offset allows linking the data object to a sub-volume (in the z direction) + // offset is measured in floats. + CFloat32CustomMemoryMatlab3D(const mxArray* _pArray, bool bUnshare, size_t iOffset) + { + // Convert from slice to offset + mwSize dims[3]; + get3DMatrixDims(_pArray, dims); + iOffset *= dims[0]; + iOffset *= dims[1]; + + //fprintf(stderr, "Passed:\narray: %p\tdata: %p\n", (void*)_pArray, (void*)mxGetData(_pArray)); + // First unshare the input array, so that we may modify it. + if (bUnshare) { +#if 0 + // Unsupported in Matlab R2014b + if (mxIsSharedArray(_pArray)) { + fprintf(stderr, "Performance note: unsharing shared array in link\n"); + } +#endif + mxUnshareArray(_pArray, false); + //fprintf(stderr, "Unshared:\narray: %p\tdata: %p\n", (void*)_pArray, (void*)mxGetData(_pArray)); + } + // Then create a (persistent) copy so the data won't be deleted + // or changed. + m_pLink = mxCreateSharedDataCopy(_pArray); + //fprintf(stderr, "SharedDataCopy:\narray: %p\tdata: %p\n", (void*)m_pLink, (void*)mxGetData(m_pLink)); + mexMakeArrayPersistent(m_pLink); + m_fPtr = (float *)mxGetData(m_pLink); + m_fPtr += iOffset; + } + virtual ~CFloat32CustomMemoryMatlab3D() { + // destroy the shared array + //fprintf(stderr, "Destroy:\narray: %p\tdata: %p\n", (void*)m_pLink, (void*)mxGetData(m_pLink)); + mxDestroyArray(m_pLink); + } +private: + mxArray* m_pLink; +}; +#endif + +//----------------------------------------------------------------------------------------- +bool +checkID(const astra::int32 & id, astra::CFloat32Data3DMemory *& pDataObj) +{ + pDataObj = dynamic_cast( + astra::CData3DManager::getSingleton().get(id) ); + return (pDataObj && pDataObj->isInitialized()); +} + +//----------------------------------------------------------------------------------------- +bool +checkDataType(const mxArray * const in) +{ + return (mex_is_scalar(in) || mxIsDouble(in) || mxIsSingle(in)); +} + +//----------------------------------------------------------------------------------------- +bool +checkStructs(const mxArray * const in) +{ + return mxIsStruct(in); +} + +//----------------------------------------------------------------------------------------- +bool +checkDataSize(const mxArray * const mArray, + const astra::CProjectionGeometry3D * const geom) +{ + mwSize dims[3]; + get3DMatrixDims(mArray, dims); + return (geom->getDetectorColCount() == dims[0] + && geom->getProjectionCount() == dims[1] + && geom->getDetectorRowCount() == dims[2]); +} + +//----------------------------------------------------------------------------------------- +bool +checkDataSize(const mxArray * const mArray, + const astra::CVolumeGeometry3D * const geom) +{ + mwSize dims[3]; + get3DMatrixDims(mArray, dims); + return (geom->getGridColCount() == dims[0] + && geom->getGridRowCount() == dims[1] + && geom->getGridSliceCount() == dims[2]); +} + +//----------------------------------------------------------------------------------------- +bool +checkDataSize(const mxArray * const mArray, + const astra::CProjectionGeometry3D * const geom, + const mwIndex & zOffset) +{ + mwSize dims[3]; + get3DMatrixDims(mArray, dims); + return (geom->getDetectorColCount() == dims[0] + && geom->getProjectionCount() == dims[1] + && (zOffset + geom->getDetectorRowCount()) <= dims[2]); +} + +//----------------------------------------------------------------------------------------- +bool +checkDataSize(const mxArray * const mArray, + const astra::CVolumeGeometry3D * const geom, + const mwIndex & zOffset) +{ + mwSize dims[3]; + get3DMatrixDims(mArray, dims); + return (geom->getGridColCount() == dims[0] + && geom->getGridRowCount() == dims[1] + && (zOffset + geom->getGridSliceCount()) <= dims[2]); +} + +//----------------------------------------------------------------------------------------- +void +updateStatistics(const std::vector & vecIn) +{ + const size_t tot_size = vecIn.size(); + for (size_t count = 0; count < tot_size; count++) + { + vecIn[count]->updateStatistics(); + } +} + +//----------------------------------------------------------------------------------------- +void +getDataPointers(const std::vector & vecIn, + std::vector & vecOut) +{ + const size_t tot_size = vecIn.size(); + vecOut.resize(tot_size); + for (size_t count = 0; count < tot_size; count++) + { + vecOut[count] = vecIn[count]->getData(); + } +} + +//----------------------------------------------------------------------------------------- +void +getDataSizes(const std::vector & vecIn, + std::vector & vecOut) +{ + const size_t tot_size = vecIn.size(); + vecOut.resize(tot_size); + for (size_t count = 0; count < tot_size; count++) + { + vecOut[count] = vecIn[count]->getSize(); + } +} + +//----------------------------------------------------------------------------------------- +astra::CFloat32Data3DMemory * +allocateDataObject(const std::string & sDataType, + const mxArray * const geometry, const mxArray * const data, + const mxArray * const unshare, const mxArray * const zIndex) +{ + astra::CFloat32Data3DMemory* pDataObject3D = NULL; + + bool bUnshare = true; + if (unshare) + { + if (!mex_is_scalar(unshare)) + { + mexErrMsgTxt("Argument 5 (read-only) must be scalar"); + return NULL; + } + // unshare the array if we're not linking read-only + bUnshare = !(bool)mxGetScalar(unshare); + } + + mwIndex iZ = 0; + if (zIndex) + { + if (!mex_is_scalar(zIndex)) + { + mexErrMsgTxt("Argument 6 (Z) must be scalar"); + return NULL; + } + iZ = (mwSignedIndex)mxGetScalar(zIndex); + } + + // SWITCH DataType + if (sDataType == "-vol") + { + // Read geometry + astra::XMLDocument* xml = struct2XML("VolumeGeometry", geometry); + if (!xml) { + return NULL; + } + astra::Config cfg; + cfg.self = xml->getRootNode(); + + astra::CVolumeGeometry3D* pGeometry = new astra::CVolumeGeometry3D(); + if (!pGeometry->initialize(cfg)) + { + mexErrMsgTxt("Geometry class not initialized. \n"); + delete pGeometry; + delete xml; + return NULL; + } + delete xml; + + // If data is specified, check dimensions + if (data && !mex_is_scalar(data)) + { + if (! (zIndex + ? checkDataSize(data, pGeometry, iZ) + : checkDataSize(data, pGeometry)) ) + { + mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); + delete pGeometry; + return NULL; + } + } + + // Initialize data object +#ifdef USE_MATLAB_UNDOCUMENTED + if (unshare) { + CFloat32CustomMemoryMatlab3D* pHandle = + new CFloat32CustomMemoryMatlab3D(data, bUnshare, iZ); + + // Initialize data object + pDataObject3D = new astra::CFloat32VolumeData3DMemory(pGeometry, pHandle); + } + else + { + pDataObject3D = new astra::CFloat32VolumeData3DMemory(pGeometry); + } +#else + pDataObject3D = new astra::CFloat32VolumeData3DMemory(pGeometry); +#endif + delete pGeometry; + } + else if (sDataType == "-sino" || sDataType == "-proj3d" || sDataType == "-sinocone") + { + // Read geometry + astra::XMLDocument* xml = struct2XML("ProjectionGeometry", geometry); + if (!xml) { + return NULL; + } + astra::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"); + astra::CProjectionGeometry3D* pGeometry = 0; + if (type == "parallel3d") { + pGeometry = new astra::CParallelProjectionGeometry3D(); + } else if (type == "parallel3d_vec") { + pGeometry = new astra::CParallelVecProjectionGeometry3D(); + } else if (type == "cone") { + pGeometry = new astra::CConeProjectionGeometry3D(); + } else if (type == "cone_vec") { + pGeometry = new astra::CConeVecProjectionGeometry3D(); + } else { + mexErrMsgTxt("Invalid geometry type.\n"); + return NULL; + } + + if (!pGeometry->initialize(cfg)) { + mexErrMsgTxt("Geometry class not initialized. \n"); + delete pGeometry; + delete xml; + return NULL; + } + delete xml; + + // If data is specified, check dimensions + if (data && !mex_is_scalar(data)) + { + if (! (zIndex + ? checkDataSize(data, pGeometry, iZ) + : checkDataSize(data, pGeometry)) ) + { + mexErrMsgTxt("The dimensions of the data do not match those specified in the geometry. \n"); + delete pGeometry; + return NULL; + } + } + + // Initialize data object +#ifdef USE_MATLAB_UNDOCUMENTED + if (unshare) + { + CFloat32CustomMemoryMatlab3D* pHandle = + new CFloat32CustomMemoryMatlab3D(data, bUnshare, iZ); + + // Initialize data object + pDataObject3D = new astra::CFloat32ProjectionData3DMemory(pGeometry, pHandle); + } + else + { + pDataObject3D = new astra::CFloat32ProjectionData3DMemory(pGeometry); + } +#else + pDataObject3D = new astra::CFloat32ProjectionData3DMemory(pGeometry); +#endif + delete pGeometry; + } + else + { + mexErrMsgTxt("Invalid datatype. Please specify '-vol' or '-proj3d'. \n"); + return NULL; + } + + // Check initialization + if (!pDataObject3D->isInitialized()) + { + mexErrMsgTxt("Couldn't initialize data object.\n"); + delete pDataObject3D; + return NULL; + } + + return pDataObject3D; +} + diff --git a/matlab/mex/mexDataManagerHelpFunctions.h b/matlab/mex/mexDataManagerHelpFunctions.h new file mode 100644 index 0000000..36b74d8 --- /dev/null +++ b/matlab/mex/mexDataManagerHelpFunctions.h @@ -0,0 +1,90 @@ +/* +----------------------------------------------------------------------- +Copyright 2013 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 MEXDATAMANAGERHELPFUNCTIONS_H_ +#define MEXDATAMANAGERHELPFUNCTIONS_H_ + +#include + +#include "astra/Globals.h" +#include "astra/AstraObjectManager.h" +#include "astra/Float32Data3DMemory.h" +#include "astra/ProjectionGeometry3D.h" +#include "astra/VolumeGeometry3D.h" + +#include "mexCopyDataHelpFunctions.h" + +bool checkID(const astra::int32 &, astra::CFloat32Data3DMemory *&); + +bool checkDataType(const mxArray * const); +bool checkStructs(const mxArray * const); + +bool checkDataSize(const mxArray * const, const astra::CProjectionGeometry3D * const); +bool checkDataSize(const mxArray * const, const astra::CVolumeGeometry3D * const); +bool checkDataSize(const mxArray * const, const astra::CProjectionGeometry3D * const, + const mwIndex & zOffset); +bool checkDataSize(const mxArray * const, const astra::CVolumeGeometry3D * const, + const mwIndex & zOffset); + +void updateStatistics(const std::vector &); + +void getDataPointers(const std::vector &, + std::vector &); +void getDataSizes(const std::vector &, + std::vector &); + +astra::CFloat32Data3DMemory * allocateDataObject(const std::string & sDataType, + const mxArray * const geometry, const mxArray * const data, + const mxArray * const unshare = NULL, const mxArray * const zOffset = NULL); + +//----------------------------------------------------------------------------------------- +template +void generic_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; + } + + // step2: get data object/s + astra::CFloat32Data3DMemory* pDataObject = NULL; + if (!checkID(mxGetScalar(prhs[1]), pDataObject)) { + mexErrMsgTxt("Data object not found or not initialized properly.\n"); + return; + } + + // create output + if (1 <= nlhs) { + plhs[0] = createEquivMexArray(pDataObject); + copyCFloat32ArrayToMex(pDataObject->getData(), plhs[0]); + } +} + +#endif /* MEXDATAMANAGERHELPFUNCTIONS_H_ */ -- cgit v1.2.3