Merge branch 'master' into volgeom3d

Conflicts:
	src/CudaBackProjectionAlgorithm3D.cpp

15 files changed, 364 insertions, 67 deletions

diff --git a/matlab/mex/astra_mex_algorithm_c.cpp b/matlab/mex/astra_mex_algorithm_c.cpp
index e4afa63..ec7aa72 100644
--- a/matlab/mex/astra_mex_algorithm_c.cpp
+++ b/matlab/mex/astra_mex_algorithm_c.cpp
@@ -81,7 +81,7 @@ void astra_mex_algorithm_create(int nlhs, mxArray* plhs[], int nrhs, const mxArr
 	// turn MATLAB struct to an XML-based Config object
 	Config* cfg = structToConfig("Algorithm", prhs[1]);
 
-	CAlgorithm* pAlg = CAlgorithmFactory::getSingleton().create(cfg->self->getAttribute("type"));
+	CAlgorithm* pAlg = CAlgorithmFactory::getSingleton().create(cfg->self.getAttribute("type"));
 	if (!pAlg) {
 		delete cfg;
 		mexErrMsgTxt("Unknown algorithm. \n");
diff --git a/matlab/mex/astra_mex_c.cpp b/matlab/mex/astra_mex_c.cpp
index 4a331f5..a9b9654 100644
--- a/matlab/mex/astra_mex_c.cpp
+++ b/matlab/mex/astra_mex_c.cpp
@@ -36,9 +36,9 @@ $Id$
 #include "mexInitFunctions.h"
 
 #include "astra/Globals.h"
-
+#ifdef ASTRA_CUDA
 #include "../cuda/2d/darthelper.h"
-
+#endif
 using namespace std;
 using namespace astra;
 
diff --git a/matlab/mex/astra_mex_data2d_c.cpp b/matlab/mex/astra_mex_data2d_c.cpp
index 9576896..909d229 100644
--- a/matlab/mex/astra_mex_data2d_c.cpp
+++ b/matlab/mex/astra_mex_data2d_c.cpp
@@ -149,7 +149,7 @@ void astra_mex_data2d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArra
 		Config* cfg = structToConfig("ProjectionGeometry", prhs[2]);
 		// FIXME: Change how the base class is created. (This is duplicated
 		// in 'change_geometry' and Projector2D.cpp.)
-		std::string type = cfg->self->getAttribute("type");
+		std::string type = cfg->self.getAttribute("type");
 		CProjectionGeometry2D* pGeometry;
 		if (type == "sparse_matrix") {
 			pGeometry = new CSparseMatrixProjectionGeometry2D();
@@ -438,7 +438,7 @@ void astra_mex_data2d_change_geometry(int nlhs, mxArray* plhs[], int nrhs, const
 		Config* cfg = structToConfig("ProjectionGeometry2D", prhs[2]);
 		// FIXME: Change how the base class is created. (This is duplicated
 		// in 'create' and Projector2D.cpp.)
-		std::string type = cfg->self->getAttribute("type");
+		std::string type = cfg->self.getAttribute("type");
 		CProjectionGeometry2D* pGeometry;
 		if (type == "sparse_matrix") {
 			pGeometry = new CSparseMatrixProjectionGeometry2D();
diff --git a/matlab/mex/astra_mex_data3d_c.cpp b/matlab/mex/astra_mex_data3d_c.cpp
index 6096adc..fe4ce37 100644
--- a/matlab/mex/astra_mex_data3d_c.cpp
+++ b/matlab/mex/astra_mex_data3d_c.cpp
@@ -346,7 +346,7 @@ void astra_mex_data3d_change_geometry(int nlhs, mxArray* plhs[], int nrhs, const
 		astra::Config* cfg = structToConfig("ProjectionGeometry3D", geometry);
 		// FIXME: Change how the base class is created. (This is duplicated
 		// in Projector3D.cpp.)
-		std::string type = cfg->self->getAttribute("type");
+		std::string type = cfg->self.getAttribute("type");
 		astra::CProjectionGeometry3D* pGeometry = 0;
 		if (type == "parallel3d") {
 			pGeometry = new astra::CParallelProjectionGeometry3D();
diff --git a/matlab/mex/astra_mex_log_c.cpp b/matlab/mex/astra_mex_log_c.cpp
index ea4621e..905612c 100644
--- a/matlab/mex/astra_mex_log_c.cpp
+++ b/matlab/mex/astra_mex_log_c.cpp
@@ -55,7 +55,7 @@ void astra_mex_log_debug(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prh
 	string filename = mexToString(prhs[1]);
 	int linenumber = (int)mxGetScalar(prhs[2]);
 	string message = mexToString(prhs[3]);
-	astra::CLogger::debug(filename.c_str(),linenumber,message.c_str());
+	astra::CLogger::debug(filename.c_str(),linenumber,"%s",message.c_str());
 }
 
 //-----------------------------------------------------------------------------------------
@@ -75,7 +75,7 @@ void astra_mex_log_info(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs
 	string filename = mexToString(prhs[1]);
 	int linenumber = (int)mxGetScalar(prhs[2]);
 	string message = mexToString(prhs[3]);
-	astra::CLogger::info(filename.c_str(),linenumber,message.c_str());
+	astra::CLogger::info(filename.c_str(),linenumber,"%s",message.c_str());
 }
 
 //-----------------------------------------------------------------------------------------
@@ -95,7 +95,7 @@ void astra_mex_log_warn(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs
 	string filename = mexToString(prhs[1]);
 	int linenumber = (int)mxGetScalar(prhs[2]);
 	string message = mexToString(prhs[3]);
-	astra::CLogger::warn(filename.c_str(),linenumber,message.c_str());
+	astra::CLogger::warn(filename.c_str(),linenumber,"%s",message.c_str());
 }
 
 //-----------------------------------------------------------------------------------------
@@ -115,7 +115,7 @@ void astra_mex_log_error(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prh
 	string filename = mexToString(prhs[1]);
 	int linenumber = (int)mxGetScalar(prhs[2]);
 	string message = mexToString(prhs[3]);
-	astra::CLogger::error(filename.c_str(),linenumber,message.c_str());
+	astra::CLogger::error(filename.c_str(),linenumber,"%s",message.c_str());
 }
 
 //-----------------------------------------------------------------------------------------
diff --git a/matlab/mex/astra_mex_projector3d_c.cpp b/matlab/mex/astra_mex_projector3d_c.cpp
index c3b547f..e25802c 100644
--- a/matlab/mex/astra_mex_projector3d_c.cpp
+++ b/matlab/mex/astra_mex_projector3d_c.cpp
@@ -137,7 +137,9 @@ void astra_mex_projector3d_get_projection_geometry(int nlhs, mxArray* plhs[], in
 
 	// step3: get projection_geometry and turn it into a MATLAB struct
 	if (1 <= nlhs) {
-		plhs[0] = configToStruct(pProjector->getProjectionGeometry()->getConfiguration());
+		Config *cfg = pProjector->getProjectionGeometry()->getConfiguration();
+		plhs[0] = configToStruct(cfg);
+		delete cfg;
 	}
 }
 
@@ -163,7 +165,9 @@ void astra_mex_projector3d_get_volume_geometry(int nlhs, mxArray* plhs[], int nr
 
 	// step3: get projection_geometry and turn it into a MATLAB struct
 	if (1 <= nlhs) {
-		plhs[0] = configToStruct(pProjector->getVolumeGeometry()->getConfiguration());
+		Config *cfg = pProjector->getVolumeGeometry()->getConfiguration();
+		plhs[0] = configToStruct(cfg);
+		delete cfg;
 	}
 }
 
diff --git a/matlab/mex/astra_mex_projector_c.cpp b/matlab/mex/astra_mex_projector_c.cpp
index 204ba8e..bf701af 100644
--- a/matlab/mex/astra_mex_projector_c.cpp
+++ b/matlab/mex/astra_mex_projector_c.cpp
@@ -160,7 +160,9 @@ void astra_mex_projector_projection_geometry(int nlhs, mxArray* plhs[], int nrhs
 
 	// step3: get projection_geometry and turn it into a MATLAB struct
 	if (1 <= nlhs) {
-		plhs[0] = configToStruct(pProjector->getProjectionGeometry()->getConfiguration());
+		Config *cfg =  pProjector->getProjectionGeometry()->getConfiguration();
+		plhs[0] = configToStruct(cfg);
+		delete cfg;
 	}
 }
 
@@ -189,8 +191,9 @@ void astra_mex_projector_volume_geometry(int nlhs, mxArray* plhs[], int nrhs, co
 
 	// step3: get projection_geometry and turn it into a MATLAB struct
 	if (1 <= nlhs) {
-		plhs[0] = configToStruct(pProjector->getVolumeGeometry()->getConfiguration());
-
+		Config *cfg = pProjector->getVolumeGeometry()->getConfiguration();
+		plhs[0] = configToStruct(cfg);
+		delete cfg;
 	}
 }
 
diff --git a/matlab/mex/mexDataManagerHelpFunctions.cpp b/matlab/mex/mexDataManagerHelpFunctions.cpp
index d482428..1794abb 100644
--- a/matlab/mex/mexDataManagerHelpFunctions.cpp
+++ b/matlab/mex/mexDataManagerHelpFunctions.cpp
@@ -285,7 +285,7 @@ allocateDataObject(const std::string & sDataType,
 		astra::Config* cfg = structToConfig("ProjectionGeometry3D", geometry);
 		// FIXME: Change how the base class is created. (This is duplicated
 		// in Projector3D.cpp.)
-		std::string type = cfg->self->getAttribute("type");
+		std::string type = cfg->self.getAttribute("type");
 		astra::CProjectionGeometry3D* pGeometry = 0;
 		if (type == "parallel3d") {
 			pGeometry = new astra::CParallelProjectionGeometry3D();
diff --git a/matlab/mex/mexHelpFunctions.cpp b/matlab/mex/mexHelpFunctions.cpp
index c0ac711..58e84d2 100644
--- a/matlab/mex/mexHelpFunctions.cpp
+++ b/matlab/mex/mexHelpFunctions.cpp
@@ -185,7 +185,7 @@ Config* structToConfig(string rootname, const mxArray* pStruct)
 }
 
 //-----------------------------------------------------------------------------------------
-bool structToXMLNode(XMLNode* node, const mxArray* pStruct) 
+bool structToXMLNode(XMLNode node, const mxArray* pStruct) 
 {
 	// loop all fields
 	int nfields = mxGetNumberOfFields(pStruct);
@@ -199,16 +199,16 @@ bool structToXMLNode(XMLNode* node, const mxArray* pStruct)
 		if (mxIsChar(pField)) {
 			string sValue = mexToString(pField);
 			if (sFieldName == "type") {
-				node->addAttribute("type", sValue);
+				node.addAttribute("type", sValue);
 			} else {
-				delete node->addChildNode(sFieldName, sValue);
+				node.addChildNode(sFieldName, sValue);
 			}
 		}
 
 		// scalar
 		else if (mxIsNumeric(pField) && mxGetM(pField)*mxGetN(pField) == 1) {
 			string sValue = mexToString(pField);
-			delete node->addChildNode(sFieldName, sValue);
+			node.addChildNode(sFieldName, sValue);
 		}
 
 		// numerical array
@@ -217,20 +217,9 @@ bool structToXMLNode(XMLNode* node, const mxArray* pStruct)
 				mexErrMsgTxt("Numeric input must be double.");
 				return false;
 			}
-			XMLNode* listbase = node->addChildNode(sFieldName);
-			listbase->addAttribute("listsize", mxGetM(pField)*mxGetN(pField));
+			XMLNode listbase = node.addChildNode(sFieldName);
 			double* pdValues = mxGetPr(pField);
-			int index = 0;
-			for (unsigned int row = 0; row < mxGetM(pField); row++) {
-				for (unsigned int col = 0; col < mxGetN(pField); col++) {
-					XMLNode* item = listbase->addChildNode("ListItem");
-					item->addAttribute("index", index);
-					item->addAttribute("value", pdValues[col*mxGetM(pField)+row]);
-					index++;
-					delete item;
-				}
-			}
-			delete listbase;
+			listbase.setContent(pdValues, mxGetN(pField), mxGetM(pField), true);
 		}
 
 		// not castable to a single string
@@ -240,9 +229,8 @@ bool structToXMLNode(XMLNode* node, const mxArray* pStruct)
 				if (!ret)
 					return false;
 			} else {
-				XMLNode* newNode = node->addChildNode(sFieldName);
+				XMLNode newNode = node.addChildNode(sFieldName);
 				bool ret = structToXMLNode(newNode, pField);
-				delete newNode;
 				if (!ret)
 					return false;
 			}
@@ -254,7 +242,7 @@ bool structToXMLNode(XMLNode* node, const mxArray* pStruct)
 }
 //-----------------------------------------------------------------------------------------
 // Options struct to xml node
-bool optionsToXMLNode(XMLNode* node, const mxArray* pOptionStruct)
+bool optionsToXMLNode(XMLNode node, const mxArray* pOptionStruct)
 {
 	// loop all fields
 	int nfields = mxGetNumberOfFields(pOptionStruct);
@@ -262,7 +250,7 @@ bool optionsToXMLNode(XMLNode* node, const mxArray* pOptionStruct)
 		std::string sFieldName = std::string(mxGetFieldNameByNumber(pOptionStruct, i));
 		const mxArray* pField = mxGetFieldByNumber(pOptionStruct, 0, i);
 
-		if (node->hasOption(sFieldName)) {
+		if (node.hasOption(sFieldName)) {
 			mexErrMsgTxt("Duplicate option");
 			return false;
 		}
@@ -270,7 +258,7 @@ bool optionsToXMLNode(XMLNode* node, const mxArray* pOptionStruct)
 		// string or scalar
 		if (mxIsChar(pField) || mexIsScalar(pField)) {
 			string sValue = mexToString(pField);
-			node->addOption(sFieldName, sValue);
+			node.addOption(sFieldName, sValue);
 		}
 		// numerical array
 		else if (mxIsNumeric(pField) && mxGetM(pField)*mxGetN(pField) > 1) {
@@ -279,21 +267,10 @@ bool optionsToXMLNode(XMLNode* node, const mxArray* pOptionStruct)
 				return false;
 			}
 
-			XMLNode* listbase = node->addChildNode("Option");
-			listbase->addAttribute("key", sFieldName);
-			listbase->addAttribute("listsize", mxGetM(pField)*mxGetN(pField));
+			XMLNode listbase = node.addChildNode("Option");
+			listbase.addAttribute("key", sFieldName);
 			double* pdValues = mxGetPr(pField);
-			int index = 0;
-			for (unsigned int row = 0; row < mxGetM(pField); row++) {
-				for (unsigned int col = 0; col < mxGetN(pField); col++) {
-					XMLNode* item = listbase->addChildNode("ListItem");
-					item->addAttribute("index", index);
-					item->addAttribute("value", pdValues[col*mxGetM(pField)+row]);
-					index++;
-					delete item;
-				}
-			}
-			delete listbase;
+			listbase.setContent(pdValues, mxGetN(pField), mxGetM(pField), true);
 		} else {
 			mexErrMsgTxt("Unsupported option type");
 			return false;
@@ -343,30 +320,33 @@ mxArray* configToStruct(astra::Config* cfg)
 }
 
 //-----------------------------------------------------------------------------------------
-mxArray* XMLNodeToStruct(astra::XMLNode* node)
+mxArray* XMLNodeToStruct(astra::XMLNode node)
 {
 	std::map<std::string, mxArray*> mList;
 	std::map<std::string, mxArray*> mOptions;
 
 	// type_attribute
-	if (node->hasAttribute("type")) {
-		mList["type"] = mxCreateString(node->getAttribute("type").c_str());
+	if (node.hasAttribute("type")) {
+		mList["type"] = mxCreateString(node.getAttribute("type").c_str());
 	}
 
-	list<XMLNode*> nodes = node->getNodes();
-	for (list<XMLNode*>::iterator it = nodes.begin(); it != nodes.end(); it++) {
-		XMLNode* subnode = (*it);
+	list<XMLNode> nodes = node.getNodes();
+	for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); it++) {
+		XMLNode subnode = (*it);
 
 		// option
-		if (subnode->getName() == "Option") {
-			mOptions[subnode->getAttribute("key")] = stringToMxArray(subnode->getAttribute("value"));
+		if (subnode.getName() == "Option") {
+			if(subnode.hasAttribute("value")){
+				mOptions[subnode.getAttribute("key")] = stringToMxArray(subnode.getAttribute("value"));
+			}else{
+				mOptions[subnode.getAttribute("key")] = stringToMxArray(subnode.getContent());
+			}
 		}
 
 		// regular content
 		else {
-			mList[subnode->getName()] = stringToMxArray(subnode->getContent());
+			mList[subnode.getName()] = stringToMxArray(subnode.getContent());
 		}
-		delete subnode;
 	}
 
 	if (mOptions.size() > 0) mList["options"] = buildStruct(mOptions);
diff --git a/matlab/mex/mexHelpFunctions.h b/matlab/mex/mexHelpFunctions.h
index f9ffcf2..3ac5bd8 100644
--- a/matlab/mex/mexHelpFunctions.h
+++ b/matlab/mex/mexHelpFunctions.h
@@ -58,13 +58,13 @@ mxArray* anyToMxArray(boost::any _any);
 
 // turn a MATLAB struct into a Config object
 astra::Config* structToConfig(string rootname, const mxArray* pStruct);
-bool structToXMLNode(astra::XMLNode* node, const mxArray* pStruct);
-bool optionsToXMLNode(astra::XMLNode* node, const mxArray* pOptionStruct);
+bool structToXMLNode(astra::XMLNode node, const mxArray* pStruct);
+bool optionsToXMLNode(astra::XMLNode node, const mxArray* pOptionStruct);
 std::map<std::string, mxArray*> parseStruct(const mxArray* pInput);
 
 // turn a Config object into a MATLAB struct
 mxArray* configToStruct(astra::Config* cfg);
-mxArray* XMLNodeToStruct(astra::XMLNode* xml);
+mxArray* XMLNodeToStruct(astra::XMLNode xml);
 mxArray* stringToMxArray(std::string input);
 mxArray* buildStruct(std::map<std::string, mxArray*> mInput);
 
diff --git a/matlab/mex/mexInitFunctions.cpp b/matlab/mex/mexInitFunctions.cpp
index d8a50d7..89a31a1 100644
--- a/matlab/mex/mexInitFunctions.cpp
+++ b/matlab/mex/mexInitFunctions.cpp
@@ -8,7 +8,7 @@ bool mexIsInitialized=false;
  *
  */
 void logCallBack(const char *msg, size_t len){
-    mexPrintf(msg);
+    mexPrintf("%s",msg);
 }
 
 /**
diff --git a/matlab/tools/astra_create_fbp_reconstruction.m b/matlab/tools/astra_create_fbp_reconstruction.m
index 5540f27..a2561b7 100644
--- a/matlab/tools/astra_create_fbp_reconstruction.m
+++ b/matlab/tools/astra_create_fbp_reconstruction.m
@@ -19,6 +19,7 @@ cfg.ProjectorId = proj_id;
 cfg.Options.GPUindex = 0;
 alg_id = astra_mex_algorithm('create', cfg);
 astra_mex_algorithm('run', alg_id);
+astra_mex_algorithm('delete', alg_id);
 
 if numel(sinogram) ~= 1
 	astra_mex_data2d('delete', sinogram_id);
diff --git a/matlab/tools/astra_create_reconstruction_cuda.m b/matlab/tools/astra_create_reconstruction_cuda.m
index 7d9e1dd..7d0421c 100644
--- a/matlab/tools/astra_create_reconstruction_cuda.m
+++ b/matlab/tools/astra_create_reconstruction_cuda.m
@@ -45,7 +45,7 @@ if strcmp(rec_type,'')
 end
 
 % configure
-cfg = astra_struct('SIRT_CUDA');
+cfg = astra_struct(rec_type);
 cfg.ProjectionGeometry = proj_geom;
 cfg.ReconstructionGeometry = vol_geom;
 cfg.ProjectionDataId = sinogram_id;
diff --git a/matlab/tools/opTomo.m b/matlab/tools/opTomo.m
new file mode 100644
index 0000000..14128d2
--- /dev/null
+++ b/matlab/tools/opTomo.m
@@ -0,0 +1,280 @@
+%OPTOMO Wrapper for ASTRA tomography projector
+%
+%   OP = OPTOMO(TYPE, PROJ_GEOM, VOL_GEOM) generates a Spot operator OP for
+%   the ASTRA forward and backprojection operations. The string TYPE
+%   determines the model used for the projections. Possible choices are:
+%       TYPE:  * using the CPU
+%                'line'   - use a line kernel
+%                'linear' - use a Joseph kernel
+%                'strip'  - use the strip kernel
+%              * using the GPU
+%                'cuda'  - use a Joseph kernel, on the GPU, currently using
+%                          'cuda' is the only option in 3D.
+%   The PROJ_GEOM and VOL_GEOM structures are projection and volume
+%   geometries as used in the ASTRA toolbox.
+%
+%   OP = OPTOMO(TYPE, PROJ_GEOM, VOL_GEOM, GPU_INDEX) also specify the
+%   index of the GPU that should be used, if multiple GPUs are present in
+%   the host system. By default GPU_INDEX is 0.
+%
+%   Note: this code depends on the Matlab toolbox 
+%   "Spot - A Linear-Operator Toolbox" which can be downloaded from
+%   http://www.cs.ubc.ca/labs/scl/spot/
+%--------------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+% 
+% Copyright: 2014-2015, CWI, Amsterdam
+% License:   Open Source under GPLv3
+% Author:    Folkert Bleichrodt
+% Contact:   F.Bleichrodt@cwi.nl
+% Website:   http://sf.net/projects/astra-toolbox
+%--------------------------------------------------------------------------
+% $Id$
+
+classdef opTomo < opSpot
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Properties
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    properties ( Access = private )
+        % multiplication function
+        funHandle
+        % ASTRA identifiers
+        sino_id
+        vol_id
+        fp_alg_id
+        bp_alg_id
+        % ASTRA IDs handle
+        astra_handle
+        % geometries
+        proj_geom;
+        vol_geom;
+    end % properties
+    
+    properties ( SetAccess = private, GetAccess = public )
+        proj_size
+        vol_size
+    end % properties
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Methods - public
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    methods
+
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        % opTomo - constructor
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        function op = opTomo(type, proj_geom, vol_geom, gpu_index)
+            
+            if nargin < 4 || isempty(gpu_index), gpu_index = 0; end
+            
+            proj_size = astra_geom_size(proj_geom);
+            vol_size  = astra_geom_size(vol_geom);
+            
+            % construct operator
+            op = op@opSpot('opTomo', prod(proj_size), prod(vol_size));
+            
+            % determine the dimension
+            is2D = ~isfield(vol_geom, 'GridSliceCount');
+            gpuEnabled = strcmpi(type, 'cuda');
+            
+            if is2D
+                % create a projector
+                proj_id = astra_create_projector(type, proj_geom, vol_geom);
+                
+                % create a function handle
+                op.funHandle = @opTomo_intrnl2D;
+                
+                % Initialize ASTRA data objects.
+                % projection data
+                sino_id = astra_mex_data2d('create', '-sino', proj_geom, 0);
+                % image data
+                vol_id  = astra_mex_data2d('create', '-vol', vol_geom, 0);
+                
+                % Setup forward and back projection algorithms.
+                if gpuEnabled
+                    fp_alg = 'FP_CUDA';
+                    bp_alg = 'BP_CUDA';
+                    proj_id = [];
+                else
+                    fp_alg = 'FP';
+                    bp_alg = 'BP';
+                    proj_id = astra_create_projector(type, proj_geom, vol_geom);
+                end
+                
+                % configuration for ASTRA fp algorithm
+                cfg_fp = astra_struct(fp_alg);
+                cfg_fp.ProjectorId      = proj_id;
+                cfg_fp.ProjectionDataId = sino_id;
+                cfg_fp.VolumeDataId     = vol_id;
+                
+                % configuration for ASTRA bp algorithm
+                cfg_bp = astra_struct(bp_alg);
+                cfg_bp.ProjectionDataId     = sino_id;
+                cfg_bp.ProjectorId          = proj_id;
+                cfg_bp.ReconstructionDataId = vol_id;
+                
+                % set GPU index
+                if gpuEnabled
+                    cfg_fp.option.GPUindex = gpu_index;
+                    cfg_bp.option.GPUindex = gpu_index;
+                end
+                
+                fp_alg_id = astra_mex_algorithm('create', cfg_fp);
+                bp_alg_id = astra_mex_algorithm('create', cfg_bp);
+                
+                % Create handle to ASTRA objects, so they will be deleted
+                % if opTomo is deleted.
+                op.astra_handle = opTomo_helper_handle([sino_id, ...
+                    vol_id, proj_id, fp_alg_id, bp_alg_id]);
+                
+                op.fp_alg_id   = fp_alg_id;
+                op.bp_alg_id   = bp_alg_id;
+                op.sino_id     = sino_id;
+                op.vol_id      = vol_id;
+            else
+                % 3D
+                % only gpu/cuda code for 3D
+                if ~gpuEnabled
+                    error(['Only type ' 39 'cuda' 39 ' is supported ' ...
+                           'for 3D geometries.'])
+                end
+                
+                % create a function handle
+                op.funHandle = @opTomo_intrnl3D;
+            end
+      
+            
+            % pass object properties
+            op.proj_size   = proj_size;
+            op.vol_size    = vol_size;
+            op.proj_geom   = proj_geom;
+            op.vol_geom    = vol_geom;
+            op.cflag       = false;
+            op.sweepflag   = false;
+
+        end % opTomo - constructor
+        
+    end % methods - public
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Methods - protected
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    methods( Access = protected )
+
+        % multiplication
+        function y = multiply(op,x,mode)
+            
+            % ASTRA cannot handle sparse vectors
+            if issparse(x)
+                x = full(x);
+            end
+            
+            % convert input to single
+            if isa(x, 'single') == false
+                x = single(x);
+            end
+            
+            % the multiplication
+            y = op.funHandle(op, x, mode);
+            
+            % make sure output is column vector
+            y = y(:);
+            
+        end % multiply
+        
+    end % methods - protected
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Methods - private
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    methods( Access = private )
+        
+        % 2D projection code
+        function y = opTomo_intrnl2D(op,x,mode)
+                       
+            if mode == 1              
+                % X is passed as a vector, reshape it into an image.             
+                x = reshape(x, op.vol_size);
+                
+                % Matlab data copied to ASTRA data
+                astra_mex_data2d('store', op.vol_id, x);
+                
+                % forward projection
+                astra_mex_algorithm('iterate', op.fp_alg_id);
+                
+                % retrieve Matlab array
+                y = astra_mex_data2d('get_single', op.sino_id);
+            else
+                % X is passed as a vector, reshape it into a sinogram.
+                x = reshape(x, op.proj_size);
+                
+                % Matlab data copied to ASTRA data
+                astra_mex_data2d('store', op.sino_id, x);
+                
+                % backprojection
+                astra_mex_algorithm('iterate', op.bp_alg_id);
+                
+                % retrieve Matlab array
+                y = astra_mex_data2d('get_single', op.vol_id);
+            end
+        end % opTomo_intrnl2D
+        
+        
+        % 3D projection code
+        function y = opTomo_intrnl3D(op,x,mode)
+            
+            if mode == 1
+                % X is passed as a vector, reshape it into an image
+                x = reshape(x, op.vol_size);
+                
+                % initialize output
+                y = zeros(op.proj_size, 'single');
+                
+                % link matlab array to ASTRA
+                vol_id  = astra_mex_data3d_c('link', '-vol', op.vol_geom, x, 0);
+                sino_id = astra_mex_data3d_c('link', '-sino', op.proj_geom, y, 1);
+                
+                % initialize fp algorithm
+                cfg = astra_struct('FP3D_CUDA');
+                cfg.ProjectionDataId = sino_id;
+                cfg.VolumeDataId     = vol_id;
+                
+                alg_id = astra_mex_algorithm('create', cfg);
+                             
+                % forward projection
+                astra_mex_algorithm('iterate', alg_id);
+                
+                % cleanup
+                astra_mex_data3d('delete', vol_id);
+                astra_mex_data3d('delete', sino_id);
+            else
+                % X is passed as a vector, reshape it into projection data
+                x = reshape(x, op.proj_size);
+                
+                % initialize output
+                y = zeros(op.vol_size,'single');
+                
+                % link matlab array to ASTRA
+                vol_id  = astra_mex_data3d_c('link', '-vol', op.vol_geom, y, 1);
+                sino_id = astra_mex_data3d_c('link', '-sino', op.proj_geom, x, 0);
+                                
+                % initialize bp algorithm
+                cfg = astra_struct('BP3D_CUDA');
+                cfg.ProjectionDataId     = sino_id;
+                cfg.ReconstructionDataId = vol_id;
+                
+                alg_id = astra_mex_algorithm('create', cfg);
+
+                % backprojection
+                astra_mex_algorithm('iterate', alg_id);
+                
+                % cleanup
+                astra_mex_data3d('delete', vol_id);
+                astra_mex_data3d('delete', sino_id);
+            end 
+        end % opTomo_intrnl3D
+        
+    end % methods - private
+ 
+end % classdef
diff --git a/matlab/tools/opTomo_helper_handle.m b/matlab/tools/opTomo_helper_handle.m
new file mode 100644
index 0000000..d9be51f
--- /dev/null
+++ b/matlab/tools/opTomo_helper_handle.m
@@ -0,0 +1,29 @@
+classdef opTomo_helper_handle < handle
+    %ASTRA.OPTOMO_HELPER_HANDLE Handle class around an astra identifier
+    %   Automatically deletes the data when object is deleted. 
+    %   Multiple id's can be passed as an array as input to 
+    %   the constructor.
+    
+    properties
+        id
+    end
+    
+    methods
+        function obj = opTomo_helper_handle(id)
+            obj.id = id;
+        end
+        function delete(obj)
+            for i = 1:numel(obj.id)
+                % delete any kind of object
+                astra_mex_data2d('delete', obj.id(i));
+                astra_mex_data3d('delete', obj.id(i));
+                astra_mex_algorithm('delete', obj.id(i));
+                astra_mex_matrix('delete', obj.id(i));
+                astra_mex_projector('delete', obj.id(i));
+                astra_mex_projector3d('delete', obj.id(i))
+            end
+        end
+    end
+    
+end
+