Add mxarray/link functionality for 3d volumes

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
 
 //-----------------------------------------------------------------------------------------
 /**
@@ -261,6 +261,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
@@ -70,6 +70,15 @@ CFloat32ProjectionData3DMemory::CFloat32ProjectionData3DMemory(CProjectionGeomet
 }
 
 //----------------------------------------------------------------------------------------
+// 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