Refactor astra_mex_data3d

(Modified) patch by Nicola Viganò

6 files changed, 827 insertions, 493 deletions

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 <mex.h>
 #include "mexHelpFunctions.h"
+#include "mexCopyDataHelpFunctions.h"
+#include "mexDataManagerHelpFunctions.h"
 
 #include <list>
 
@@ -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<CFloat32Data3DMemory*>(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<mxDOUBLE_CLASS>(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<CFloat32Data3DMemory*>(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<mxSINGLE_CLASS>(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<CFloat32Data3DMemory*>(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 <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+$Id$
+*/
+
+#include "mexCopyDataHelpFunctions.h"
+
+#include "mexHelpFunctions.h"
+
+#define HAVE_OMP
+
+#define ROUND_DOWN(x, s) ((x) & ~((s)-1))
+
+#ifdef HAVE_OMP
+#	include <omp.h>
+#endif
+
+#if defined(__SSE2__)
+# include <emmintrin.h>
+
+# 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 <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+$Id$
+*/
+
+#ifndef MEXCOPYDATAHELPFUNCTIONS_H_
+#define MEXCOPYDATAHELPFUNCTIONS_H_
+
+#include <mex.h>
+
+#include "astra/Globals.h"
+
+#include <vector>
+
+void copyMexToCFloat32Array(const mxArray * const, astra::float32 * const,
+		const size_t &);
+void copyCFloat32ArrayToMex(const astra::float32 * const, mxArray * const);
+
+template<mxClassID MType, class AType>
+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 <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+$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::CFloat32Data3DMemory *>(
+			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<astra::CFloat32Data3DMemory *> & 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<astra::CFloat32Data3DMemory *> & vecIn,
+		std::vector<astra::float32 *> & 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<astra::CFloat32Data3DMemory *> & vecIn,
+		std::vector<size_t> & 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 <http://www.gnu.org/licenses/>.
+
+-----------------------------------------------------------------------
+$Id$
+*/
+
+#ifndef MEXDATAMANAGERHELPFUNCTIONS_H_
+#define MEXDATAMANAGERHELPFUNCTIONS_H_
+
+#include <mex.h>
+
+#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<astra::CFloat32Data3DMemory *> &);
+
+void getDataPointers(const std::vector<astra::CFloat32Data3DMemory *> &,
+		std::vector<astra::float32 *> &);
+void getDataSizes(const std::vector<astra::CFloat32Data3DMemory *> &,
+		std::vector<size_t> &);
+
+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<mxClassID datatype>
+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<datatype>(pDataObject);
+		copyCFloat32ArrayToMex(pDataObject->getData(), plhs[0]);
+	}
+}
+
+#endif /* MEXDATAMANAGERHELPFUNCTIONS_H_ */