summaryrefslogtreecommitdiffstats
path: root/matlab
diff options
context:
space:
mode:
Diffstat (limited to 'matlab')
-rw-r--r--matlab/mex/astra_mex_algorithm_c.cpp2
-rw-r--r--matlab/mex/astra_mex_data2d_c.cpp4
-rw-r--r--matlab/mex/astra_mex_data3d_c.cpp2
-rw-r--r--matlab/mex/mexDataManagerHelpFunctions.cpp2
-rw-r--r--matlab/mex/mexHelpFunctions.cpp68
-rw-r--r--matlab/mex/mexHelpFunctions.h6
-rw-r--r--matlab/tools/opTomo.m280
-rw-r--r--matlab/tools/opTomo_helper_handle.m29
8 files changed, 339 insertions, 54 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_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/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..87a9672 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,29 @@ 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") {
+ mOptions[subnode.getAttribute("key")] = stringToMxArray(subnode.getAttribute("value"));
}
// 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/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
+
generated by cgit v1.2.3 (git 2.25.1) at 2024-09-21 18:02:42 +0000