From 01e94c82d907b8d6aa155affc01160396e794b31 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 22 Apr 2014 14:16:06 +0000 Subject: Add mxarray/link functionality for 3d volumes --- include/astra/Float32Data3DMemory.h | 22 +++ include/astra/Float32ProjectionData3DMemory.h | 27 +++ include/astra/Float32VolumeData3DMemory.h | 25 +++ matlab/mex/astra_mex_data3d_c.cpp | 231 +++++++++++++++++++++++++- src/Float32Data3DMemory.cpp | 57 ++++++- src/Float32ProjectionData3DMemory.cpp | 18 ++ src/Float32VolumeData3DMemory.cpp | 17 ++ 7 files changed, 388 insertions(+), 9 deletions(-) diff --git a/include/astra/Float32Data3DMemory.h b/include/astra/Float32Data3DMemory.h index 3a445b6..f72e0ee 100644 --- a/include/astra/Float32Data3DMemory.h +++ b/include/astra/Float32Data3DMemory.h @@ -147,6 +147,25 @@ protected: */ bool _initialize(int _iWidth, int _iHeight, int _iDepth, float32 _fScalar); + /** Initialization. Initializes an instance of the CFloat32Data3DMemory class with pre-allocated memory. + * Can only be called by derived classes. + * + * Initializes an instance of the CFloat32Data3DMemory class. Memory + * is pre-allocated and passed via the abstract CFloat32CustomMemory handle + * class. The handle will be deleted when the memory can be freed. + * You should override the destructor to provide custom behaviour on free. + * If the object has been initialized before, the + * object is reinitialized and memory is freed and reallocated if necessary. + * This function does not set m_bInitialized to true if everything is ok. + * + * @param _iWidth width of the 2D data (x-axis), must be > 0 + * @param _iHeight height of the 2D data (y-axis), must be > 0 + * @param _iDepth depth of the 2D data (z-axis), must be > 0 + * @param _pCustomMemory the custom memory handle + */ + + bool _initialize(int _iWidth, int _iHeight, int _iDepth, CFloat32CustomMemory* _pCustomMemory); + public: /** Default constructor. Sets all numeric member variables to 0 and all pointer member variables to NULL. @@ -270,6 +289,9 @@ public: * @return l-value */ virtual CFloat32Data3D& clampMax(float32& _fMax); + +private: + CFloat32CustomMemory* m_pCustomMemory; }; diff --git a/include/astra/Float32ProjectionData3DMemory.h b/include/astra/Float32ProjectionData3DMemory.h index 8b61d45..fb54425 100644 --- a/include/astra/Float32ProjectionData3DMemory.h +++ b/include/astra/Float32ProjectionData3DMemory.h @@ -98,6 +98,19 @@ public: */ CFloat32ProjectionData3DMemory(CProjectionGeometry3D* _pGeometry, float32 _fScalar); + /** Constructor. Create an instance of the CFloat32ProjectionData3DMemory class with pre-allocated memory. + * + * Creates an instance of the CFloat32ProjectionData3DMemory class. Memory + * is pre-allocated and passed via the abstract CFloat32CustomMemory handle + * class. The handle will be deleted when the memory can be freed. + * You should override the destructor to provide custom behaviour on free. + * + * @param _pGeometry Projection Geometry object. This object will be HARDCOPIED into this class. + * @param _pCustomMemory custom memory handle + * + */ + CFloat32ProjectionData3DMemory(CProjectionGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory); + /** * Destructor. */ @@ -140,6 +153,20 @@ public: */ bool initialize(CProjectionGeometry3D* _pGeometry, const float32* _pfData); + /** Initialization. Initializes an instance of the CFloat32ProjectionData3DMemory class with pre-allocated memory. + * + * Memory is pre-allocated and passed via the abstract CFloat32CustomMemory handle + * class. The handle will be deleted when the memory can be freed. + * You should override the destructor to provide custom behaviour on free. + * + * @param _pGeometry Projection Geometry object. This object will be HARDCOPIED into this class. + * @param _pCustomMemory custom memory handle + * + */ + bool initialize(CProjectionGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory); + + + /** Fetch a COPY of a projection of the data. Note that if you update the 2D data slice, the data in the * 3D data object will remain unaltered. To copy the data back in the 3D-volume you must return the data by calling 'returnProjection'. * diff --git a/include/astra/Float32VolumeData3DMemory.h b/include/astra/Float32VolumeData3DMemory.h index 51df93e..b3f3891 100644 --- a/include/astra/Float32VolumeData3DMemory.h +++ b/include/astra/Float32VolumeData3DMemory.h @@ -84,6 +84,19 @@ public: */ CFloat32VolumeData3DMemory(CVolumeGeometry3D* _pGeometry, float32 _fScalar); + /** Constructor. Create an instance of the CFloat32VolumeData3DMemory class with pre-allocated memory. + * + * Creates an instance of the CFloat32VolumeData3DMemory class. Memory + * is pre-allocated and passed via the abstract CFloat32CustomMemory handle + * class. The handle will be deleted when the memory can be freed. + * You should override the destructor to provide custom behaviour on free. + * + * @param _pGeometry Volume Geometry object. This object will be HARDCOPIED into this class. + * @param _pCustomMemory custom memory handle + * + */ + CFloat32VolumeData3DMemory(CVolumeGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory); + /** Destructor. */ virtual ~CFloat32VolumeData3DMemory(); @@ -122,6 +135,18 @@ public: */ bool initialize(CVolumeGeometry3D* _pGeometry, float32 _fScalar); + /** Initialization. Initializes an instance of the CFloat32VolumeData3DMemory class with pre-allocated memory. + * + * Memory is pre-allocated and passed via the abstract CFloat32CustomMemory handle + * class. The handle will be deleted when the memory can be freed. + * You should override the destructor to provide custom behaviour on free. + * + * @param _pGeometry Volume Geometry object. This object will be HARDCOPIED into this class. + * @param _pCustomMemory custom memory handle + * + */ + bool initialize(CVolumeGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory); + /** Which type is this class? * * @return DataType: VOLUME diff --git a/matlab/mex/astra_mex_data3d_c.cpp b/matlab/mex/astra_mex_data3d_c.cpp index 1af8844..dfb3003 100644 --- a/matlab/mex/astra_mex_data3d_c.cpp +++ b/matlab/mex/astra_mex_data3d_c.cpp @@ -53,7 +53,7 @@ $Id$ using namespace std; using namespace astra; - +#define USE_MATLAB_UNDOCUMENTED //----------------------------------------------------------------------------------------- /** @@ -260,6 +260,231 @@ void astra_mex_data3d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArra } +//----------------------------------------------------------------------------------------- + +/** id = astra_mex_data3d('link', datatype, geometry, data); + * + * Create a new data 3d object in the astra-library. + * type: + * geom: MATLAB struct with the geometry for the data + * data: A MATLAB matrix containing the data. + * This matrix will be edited _in-place_! + * id: identifier of the 3d data object as it is now stored in the astra-library. + */ + +#ifdef USE_MATLAB_UNDOCUMENTED +extern "C" { +mxArray *mxCreateSharedDataCopy(const mxArray *pr); +bool mxUnshareArray(const mxArray *pr, const bool noDeepCopy); +mxArray *mxUnreference(const mxArray *pr); +} + +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) { + 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 + // 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; + + if (mex_is_scalar(_pArray) || !mxIsSingle(_pArray)) { + mexErrMsgTxt("Data must be a single array."); + 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); + } + 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; + } + + //pDataObject3D->updateStatistics(); + + // step4: store data object + int iIndex = CData3DManager::getSingleton().store(pDataObject3D); + + // step5: return data id + if (1 <= nlhs) { + plhs[0] = mxCreateDoubleScalar(iIndex); + } + +} + + +#endif + + //----------------------------------------------------------------------------------------- /** * [id] = astra_mex_io_data('create_cache', config); @@ -989,6 +1214,10 @@ void mexFunction(int nlhs, mxArray* plhs[], // 3D data if (sMode == std::string("create")) { astra_mex_data3d_create(nlhs, plhs, nrhs, prhs); +#ifdef USE_MATLAB_UNDOCUMENTED + } else if (sMode == "link") { + astra_mex_data3d_link(nlhs, plhs, nrhs, prhs); +#endif } else if (sMode == std::string("create_cache")) { astra_mex_data3d_create_cache(nlhs, plhs, nrhs, prhs); } else if (sMode == std::string("get")) { diff --git a/src/Float32Data3DMemory.cpp b/src/Float32Data3DMemory.cpp index 386770c..4522864 100644 --- a/src/Float32Data3DMemory.cpp +++ b/src/Float32Data3DMemory.cpp @@ -157,6 +157,36 @@ bool CFloat32Data3DMemory::_initialize(int _iWidth, int _iHeight, int _iDepth, f return true; } //---------------------------------------------------------------------------------------- +// Initializes an instance of the CFloat32Data3DMemory class with pre-allocated memory +bool CFloat32Data3DMemory::_initialize(int _iWidth, int _iHeight, int _iDepth, CFloat32CustomMemory* _pCustomMemory) +{ + // basic checks + ASTRA_ASSERT(_iWidth > 0); + ASTRA_ASSERT(_iHeight > 0); + ASTRA_ASSERT(_iDepth > 0); + ASTRA_ASSERT(_pCustomMemory != NULL); + + if (m_bInitialized) { + _unInit(); + } + + // calculate size + m_iWidth = _iWidth; + m_iHeight = _iHeight; + m_iDepth = _iDepth; + m_iSize = (size_t)m_iWidth * m_iHeight * m_iDepth; + + // allocate memory for the data, but do not fill it + m_pCustomMemory = _pCustomMemory; + m_pfData = NULL; + m_ppfDataRowInd = NULL; + m_pppfDataSliceInd = NULL; + _allocateData(); + + // initialization complete + return true; +} +//---------------------------------------------------------------------------------------- //---------------------------------------------------------------------------------------- @@ -172,14 +202,18 @@ void CFloat32Data3DMemory::_allocateData() ASTRA_ASSERT(m_ppfDataRowInd == NULL); ASTRA_ASSERT(m_pppfDataSliceInd == NULL); - // allocate contiguous block + if (!m_pCustomMemory) { + // allocate contiguous block #ifdef _MSC_VER - m_pfData = (float32*)_aligned_malloc(m_iSize * sizeof(float32), 16); + m_pfData = (float32*)_aligned_malloc(m_iSize * sizeof(float32), 16); #else - int ret = posix_memalign((void**)&m_pfData, 16, m_iSize * sizeof(float32)); - ASTRA_ASSERT(ret == 0); + int ret = posix_memalign((void**)&m_pfData, 16, m_iSize * sizeof(float32)); + ASTRA_ASSERT(ret == 0); #endif - ASTRA_ASSERT(((size_t)m_pfData & 15) == 0); + ASTRA_ASSERT(((size_t)m_pfData & 15) == 0); + } else { + m_pfData = m_pCustomMemory->m_fPtr; + } // create array of pointers to each row of the data block m_ppfDataRowInd = new float32*[m_iHeight*m_iDepth]; @@ -209,12 +243,18 @@ void CFloat32Data3DMemory::_freeData() delete[] m_pppfDataSliceInd; // free memory for index table delete[] m_ppfDataRowInd; - // free memory for data block + + if (!m_pCustomMemory) { + // free memory for data block #ifdef _MSC_VER - _aligned_free(m_pfData); + _aligned_free(m_pfData); #else - free(m_pfData); + free(m_pfData); #endif + } else { + delete m_pCustomMemory; + m_pCustomMemory = 0; + } } //---------------------------------------------------------------------------------------- @@ -229,6 +269,7 @@ void CFloat32Data3DMemory::_clear() m_pfData = NULL; m_ppfDataRowInd = NULL; m_pppfDataSliceInd = NULL; + m_pCustomMemory = NULL; //m_fGlobalMin = 0.0f; //m_fGlobalMax = 0.0f; diff --git a/src/Float32ProjectionData3DMemory.cpp b/src/Float32ProjectionData3DMemory.cpp index 4d23688..6d89f95 100644 --- a/src/Float32ProjectionData3DMemory.cpp +++ b/src/Float32ProjectionData3DMemory.cpp @@ -69,6 +69,15 @@ CFloat32ProjectionData3DMemory::CFloat32ProjectionData3DMemory(CProjectionGeomet m_bInitialized = initialize(_pGeometry, _fScalar); } +//---------------------------------------------------------------------------------------- +// Create an instance of the CFloat32ProjectionData2D class with pre-allocated data +CFloat32ProjectionData3DMemory::CFloat32ProjectionData3DMemory(CProjectionGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory) +{ + m_bInitialized = false; + m_bInitialized = initialize(_pGeometry, _pCustomMemory); +} + + //---------------------------------------------------------------------------------------- // Initialization bool CFloat32ProjectionData3DMemory::initialize(CProjectionGeometry3D* _pGeometry) @@ -104,6 +113,15 @@ bool CFloat32ProjectionData3DMemory::initialize(CProjectionGeometry3D* _pGeometr return m_bInitialized; } +//---------------------------------------------------------------------------------------- +// Initialization +bool CFloat32ProjectionData3DMemory::initialize(CProjectionGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory) +{ + m_pGeometry = _pGeometry->clone(); + m_bInitialized = _initialize(m_pGeometry->getDetectorColCount(), m_pGeometry->getProjectionCount(), m_pGeometry->getDetectorRowCount(), _pCustomMemory); + return m_bInitialized; +} + //---------------------------------------------------------------------------------------- // Destructor diff --git a/src/Float32VolumeData3DMemory.cpp b/src/Float32VolumeData3DMemory.cpp index aa3832b..96119f5 100644 --- a/src/Float32VolumeData3DMemory.cpp +++ b/src/Float32VolumeData3DMemory.cpp @@ -67,6 +67,14 @@ CFloat32VolumeData3DMemory::CFloat32VolumeData3DMemory(CVolumeGeometry3D* _pGeom m_bInitialized = false; m_bInitialized = initialize(_pGeometry, _fScalar); } +//---------------------------------------------------------------------------------------- +// Create an instance of the CFloat32VolumeData2D class with pre-allocated data +CFloat32VolumeData3DMemory::CFloat32VolumeData3DMemory(CVolumeGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory) +{ + m_bInitialized = false; + m_bInitialized = initialize(_pGeometry, _pCustomMemory); +} + //---------------------------------------------------------------------------------------- // Destructor @@ -105,6 +113,15 @@ bool CFloat32VolumeData3DMemory::initialize(CVolumeGeometry3D* _pGeometry, float m_bInitialized = _initialize(m_pGeometry->getGridColCount(), m_pGeometry->getGridRowCount(), m_pGeometry->getGridSliceCount(), _fScalar); return m_bInitialized; } +//---------------------------------------------------------------------------------------- +// Initialization +bool CFloat32VolumeData3DMemory::initialize(CVolumeGeometry3D* _pGeometry, CFloat32CustomMemory* _pCustomMemory) +{ + m_pGeometry = _pGeometry->clone(); + m_bInitialized = _initialize(m_pGeometry->getGridColCount(), m_pGeometry->getGridRowCount(), m_pGeometry->getGridSliceCount(), _pCustomMemory); + return m_bInitialized; +} + //---------------------------------------------------------------------------------------- // Fetch a slice -- cgit v1.2.3