Update version to 1.3

18 files changed, 2473 insertions, 0 deletions

diff --git a/matlab/algorithms/DART/DARTalgorithm.m b/matlab/algorithms/DART/DARTalgorithm.m
new file mode 100644
index 0000000..b7e63b5
--- /dev/null
+++ b/matlab/algorithms/DART/DARTalgorithm.m
@@ -0,0 +1,229 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef DARTalgorithm < matlab.mixin.Copyable
+
+	% Algorithm class for Discrete Algebraic Reconstruction Technique (DART).
+	
+	% todo: reset()
+	% todo: fixed random seed
+	% todo: initialize from settings (?)
+
+	%----------------------------------------------------------------------
+	properties (GetAccess=public, SetAccess=public)
+	
+		tomography		= TomographyDefault();		% POLICY: Tomography object. 
+		segmentation	= SegmentationDefault();	% POLICY: Segmentation object. 
+		smoothing		= SmoothingDefault();		% POLICY: Smoothing object.
+		masking			= MaskingDefault();			% POLICY: Masking object.
+		output			= OutputDefault();			% POLICY: Output object.
+		statistics		= StatisticsDefault();		% POLICY: Statistics object.
+
+		base			= struct();					% DATA(set): base structure, should contain: 'sinogram', 'proj_geom', 'phantom' (optional).
+		
+		memory			= 'yes';					% SETTING: reduce memory usage? (disables some features)
+		
+		testdata = struct();
+		
+	end
+	%----------------------------------------------------------------------
+	properties (GetAccess=public, SetAccess=private)
+		V0					= [];			% DATA(get): Initial reconstruction.
+		V					= [];			% DATA(get): Reconstruction.
+		S					= [];			% DATA(get): Segmentation.
+		R					= [];			% DATA(get): Residual projection data.
+		Mask				= [];			% DATA(get): Reconstruction Mask.
+		stats				= struct();		% Structure containing various statistics.
+		iterationcount		= 0;			% Number of performed iterations.	
+		start_tic			= 0;		  
+		initialized			= 0;			% Is initialized?
+	end
+	%----------------------------------------------------------------------
+	properties (Access=private)		
+		adaptparam_name		= {};
+		adaptparam_values	= {};
+		adaptparam_iters	= {};
+	end
+	
+	%----------------------------------------------------------------------
+	methods
+		
+		%------------------------------------------------------------------
+		function this = DARTalgorithm(base)
+			% Constructor
+			% >> D = DARTalgorithm(base); base is a matlab struct (or the path towards one) 
+			%					          that should contain 'sinogram', 'proj_geom'.
+
+			if ischar(base)
+				this.base = load(base);
+			else
+				this.base = base;
+			end
+		end	
+	
+		
+		%------------------------------------------------------------------
+		function D = deepcopy(this)
+			% Create a deep copy of this object.
+			% >> D2 = D.deepcopy();
+			
+			D = copy(this);
+			props = properties(this);
+			for i = 1:length(props)
+				if isa(this.(props{i}), 'handle')
+					D.(props{i}) = copy(this.(props{i}));
+				end
+			end			
+	
+		end
+		
+		%------------------------------------------------------------------
+		function this = initialize(this)
+			% Initializes this object.
+			% >> D.initialize();
+			
+			% Initialize tomography part
+			this.tomography.initialize(this);
+
+			% Create an Initial Reconstruction
+			if isfield(this.base, 'V0')
+				this.V0 = this.base.V0;
+			else
+				this.output.pre_initial_iteration(this);
+				this.V0 = this.tomography.createInitialReconstruction(this, this.base.sinogram);
+				this.output.post_initial_iteration(this);
+			end
+			this.V = this.V0;
+			if strcmp(this.memory,'yes')
+				this.base.V0 = [];
+				this.V0 = [];
+			end
+			this.initialized = 1;
+		end		
+		
+		%------------------------------------------------------------------
+		% iterate
+		function this = iterate(this, iters)
+			% Perform several iterations of the DART algorithm.
+			% >> D.iterate(iterations);			
+			
+			this.start_tic = tic;
+			
+			for iteration = 1:iters
+				
+				this.iterationcount = this.iterationcount + 1;				
+		
+				%----------------------------------------------------------
+				% Initial Output
+				this.output.pre_iteration(this);
+
+				%----------------------------------------------------------
+				% Update Adaptive Parameters
+				for i = 1:numel(this.adaptparam_name)
+					
+					for j = 1:numel(this.adaptparam_iters{i})
+						if this.iterationcount == this.adaptparam_iters{i}(j)
+							new_value = this.adaptparam_values{i}(j);
+							eval(['this.' this.adaptparam_name{i} ' = ' num2str(new_value) ';']);
+							disp(['value ' this.adaptparam_name{i} ' changed to ' num2str(new_value) ]);
+						end						
+					end
+				
+				end
+
+				%----------------------------------------------------------
+				% Segmentation
+				this.segmentation.estimate_grey_levels(this, this.V);
+				this.S = this.segmentation.apply(this, this.V);
+
+				%----------------------------------------------------------
+				% Select Update and Fixed Pixels
+				this.Mask = this.masking.apply(this, this.S);
+
+				this.V = (this.V .* this.Mask) + (this.S .* (1 - this.Mask));
+				%this.V(this.Mask == 0) = this.S(this.Mask == 0); 
+
+				F = this.V;
+				F(this.Mask == 1) = 0;
+
+				%----------------------------------------------------------
+				% Create Residual Projection Difference
+				%this.testdata.F{iteration} = F;
+				this.R = this.base.sinogram - this.tomography.createForwardProjection(this, F);
+				%this.testdata.R{iteration} = this.R;
+				
+				%----------------------------------------------------------
+				% ART Loose Part
+				%this.testdata.V1{iteration} = this.V;
+				%this.testdata.Mask{iteration} = this.Mask;
+				
+				%X = zeros(size(this.V));
+				%Y = this.tomography.createReconstruction(this, this.R, X, this.Mask);
+				%this.V(this.Mask) = Y(this.Mask);
+				this.V = this.tomography.createReconstruction(this, this.R, this.V, this.Mask);
+				
+				%this.testdata.V2{iteration} = this.V;
+				
+				%----------------------------------------------------------
+				% Blur
+				this.V = this.smoothing.apply(this, this.V);
+				%this.testdata.V3{iteration} = this.V;
+				
+				%----------------------------------------------------------
+				% Calculate Statistics
+				this.stats = this.statistics.apply(this);
+
+				%----------------------------------------------------------
+				% Output
+				this.output.post_iteration(this);
+			
+			end % end iteration loop
+			
+			%test = this.testdata;
+			%save('testdata.mat','test');
+		
+		end				
+
+		%------------------------------------------------------------------
+		% get data		
+		function data = getdata(this, string)
+			if numel(this.(string)) == 1
+				data = astra_mex_data2d('get',this.(string));
+			else
+				data = this.(string);
+			end
+		end
+		
+		%------------------------------------------------------------------
+		% add adaptive parameter
+		function this = adaptiveparameter(this, name, values, iterations)
+			this.adaptparam_name{end+1} = name;
+			this.adaptparam_values{end+1} = values;
+			this.adaptparam_iters{end+1} = iterations;
+		end	
+		
+		%------------------------------------------------------------------
+		function settings = getsettings(this)
+			% Returns a structure containing all settings of this object.
+			% >> settings = tomography.getsettings();		
+			
+			settings.tomography = this.tomography.getsettings();
+			settings.smoothing = this.smoothing.getsettings();
+			settings.masking = this.masking.getsettings();
+			settings.segmentation = this.segmentation.getsettings();
+		end		
+		%------------------------------------------------------------------
+		
+	end % methods
+		
+end % class
+
+
diff --git a/matlab/algorithms/DART/IterativeTomography.m b/matlab/algorithms/DART/IterativeTomography.m
new file mode 100644
index 0000000..b30640e
--- /dev/null
+++ b/matlab/algorithms/DART/IterativeTomography.m
@@ -0,0 +1,455 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef IterativeTomography < matlab.mixin.Copyable
+
+	% Algorithm class for 2D Iterative Tomography.
+
+	%----------------------------------------------------------------------
+	properties (SetAccess=public, GetAccess=public)				
+		superresolution		= 1;				% SETTING: Volume upsampling factor.
+		proj_type			= 'linear';			% SETTING: Projector type, only when gpu='no'.
+		method				= 'SIRT_CUDA';		% SETTING: Iterative method (see ASTRA toolbox documentation).
+		gpu					= 'yes';			% SETTING: Use gpu? {'yes', 'no'}
+		gpu_core			= 0;				% SETTING: Which gpu core to use? Only when gpu='yes'.
+		inner_circle		= 'yes';			% SETTING: Do roi only? {'yes', 'no'}
+		image_size			= [];				% SETTING: Overwrite default reconstruction size.  Only if no vol_geom is specified.
+		use_minc			= 'no';				% SETTING: Use minimum constraint. {'no', 'yes'}
+	end
+	%----------------------------------------------------------------------
+	properties (SetAccess=public, GetAccess=public)	
+		sinogram			= [];				% DATA: Projection data.
+		proj_geom			= [];				% DATA: Projection geometry.
+		V					= [];				% DATA: Volume data.  Also used to set initial estimate (optional).
+		vol_geom			= [];				% DATA: Volume geometry.
+	end
+	%----------------------------------------------------------------------
+	properties (SetAccess=private, GetAccess=public)
+		initialized     	= 0;				% Is this object initialized?
+	end
+	%----------------------------------------------------------------------
+	properties (SetAccess=protected, GetAccess=protected)
+		proj_geom_sr		= [];				% PROTECTED: geometry of sinogram (with super-resolution)
+        proj_id             = [];				% PROTECTED: astra id of projector (when gpu='no')
+        proj_id_sr          = [];				% PROTECTED: astra id of super-resolution projector (when gpu='no')	
+		cfg_base			= struct();			% PROTECTED: base configuration structure for the reconstruction algorithm.
+    end
+	%----------------------------------------------------------------------
+
+	methods (Access=public)
+		
+		%------------------------------------------------------------------
+		function this = IterativeTomography(varargin)
+			% Constructor
+			% >> tomography = IterativeTomography(); 
+			% >> tomography = IterativeTomography(base); base struct should contain fields 'sinogram' and 'proj_geom'.
+			% >> tomography = IterativeTomography(base_filename); mat-file should contain 'sinogram' and 'proj_geom'.
+			% >> tomography = IterativeTomography(proj_geom, vol_geom);
+			% >> tomography = IterativeTomography(sinogram, proj_geom);
+			
+			% Input: IterativeTomography(base)
+			if nargin == 1
+				if ischar(varargin{1})
+					this.sinogram = load(varargin{1},'sinogram');
+					this.proj_geom = load(varargin{1},'proj_geom');				
+				else
+					this.sinogram = varargin{1}.sinogram;
+					this.proj_geom = varargin{1}.proj_geom;
+				end
+			end
+			% Input: IterativeTomography(proj_geom, vol_geom)
+			if nargin == 2 && isstruct(varargin{1}) && isstruct(varargin{2})
+				this.proj_geom = varargin{1};
+				this.vol_geom =	varargin{2};
+			% Input: IterativeTomography(sinogram, proj_geom)
+			elseif nargin == 2
+				this.sinogram = varargin{1};
+				this.proj_geom = varargin{2};
+			end		
+
+		end
+
+		%------------------------------------------------------------------
+		function delete(this)
+			% Destructor
+			% >> clear tomography;
+			if strcmp(this.gpu,'no') && numel(this.proj_id) > 0
+				astra_mex_projector('delete', this.proj_id, this.proj_id_sr);
+			end
+		end
+			
+		%------------------------------------------------------------------
+		function settings = getsettings(this)
+			% Returns a structure containing all settings of this object.
+			% >> settings = tomography.getsettings();
+			settings.superresolution	= this.superresolution;
+			settings.proj_type			= this.proj_type;
+			settings.method				= this.method;
+			settings.gpu				= this.gpu;
+			settings.gpu_core			= this.gpu_core;
+			settings.inner_circle		= this.inner_circle;
+			settings.image_size			= this.image_size;	
+			settings.use_minc			= this.use_minc;
+		end		
+		
+		%------------------------------------------------------------------
+		function ok = initialize(this)
+			% Initialize this object.  Returns 1 if succesful.
+			% >> tomography.initialize();
+			
+			% create projection geometry with super-resolution
+			this.proj_geom_sr = astra_geom_superresolution(this.proj_geom, this.superresolution);			
+			
+			% if no volume geometry is specified by the user: create volume geometry
+			if numel(this.vol_geom) == 0
+				if numel(this.image_size) < 2
+					this.image_size(1) = this.proj_geom.DetectorCount;
+					this.image_size(2) = this.proj_geom.DetectorCount;
+				end
+				this.vol_geom = astra_create_vol_geom(this.image_size(1) * this.superresolution, this.image_size(2) * this.superresolution, ...
+													 -this.image_size(1)/2, this.image_size(1)/2, -this.image_size(2)/2, this.image_size(2)/2);
+			else
+				this.image_size(1) = this.vol_geom.GridRowCount;
+				this.image_size(2) = this.vol_geom.GridColCount;
+			end
+			
+			% create projector
+			if strcmp(this.gpu, 'no')
+				this.proj_id = astra_create_projector(this.proj_type, this.proj_geom, this.vol_geom);
+				this.proj_id_sr = astra_create_projector(this.proj_type, this.proj_geom_sr, this.vol_geom);			
+			end
+			
+			% create reconstruction configuration
+			this.cfg_base = astra_struct(upper(this.method));
+			if strcmp(this.gpu,'no')
+				this.cfg_base.ProjectorId = this.proj_id;
+				this.cfg_base.ProjectionGeometry = this.proj_geom;
+				this.cfg_base.ReconstructionGeometry = this.vol_geom;
+				this.cfg_base.option.ProjectionOrder = 'random';
+			end
+			this.cfg_base.option.DetectorSuperSampling = this.superresolution;
+			%this.cfg_base.option.DetectorSuperSampling = 8;
+			if strcmp(this.gpu,'yes')
+				this.cfg_base.option.GPUindex = this.gpu_core;
+			end
+			this.cfg_base.option.UseMinConstraint = this.use_minc;
+		
+			this.initialized = 1;
+			ok = this.initialized;
+		end
+		
+		%------------------------------------------------------------------
+		function P = project(this, volume)
+			% Compute forward projection.  
+			% Stores sinogram in tomography.sinogram if it is still empty.
+			% >> P = tomography.project();   projects 'tomography.V'.
+			% >> P = tomography.project(volume);   projects 'volume'.
+			
+			if ~this.initialized 
+				this.initialize();
+			end
+			
+			if nargin == 1 % tomography.project();
+				P = this.project_c(this.V);
+			
+			elseif nargin == 2 % tomography.project(volume); 
+				P = this.project_c(volume);
+			end
+			
+			% store
+			if numel(this.sinogram) == 0
+				this.sinogram = P;
+			end
+		end
+		
+		%------------------------------------------------------------------
+		function V = reconstruct(this, iterations)
+			% Compute reconstruction.
+			% Uses tomography.sinogram
+			% Initial solution (if available) should be stored in tomography.V  
+			% >> V = tomography.reconstruct(iterations);
+			
+			if ~this.initialized 
+				this.initialize();
+			end			
+			
+			this.V = this.reconstruct_c(this.sinogram, this.V, [], iterations);
+			if strcmp(this.inner_circle,'yes')
+				this.V = this.selectROI(this.V);
+			end
+			V = this.V;
+		end
+
+		%------------------------------------------------------------------
+		function I = reconstructMask(this, mask, iterations)
+			% Compute reconstruction with mask.
+			% Uses tomography.sinogram  
+			% Initial solution (if available) should be stored in tomography.V  
+			% >> V = tomography.reconstructMask(mask, iterations);			
+			
+			if ~this.initialized 
+				this.initialize();
+			end
+			
+			if strcmp(this.inner_circle,'yes')
+				mask = this.selectROI(mask);
+			end
+			I = this.reconstruct_c(this.sinogram, this.V, mask, iterations);
+		end	
+		%------------------------------------------------------------------
+	
+	end
+	
+	%----------------------------------------------------------------------
+	methods (Access = protected)
+	
+		%------------------------------------------------------------------
+		% Protected function: create FP
+		function sinogram = project_c(this, volume)
+			
+			if this.initialized == 0
+				error('IterativeTomography not initialized');
+			end
+            
+			% data is stored in astra memory
+			if numel(volume) == 1
+				
+				if strcmp(this.gpu, 'yes')
+					sinogram_tmp = astra_create_sino_cuda(volume, this.proj_geom_sr, this.vol_geom, this.gpu_core);
+				else
+					sinogram_tmp = astra_create_sino(volume, this.proj_id);
+				end
+				
+				% sinogram downsampling
+				if this.superresolution > 1
+					sinogram_data = astra_mex_data2d('get', sinogram_tmp);
+					astra_mex_data2d('delete', sinogram_tmp);
+					sinogram_data = downsample_sinogram(sinogram_data, this.superresolution);	
+					%if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
+					%	sinogram = sinogram / this.superresolution;
+					%end
+					sinogram = astra_mex_data2d('create','sino', this.proj_geom, sinogram_data);
+				else
+					sinogram = sinogram_tmp;
+				end
+				
+			% data is stored in matlab memory	
+			else
+			
+				% 2D and 3D slice by slice
+				sinogram_tmp = zeros([astra_geom_size(this.proj_geom_sr), size(volume,3)]);
+				sinogram_tmp2 = zeros([astra_geom_size(this.proj_geom), size(volume,3)]);
+				for slice = 1:size(volume,3)
+					
+					if strcmp(this.gpu, 'yes')
+						[tmp_id, sinogram_tmp2(:,:,slice)] = astra_create_sino_sampling(volume(:,:,slice), this.proj_geom, this.vol_geom, this.gpu_core, this.superresolution);
+						astra_mex_data2d('delete', tmp_id);
+					else
+						[tmp_id, tmp] = astra_create_sino(volume(:,:,slice), this.proj_id_sr);
+						sinogram_tmp2(:,:,slice) = downsample_sinogram(tmp, this.superresolution) * (this.superresolution^2);
+						astra_mex_data2d('delete', tmp_id);
+					end
+					
+				end
+				
+				% sinogram downsampling
+				if strcmp(this.gpu, 'yes')
+					%sinogram = downsample_sinogram(sinogram_tmp, this.superresolution);
+					sinogram2 = sinogram_tmp2;
+					if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
+						sinogram2 = sinogram2 / this.superresolution;
+					elseif strcmp(this.proj_geom.type,'parallel')
+						sinogram2 = sinogram2 / (this.superresolution * this.superresolution);
+					end
+					sinogram = sinogram2;
+				else
+					sinogram = sinogram_tmp2;
+				end
+				
+			end
+			
+		end		
+
+		%------------------------------------------------------------------
+		% Protected function: reconstruct
+		function V = reconstruct_c(this, sinogram, V0, mask, iterations)
+			
+			if this.initialized == 0
+				error('IterativeTomography not initialized');
+			end			
+			
+			% data is stored in astra memory
+			if numel(sinogram) == 1
+				V = this.reconstruct_c_astra(sinogram, V0, mask, iterations);
+			
+			% data is stored in matlab memory
+			else
+				V = this.reconstruct_c_matlab(sinogram, V0, mask, iterations);
+			end
+			
+		end
+		
+		%------------------------------------------------------------------
+		% Protected function: reconstruct (data in matlab)	
+		function V = reconstruct_c_matlab(this, sinogram, V0, mask, iterations)
+			
+			if this.initialized == 0
+				error('IterativeTomography not initialized');
+			end						
+			
+			% parse method
+			method2 = upper(this.method);
+			if strcmp(method2, 'SART') || strcmp(method2, 'SART_CUDA')
+				iterations = iterations * size(sinogram,1);
+			elseif strcmp(method2, 'ART')
+				iterations = iterations * numel(sinogram);
+			end
+			
+			% create data objects
+			V = zeros(this.vol_geom.GridRowCount, this.vol_geom.GridColCount, size(sinogram,3));
+			reconstruction_id = astra_mex_data2d('create', '-vol', this.vol_geom);
+			sinogram_id = astra_mex_data2d('create', '-sino', this.proj_geom);
+			if numel(mask) > 0
+				mask_id = astra_mex_data2d('create', '-vol', this.vol_geom);
+			end
+			
+			% algorithm configuration
+			cfg = this.cfg_base;
+			cfg.ProjectionDataId = sinogram_id;
+			cfg.ReconstructionDataId = reconstruction_id;
+			if numel(mask) > 0
+				cfg.option.ReconstructionMaskId = mask_id;
+			end
+			alg_id = astra_mex_algorithm('create', cfg);
+			
+			% loop slices
+			for slice = 1:size(sinogram,3)
+				
+				% fetch slice of initial reconstruction
+				if numel(V0) > 0
+					astra_mex_data2d('store', reconstruction_id, V0(:,:,slice));
+				else
+					astra_mex_data2d('store', reconstruction_id, 0);
+				end
+				
+				% fetch slice of sinogram
+				astra_mex_data2d('store', sinogram_id, sinogram(:,:,slice));
+				
+				% fecth slice of mask
+				if numel(mask) > 0
+					astra_mex_data2d('store', mask_id, mask(:,:,slice));    
+				end
+				
+				% iterate
+				astra_mex_algorithm('iterate', alg_id, iterations);
+				
+				% fetch data
+				V(:,:,slice) = astra_mex_data2d('get', reconstruction_id);
+				
+			end				
+			
+			% correct attenuation factors for super-resolution
+			if this.superresolution > 1 && strcmp(this.gpu,'yes')
+				if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
+					if numel(mask) > 0
+						V(mask > 0) = V(mask > 0) ./ this.superresolution;
+					else
+						V = V ./ this.superresolution;
+					end
+				end		
+			end
+			
+			% garbage collection
+			astra_mex_algorithm('delete', alg_id);
+			astra_mex_data2d('delete', sinogram_id, reconstruction_id);
+			if numel(mask) > 0
+				astra_mex_data2d('delete', mask_id);
+			end	
+			
+		end
+		
+		%------------------------------------------------------------------
+		% Protected function: reconstruct (data in astra)	
+		function V = reconstruct_c_astra(this, sinogram, V0, mask, iterations)
+			
+			if this.initialized == 0
+				error('IterativeTomography not initialized');
+			end			
+		
+			if numel(V0) > 1 || numel(mask) > 1 || numel(sinogram) > 1
+				error('Not all required data is stored in the astra memory');
+			end			
+			
+			if numel(V0) == 0
+				V0 = astra_mex_data2d('create', '-vol', this.vol_geom, 0);
+			end
+			
+			% parse method
+			method2 = upper(this.method);
+			if strcmp(method2, 'SART') || strcmp(method2, 'SART_CUDA')
+				iterations = iterations * astra_geom_size(this.proj_geom, 1);
+			elseif strcmp(method2, 'ART')
+				s = astra_geom_size(this.proj_geom);
+				iterations = iterations * s(1) * s(2);
+			end	
+			
+			% algorithm configuration
+			cfg = this.cfg_base;
+			cfg.ProjectionDataId = sinogram;
+			cfg.ReconstructionDataId = V0;
+			if numel(mask) > 0
+				cfg.option.ReconstructionMaskId = mask;
+			end
+			alg_id = astra_mex_algorithm('create', cfg);				
+			
+			% iterate
+			astra_mex_algorithm('iterate', alg_id, iterations);	
+			
+			% fetch data
+			V = V0;
+			
+			% correct attenuation factors for super-resolution
+			if this.superresolution > 1
+				if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
+					if numel(mask) > 0
+						astra_data_op_masked('$1./s1', [V V], [this.superresolution this.superresolution], mask, this.gpu_core);
+					else
+						astra_data_op('$1./s1', [V V], [this.superresolution this.superresolution], this.gpu_core);
+					end
+				end	
+			end
+			
+			% garbage collection
+			astra_mex_algorithm('delete', alg_id);
+			
+		end
+		
+		%------------------------------------------------------------------
+		function V_out = selectROI(~, V_in)
+			
+			if numel(V_in) == 1
+				cfg = astra_struct('RoiSelect_CUDA');
+				cfg.DataId = V_in;
+				alg_id = astra_mex_algorithm('create',cfg);
+				astra_mex_algorithm('run', alg_id);
+				astra_mex_algorithm('delete', alg_id);
+				V_out = V_in;
+			else
+				V_out = ROIselectfull(V_in, min([size(V_in,1), size(V_in,2)]));
+			end
+			
+		end
+		%------------------------------------------------------------------
+		
+	end
+		
+end
+
diff --git a/matlab/algorithms/DART/IterativeTomography3D.m b/matlab/algorithms/DART/IterativeTomography3D.m
new file mode 100644
index 0000000..3f89457
--- /dev/null
+++ b/matlab/algorithms/DART/IterativeTomography3D.m
@@ -0,0 +1,433 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef IterativeTomography3D < matlab.mixin.Copyable
+
+	% Algorithm class for 3D Iterative Tomography.
+
+	%----------------------------------------------------------------------
+	properties (SetAccess=public, GetAccess=public)				
+		superresolution		= 1;				% SETTING: Volume upsampling factor.
+		proj_type			= 'linear';			% SETTING: Projector type, only when gpu='no'.
+		method				= 'SIRT3D_CUDA';		% SETTING: Iterative method (see ASTRA toolbox documentation).
+		gpu					= 'yes';			% SETTING: Use gpu? {'yes', 'no'}
+		gpu_core			= 0;				% SETTING: Which gpu core to use? Only when gpu='yes'.
+		inner_circle		= 'yes';			% SETTING: Do roi only? {'yes', 'no'}
+		image_size			= [];				% SETTING: Overwrite default reconstruction size.  Only if no vol_geom is specified.
+		use_minc			= 'no';				% SETTING: Use minimum constraint. {'no', 'yes'}
+        maxc                = +Inf;             % SETTING: Maximum constraint. +Inf means off.
+	end
+	%----------------------------------------------------------------------
+	properties (SetAccess=public, GetAccess=public)	
+		sinogram			= [];				% DATA: Projection data.
+		proj_geom			= [];				% DATA: Projection geometry.
+		V					= [];				% DATA: Volume data.  Also used to set initial estimate (optional).
+		vol_geom			= [];				% DATA: Volume geometry.
+	end
+	%----------------------------------------------------------------------
+	properties (SetAccess=private, GetAccess=public)
+		initialized     	= 0;				% Is this object initialized?
+	end
+	%----------------------------------------------------------------------
+	properties (SetAccess=protected, GetAccess=protected)
+		proj_geom_sr		= [];				% PROTECTED: geometry of sinogram (with super-resolution)
+        proj_id             = [];				% PROTECTED: astra id of projector (when gpu='no')
+        proj_id_sr          = [];				% PROTECTED: astra id of super-resolution projector (when gpu='no')	
+		cfg_base			= struct();			% PROTECTED: base configuration structure for the reconstruction algorithm.
+    end
+	%----------------------------------------------------------------------
+
+	methods (Access=public)
+		
+		%------------------------------------------------------------------
+		function this = IterativeTomography(varargin)
+			% Constructor
+			% >> tomography = IterativeTomography(); 
+			% >> tomography = IterativeTomography(base); base struct should contain fields 'sinogram' and 'proj_geom'.
+			% >> tomography = IterativeTomography(base_filename); mat-file should contain 'sinogram' and 'proj_geom'.
+			% >> tomography = IterativeTomography(proj_geom, vol_geom);
+			% >> tomography = IterativeTomography(sinogram, proj_geom);
+			
+			% Input: IterativeTomography(base)
+			if nargin == 1
+				if ischar(varargin{1})
+					this.sinogram = load(varargin{1},'sinogram');
+					this.proj_geom = load(varargin{1},'proj_geom');				
+				else
+					this.sinogram = varargin{1}.sinogram;
+					this.proj_geom = varargin{1}.proj_geom;
+				end
+			end
+			% Input: IterativeTomography(proj_geom, vol_geom)
+			if nargin == 2 && isstruct(varargin{1}) && isstruct(varargin{2})
+				this.proj_geom = varargin{1};
+				this.vol_geom =	varargin{2};
+			% Input: IterativeTomography(sinogram, proj_geom)
+			elseif nargin == 2
+				this.sinogram = varargin{1};
+				this.proj_geom = varargin{2};
+			end		
+
+		end
+
+		%------------------------------------------------------------------
+		function delete(this)
+			% Destructor
+			% >> clear tomography;
+			if strcmp(this.gpu,'no') && numel(this.proj_id) > 0
+				astra_mex_projector('delete', this.proj_id, this.proj_id_sr);
+			end
+		end
+			
+		%------------------------------------------------------------------
+		function settings = getsettings(this)
+			% Returns a structure containing all settings of this object.
+			% >> settings = tomography.getsettings();
+			settings.superresolution	= this.superresolution;
+			settings.proj_type			= this.proj_type;
+			settings.method				= this.method;
+			settings.gpu				= this.gpu;
+			settings.gpu_core			= this.gpu_core;
+			settings.inner_circle		= this.inner_circle;
+			settings.image_size			= this.image_size;	
+			settings.use_minc			= this.use_minc;
+            settings.maxc               = this.maxc;
+		end		
+		
+		%------------------------------------------------------------------
+		function ok = initialize(this)
+			% Initialize this object.  Returns 1 if succesful.
+			% >> tomography.initialize();
+			
+% 			% create projection geometry with super-resolution
+% 			this.proj_geom_sr = astra_geom_superresolution(this.proj_geom, this.superresolution);		
+			
+			% if no volume geometry is specified by the user: create volume geometry
+			if numel(this.vol_geom) == 0
+				if numel(this.image_size) < 2
+					this.image_size(1) = this.proj_geom.DetectorRowCount;
+					this.image_size(2) = this.proj_geom.DetectorColCount;
+				end
+				this.vol_geom = astra_create_vol_geom(this.proj_geom.DetectorColCount, this.proj_geom.DetectorColCount, this.proj_geom.DetectorRowCount);
+			else
+				this.image_size(1) = this.vol_geom.GridRowCount;
+				this.image_size(2) = this.vol_geom.GridColCount;
+			end
+			
+			% create projector
+			if strcmp(this.gpu, 'no')
+				this.proj_id = astra_create_projector(this.proj_type, this.proj_geom, this.vol_geom);
+				this.proj_id_sr = astra_create_projector(this.proj_type, this.proj_geom_sr, this.vol_geom);			
+			end
+			
+			% create reconstruction configuration
+			this.cfg_base = astra_struct(upper(this.method));
+			if strcmp(this.gpu,'no')
+				this.cfg_base.ProjectorId = this.proj_id;
+				this.cfg_base.ProjectionGeometry = this.proj_geom;
+				this.cfg_base.ReconstructionGeometry = this.vol_geom;
+				this.cfg_base.option.ProjectionOrder = 'random';
+			end
+			this.cfg_base.option.DetectorSuperSampling = this.superresolution;
+			%this.cfg_base.option.DetectorSuperSampling = 8;
+			if strcmp(this.gpu,'yes')
+				this.cfg_base.option.GPUindex = this.gpu_core;
+			end
+			this.cfg_base.option.UseMinConstraint = this.use_minc;
+            if ~isinf(this.maxc)
+              this.cfg_base.option.UseMaxConstraint = 'yes';
+              this.cfg_base.option.MaxConstraintValue = this.maxc;
+            end
+		
+			this.initialized = 1;
+			ok = this.initialized;
+		end
+		
+		%------------------------------------------------------------------
+		function P = project(this, volume)
+			% Compute forward projection.  
+			% Stores sinogram in tomography.sinogram if it is still empty.
+			% >> P = tomography.project();   projects 'tomography.V'.
+			% >> P = tomography.project(volume);   projects 'volume'.
+			
+			if ~this.initialized 
+				this.initialize();
+			end
+			
+			if nargin == 1 % tomography.project();
+				P = this.project_c(this.V);
+			
+			elseif nargin == 2 % tomography.project(volume); 
+				P = this.project_c(volume);
+			end
+			
+			% store
+			if numel(this.sinogram) == 0
+				this.sinogram = P;
+			end
+		end
+		
+		%------------------------------------------------------------------
+		function V = reconstruct(this, iterations)
+			% Compute reconstruction.
+			% Uses tomography.sinogram
+			% Initial solution (if available) should be stored in tomography.V  
+			% >> V = tomography.reconstruct(iterations);
+			
+			if ~this.initialized 
+				this.initialize();
+			end			
+			
+			this.V = this.reconstruct_c(this.sinogram, this.V, [], iterations);
+			if strcmp(this.inner_circle,'yes')
+				this.V = this.selectROI(this.V);
+			end
+			V = this.V;
+		end
+
+		%------------------------------------------------------------------
+		function I = reconstructMask(this, mask, iterations)
+			% Compute reconstruction with mask.
+			% Uses tomography.sinogram  
+			% Initial solution (if available) should be stored in tomography.V  
+			% >> V = tomography.reconstructMask(mask, iterations);			
+			
+			if ~this.initialized 
+				this.initialize();
+			end
+			
+			if strcmp(this.inner_circle,'yes')
+				mask = this.selectROI(mask);
+			end
+			I = this.reconstruct_c(this.sinogram, this.V, mask, iterations);
+		end	
+		%------------------------------------------------------------------
+	
+	end
+	
+	%----------------------------------------------------------------------
+	methods (Access = protected)
+	
+		%------------------------------------------------------------------
+		% Protected function: create FP
+		function sinogram = project_c(this, volume)
+			
+			if this.initialized == 0
+				error('IterativeTomography not initialized');
+			end
+            
+			% data is stored in astra memory
+			if numel(volume) == 1
+				
+				if strcmp(this.gpu, 'yes')
+					sinogram_tmp = astra_create_sino_cuda(volume, this.proj_geom_sr, this.vol_geom, this.gpu_core);
+				else
+					sinogram_tmp = astra_create_sino(volume, this.proj_id);
+				end
+				
+				% sinogram downsampling
+				if this.superresolution > 1
+					sinogram_data = astra_mex_data2d('get', sinogram_tmp);
+					astra_mex_data2d('delete', sinogram_tmp);
+					sinogram_data = downsample_sinogram(sinogram_data, this.superresolution);	
+					%if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
+					%	sinogram = sinogram / this.superresolution;
+					%end
+					sinogram = astra_mex_data2d('create','sino', this.proj_geom, sinogram_data);
+				else
+					sinogram = sinogram_tmp;
+				end
+				
+			% data is stored in matlab memory	
+			else
+			
+				[tmp_id, sinogram] = astra_create_sino3d_cuda(volume, this.proj_geom, this.vol_geom);
+				astra_mex_data3d('delete', tmp_id);
+				
+			end
+			
+		end		
+
+		%------------------------------------------------------------------
+		% Protected function: reconstruct
+		function V = reconstruct_c(this, sinogram, V0, mask, iterations)
+			
+			if this.initialized == 0
+				error('IterativeTomography not initialized');
+			end			
+			
+			% data is stored in astra memory
+			if numel(sinogram) == 1
+				V = this.reconstruct_c_astra(sinogram, V0, mask, iterations);
+			
+			% data is stored in matlab memory
+			else
+				V = this.reconstruct_c_matlab(sinogram, V0, mask, iterations);
+			end
+			
+		end
+		
+		%------------------------------------------------------------------
+		% Protected function: reconstruct (data in matlab)	
+		function V = reconstruct_c_matlab(this, sinogram, V0, mask, iterations)
+			
+			if this.initialized == 0
+				error('IterativeTomography not initialized');
+			end						
+			
+			% parse method
+			method2 = upper(this.method);
+			if strcmp(method2, 'SART') || strcmp(method2, 'SART_CUDA')
+				iterations = iterations * size(sinogram,1);
+			elseif strcmp(method2, 'ART')
+				iterations = iterations * numel(sinogram);
+			end
+			
+			% create data objects
+% 			V = zeros(this.vol_geom.GridRowCount, this.vol_geom.GridColCount, size(sinogram,3));
+			reconstruction_id = astra_mex_data3d('create', '-vol', this.vol_geom);
+			sinogram_id = astra_mex_data3d('create', '-proj3d', this.proj_geom);
+			if numel(mask) > 0
+				mask_id = astra_mex_data3d('create', '-vol', this.vol_geom);
+			end
+			
+			% algorithm configuration
+			cfg = this.cfg_base;
+			cfg.ProjectionDataId = sinogram_id;
+			cfg.ReconstructionDataId = reconstruction_id;
+			if numel(mask) > 0
+				cfg.option.ReconstructionMaskId = mask_id;
+			end
+			alg_id = astra_mex_algorithm('create', cfg);
+			
+% 			% loop slices
+% 			for slice = 1:size(sinogram,3)
+				
+				% fetch slice of initial reconstruction
+				if numel(V0) > 0
+					astra_mex_data3d('store', reconstruction_id, V0);
+				else
+					astra_mex_data3d('store', reconstruction_id, 0);
+				end
+				
+				% fetch slice of sinogram
+				astra_mex_data3d('store', sinogram_id, sinogram);
+				
+				% fecth slice of mask
+				if numel(mask) > 0
+					astra_mex_data3d('store', mask_id, mask);    
+				end
+				
+				% iterate
+				astra_mex_algorithm('iterate', alg_id, iterations);
+				
+				% fetch data
+				V = astra_mex_data3d('get', reconstruction_id);
+				
+%			end
+			
+			% correct attenuation factors for super-resolution
+			if this.superresolution > 1 && strcmp(this.gpu,'yes')
+				if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
+					if numel(mask) > 0
+						V(mask > 0) = V(mask > 0) ./ this.superresolution;
+					else
+						V = V ./ this.superresolution;
+					end
+				end		
+			end
+			
+			% garbage collection
+			astra_mex_algorithm('delete', alg_id);
+			astra_mex_data3d('delete', sinogram_id, reconstruction_id);
+			if numel(mask) > 0
+				astra_mex_data3d('delete', mask_id);
+			end	
+			
+		end
+		
+		%------------------------------------------------------------------
+		% Protected function: reconstruct (data in astra)	
+		function V = reconstruct_c_astra(this, sinogram, V0, mask, iterations)
+			
+			if this.initialized == 0
+				error('IterativeTomography not initialized');
+			end			
+		
+			if numel(V0) > 1 || numel(mask) > 1 || numel(sinogram) > 1
+				error('Not all required data is stored in the astra memory');
+			end			
+			
+			if numel(V0) == 0
+				V0 = astra_mex_data2d('create', '-vol', this.vol_geom, 0);
+			end
+			
+			% parse method
+			method2 = upper(this.method);
+			if strcmp(method2, 'SART') || strcmp(method2, 'SART_CUDA')
+				iterations = iterations * astra_geom_size(this.proj_geom, 1);
+			elseif strcmp(method2, 'ART')
+				s = astra_geom_size(this.proj_geom);
+				iterations = iterations * s(1) * s(2);
+			end	
+			
+			% algorithm configuration
+			cfg = this.cfg_base;
+			cfg.ProjectionDataId = sinogram;
+			cfg.ReconstructionDataId = V0;
+			if numel(mask) > 0
+				cfg.option.ReconstructionMaskId = mask;
+			end
+			alg_id = astra_mex_algorithm('create', cfg);				
+			
+			% iterate
+			astra_mex_algorithm('iterate', alg_id, iterations);	
+			
+			% fetch data
+			V = V0;
+			
+			% correct attenuation factors for super-resolution
+			if this.superresolution > 1
+				if strcmp(this.proj_geom.type,'fanflat_vec') || strcmp(this.proj_geom.type,'fanflat')
+					if numel(mask) > 0
+						astra_data_op_masked('$1./s1', [V V], [this.superresolution this.superresolution], mask, this.gpu_core);
+					else
+						astra_data_op('$1./s1', [V V], [this.superresolution this.superresolution], this.gpu_core);
+					end
+				end	
+			end
+			
+			% garbage collection
+			astra_mex_algorithm('delete', alg_id);
+			
+		end
+		
+		%------------------------------------------------------------------
+		function V_out = selectROI(~, V_in)
+			
+			if numel(V_in) == 1
+				cfg = astra_struct('RoiSelect_CUDA');
+				cfg.DataId = V_in;
+				alg_id = astra_mex_algorithm('create',cfg);
+				astra_mex_algorithm('run', alg_id);
+				astra_mex_algorithm('delete', alg_id);
+				V_out = V_in;
+			else
+				V_out = ROIselectfull(V_in, min([size(V_in,1), size(V_in,2)]));
+			end
+			
+		end
+		%------------------------------------------------------------------
+		
+	end
+		
+end
+
diff --git a/matlab/algorithms/DART/Kernels.m b/matlab/algorithms/DART/Kernels.m
new file mode 100644
index 0000000..b5e3134
--- /dev/null
+++ b/matlab/algorithms/DART/Kernels.m
@@ -0,0 +1,68 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef Kernels
+	%KERNELS Summary of this class goes here
+	%   Detailed explanation goes here
+	
+	properties
+		
+	end
+	
+	methods(Static)
+
+		function K = BinaryPixelKernel(radius, conn)
+			
+			if nargin < 2
+				conn = 8;
+			end
+			
+			% 2D, 4conn
+			if conn == 4
+				K = [0 1 0; 1 1 1; 0 1 0];
+				for i = 2:radius
+					K = conv2(K,K);
+				end
+				K = double(K >= 1);
+			
+			% 2D, 8conn
+			elseif conn == 8
+				K = ones(2*radius+1, 2*radius+1);
+				
+			% 3D, 6conn	
+			elseif conn == 6
+				K = zeros(3,3,3);
+				K(:,:,1) = [0 0 0; 0 1 0; 0 0 0];
+				K(:,:,2) = [0 1 0; 1 1 1; 0 1 0];
+				K(:,:,3) = [0 0 0; 0 1 0; 0 0 0];
+				for i = 2:radius
+					K = convn(K,K);
+				end
+				K = double(K >= 1);				
+				
+			% 2D, 27conn
+			elseif conn == 26
+				K = ones(2*radius+1, 2*radius+1, 2*radius+1);
+				
+			else
+				disp('Invalid conn')
+			end
+		end
+
+		
+
+	
+		
+		
+	end
+	
+end
+
diff --git a/matlab/algorithms/DART/MaskingDefault.m b/matlab/algorithms/DART/MaskingDefault.m
new file mode 100644
index 0000000..6bd25a5
--- /dev/null
+++ b/matlab/algorithms/DART/MaskingDefault.m
@@ -0,0 +1,212 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef MaskingDefault < matlab.mixin.Copyable
+
+	% Default policy class for masking for DART.
+	
+	%----------------------------------------------------------------------
+	properties (Access=public)
+		
+		radius			= 1;			% SETTING: Radius of masking kernel.
+		conn			= 8;			% SETTING: Connectivity window. For 2D: 4 or 8.  For 3D: 6 or 26.
+		edge_threshold	= 1;			% SETTING: Number of pixels in the window that should be different.
+		random			= 0.1;			% SETTING: Percentage of random points.  Between 0 and 1.
+		gpu				= 'yes';		% SETTING: Use gpu? {'yes', 'no'}
+		gpu_core		= 0;			% SETTING: Which gpu core to use, only when gpu='yes'.
+				
+	end
+	
+	%----------------------------------------------------------------------
+	methods (Access=public)
+		
+		%------------------------------------------------------------------
+		function settings = getsettings(this)
+			% Returns a structure containing all settings of this object.
+			% >> settings = DART.masking.getsettings();				
+			settings.radius				= this.radius;
+			settings.conn				= this.conn;
+			settings.edge_threshold		= this.edge_threshold;
+			settings.random				= this.random;
+		end
+		
+		%------------------------------------------------------------------
+		function Mask = apply(this, ~, S)
+			% Applies masking.
+			% >> Mask = DART.segmentation.apply(DART, S);	
+			
+			% 2D, one slice
+			if size(S,3) == 1
+				if strcmp(this.gpu,'yes')
+					Mask = this.apply_2D_gpu(S);
+				else
+					Mask = this.apply_2D(S);
+				end
+						
+			% 3D, slice by slice
+			elseif this.conn == 4 || this.conn == 8
+				Mask = zeros(size(S));
+				for slice = 1:size(S,3)
+					if strcmp(this.gpu,'yes')
+						Mask(:,:,slice) = this.apply_2D_gpu(S(:,:,slice));
+					else
+						Mask(:,:,slice) = this.apply_2D(S(:,:,slice));						
+					end
+				end
+			
+			% 3D, full
+			else
+				if strcmp(this.gpu,'yes')
+					Mask = this.apply_3D_gpu(S);	
+				else
+					Mask = this.apply_3D(S);	
+				end
+			end
+			
+		end
+		
+	end
+		
+	%----------------------------------------------------------------------
+	methods (Access=protected)
+		
+		%------------------------------------------------------------------
+		function Mask = apply_2D_gpu(this, S)
+		
+			vol_geom = astra_create_vol_geom(size(S));
+			data_id = astra_mex_data2d('create', '-vol', vol_geom, S);
+			mask_id = astra_mex_data2d('create', '-vol', vol_geom, 0);
+
+			cfg = astra_struct('DARTMASK_CUDA');
+			cfg.SegmentationDataId = data_id;
+			cfg.MaskDataId = mask_id;
+			cfg.option.GPUindex = this.gpu_core;
+			cfg.option.Connectivity = this.conn;
+			cfg.option.Threshold = this.edge_threshold;
+			cfg.option.Radius = this.radius;
+			
+			alg_id = astra_mex_algorithm('create',cfg);	
+			astra_mex_algorithm('iterate',alg_id,1);
+			Mask = astra_mex_data2d('get', mask_id);
+		
+			astra_mex_algorithm('delete', alg_id);
+			astra_mex_data2d('delete', data_id, mask_id);
+			
+			% random
+			RandomField = double(rand(size(S)) < this.random);
+
+			% combine
+			Mask = or(Mask, RandomField);			
+			
+		end			
+		
+		%------------------------------------------------------------------
+		function Mask = apply_2D(this, S)
+		
+			r = this.radius;
+			w = 2 * r + 1;
+			
+			kernel = Kernels.BinaryPixelKernel(r, this.conn);
+			
+			% edges
+			Xlarge = zeros(size(S,1)+w-1, size(S,2)+w-1); 
+			Xlarge(1+r:end-r, 1+r:end-r) = S;
+
+			Edges = zeros(size(S));
+			for s = -r:r
+				for t = -r:r
+					if kernel(s+r+1, t+r+1) == 0
+						continue
+					end
+					Temp = abs(Xlarge(1+r:end-r, 1+r:end-r) - Xlarge(1+r+s:end-r+s, 1+r+t:end-r+t));
+					Edges(Temp > eps) = Edges(Temp > eps) + 1;
+				end
+			end
+
+			Edges = Edges > this.edge_threshold;
+			
+			% random
+			RandomField = double(rand(size(S)) < this.random);
+
+			% combine
+			Mask = or(Edges, RandomField);
+			
+		end		
+		
+		%------------------------------------------------------------------
+		function Mask = apply_3D(this, S)
+		
+			r = this.radius;
+			w = 2 * r + 1;
+			
+			kernel = Kernels.BinaryPixelKernel(r, this.conn);			
+			
+			% edges
+			Xlarge = zeros(size(S,1)+w-1, size(S,2)+w-1, size(S,3)+w-1); 
+			Xlarge(1+r:end-r, 1+r:end-r, 1+r:end-r) = S;
+
+			Edges = zeros(size(S));
+			for s = -r:r
+				for t = -r:r
+					for u = -r:r
+						if kernel(s+r+1, t+r+1, u+r+1) == 0
+							continue
+						end
+						Temp = abs(Xlarge(1+r:end-r, 1+r:end-r, 1+r:end-r) - Xlarge(1+r+s:end-r+s, 1+r+t:end-r+t, 1+r+u:end-r+u));
+						Edges(Temp > eps) = 1;
+					end
+				end
+			end
+			
+			clear Xlarge;
+
+			% random
+			RandomField = double(rand(size(S)) < this.random);
+
+			% combine
+			Mask = or(Edges, RandomField);
+			
+		end		
+		
+		%------------------------------------------------------------------
+		function Mask = apply_3D_gpu(this, S)
+	
+			vol_geom = astra_create_vol_geom(size(S));
+			data_id = astra_mex_data3d('create', '-vol', vol_geom, S);
+
+			cfg = astra_struct('DARTMASK3D_CUDA');
+			cfg.SegmentationDataId = data_id;
+			cfg.MaskDataId = data_id;
+			cfg.option.GPUindex = this.gpu_core;
+			cfg.option.Connectivity = this.conn;
+			cfg.option.Threshold = this.edge_threshold;
+			cfg.option.Radius = this.radius;
+			
+			alg_id = astra_mex_algorithm('create',cfg);	
+			astra_mex_algorithm('iterate',alg_id,1);
+			Mask = astra_mex_data3d('get', data_id);
+		
+			astra_mex_algorithm('delete', alg_id);
+			astra_mex_data3d('delete', data_id);
+			
+			% random
+			RandomField = double(rand(size(S)) < this.random);
+
+			% combine
+			Mask = or(Mask, RandomField);			
+			
+		end					
+		
+		%------------------------------------------------------------------
+	end
+	
+end
+
diff --git a/matlab/algorithms/DART/MaskingGPU.m b/matlab/algorithms/DART/MaskingGPU.m
new file mode 100644
index 0000000..c4ef2b7
--- /dev/null
+++ b/matlab/algorithms/DART/MaskingGPU.m
@@ -0,0 +1,94 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef MaskingGPU < matlab.mixin.Copyable
+
+	% Policy class for masking for DART with GPU accelerated code (deprecated).
+	
+	%----------------------------------------------------------------------
+	properties (Access=public)
+		
+		radius			= 1;			% SETTING: Radius of masking kernel.
+		conn			= 8;			% SETTING: Connectivity window. For 2D: 4 or 8.  For 3D: 6 or 26.
+		edge_threshold	= 1;			% SETTING: Number of pixels in the window that should be different.
+		gpu_core		= 0;			% SETTING:
+		random			= 0.1;			% SETTING: Percentage of random points.  Between 0 and 1.
+		
+	end
+	
+	%----------------------------------------------------------------------
+	methods (Access=public)
+		
+		%------------------------------------------------------------------
+		function settings = getsettings(this)
+			% Returns a structure containing all settings of this object.
+			% >> settings = DART.masking.getsettings();				
+			settings.radius				= this.radius;
+			settings.conn				= this.conn;
+			settings.edge_threshold		= this.edge_threshold;
+			settings.random				= this.random;
+		end
+		
+		%------------------------------------------------------------------
+		function Mask = apply(this, ~, V_in)
+			% Applies masking.
+			% >> Mask = DART.segmentation.apply(DART, V_in);	
+			
+			% 2D, one slice
+			if size(V_in,3) == 1
+				Mask = this.apply_2D(V_in);
+						
+			% 3D, slice by slice
+			elseif this.conn == 4 || this.conn == 8
+				Mask = zeros(size(V_in));
+				for slice = 1:size(V_in,3)
+					Mask(:,:,slice) = this.apply_2D(V_in(:,:,slice)); 
+				end
+			
+			% 3D, full
+			else
+				error('Full 3D masking on GPU not implemented.')
+			end
+			
+		end
+		
+	end
+		
+	%----------------------------------------------------------------------
+	methods (Access=protected)
+		
+		%------------------------------------------------------------------
+		function Mask = apply_2D(this, S)
+		
+			vol_geom = astra_create_vol_geom(size(S));
+			data_id = astra_mex_data2d('create', '-vol', vol_geom, S);
+			mask_id = astra_mex_data2d('create', '-vol', vol_geom, 0);
+
+			cfg = astra_struct('DARTMASK_CUDA');
+			cfg.SegmentationDataId = data_id;
+			cfg.MaskDataId = mask_id;
+			cfg.option.GPUindex = this.gpu_core;
+			%cfg.option.Connectivity = this.conn;
+			
+			alg_id = astra_mex_algorithm('create',cfg);	
+			astra_mex_algorithm('iterate',alg_id,1);
+			Mask = astra_mex_data2d('get', mask_id);
+		
+			astra_mex_algorithm('delete', alg_id);
+			astra_mex_data2d('delete', data_id, mask_id);
+			
+		end	
+	end
+
+
+	
+end
+
diff --git a/matlab/algorithms/DART/OutputDefault.m b/matlab/algorithms/DART/OutputDefault.m
new file mode 100644
index 0000000..a34a430
--- /dev/null
+++ b/matlab/algorithms/DART/OutputDefault.m
@@ -0,0 +1,173 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef OutputDefault < matlab.mixin.Copyable
+	 
+	% Default policy class for output for DART.
+	
+	properties (Access=public)
+		
+		directory				= '';		% SETTING: Directory to save output.
+		pre						= '';		% SETTING: Prefix of output.
+
+		save_images				= 'no';		% SETTING: Save the images. 'no', 'yes' (='S') OR {'S', 'I', 'Mask', 'P'}
+		save_results			= 'no';		% SETTING: Save the results. 'yes', 'no' OR {'base', 'stats', 'settings', 'S', 'V', 'V0'}
+		save_object				= 'no';		% SETTING: Save the DART object. {'no','yes'}
+
+		save_interval			= 1;		% SETTING: # DART iteration between saves.
+		save_images_interval	= [];		% SETTING: Overwrite interval for save_images.
+		save_results_interval	= [];		% SETTING: Overwrite interval for save_results.
+		save_object_interval	= [];		% SETTING: Overwrite interval for save_object.
+
+		slices					= 1;		% SETTING: In case of 3D, which slices to save?		
+		
+		verbose					= 'no';		% SETTING: Verbose? {'no','yes'}
+		
+	end
+	
+	methods (Access=public)
+		
+		%------------------------------------------------------------------
+		function pre_initial_iteration(this, ~)
+			if strcmp(this.verbose,'yes')
+				tic
+				fprintf(1, 'initial iteration...');
+			end
+		end
+
+		%------------------------------------------------------------------
+		function post_initial_iteration(this, ~)
+			if strcmp(this.verbose,'yes')
+				t = toc;
+				fprintf(1, 'done in %f s.\n', t);
+			end
+		end
+		
+		%------------------------------------------------------------------
+		function pre_iteration(this, DART)
+			if strcmp(this.verbose,'yes')
+				tic;
+				fprintf(1, '%s dart iteration %d...', this.pre, DART.iterationcount);
+			end
+		end
+		
+		%------------------------------------------------------------------
+		function post_iteration(this, DART)
+
+			% print output
+			if strcmp(this.verbose,'yes')
+				t = toc;
+				s = DART.statistics.tostring(DART.stats);
+				fprintf(1, 'done in %0.2fs %s.\n', t, s);
+			end
+			
+			% save DART object	
+			if do_object(this, DART)
+				save(sprintf('%s%sobject_%i.mat', this.directory, this.pre, DART.iterationcount), '-v7.3', 'DART');					
+			end
+			
+			% save .mat
+			if do_results(this, DART)
+				base = DART.base;
+				stats = DART.stats;
+				S = DART.S;
+				V = DART.V;
+				V0 = DART.V0;
+				settings = DART.getsettings();
+				if ~iscell(this.save_results) 
+					save(sprintf('%s%sresults_%i.mat', this.directory, this.pre, DART.iterationcount), '-v7.3', 'base', 'stats', 'S', 'V', 'V0', 'settings');
+				else
+					string = [];
+					for i = 1:numel(this.save_results)
+						string = [string this.save_results{i} '|'];
+					end
+					save(sprintf('%s%sresults_%i.mat', this.directory, this.pre, DART.iterationcount), '-v7.3', '-regexp', string(1:end-1));
+				end
+			end	
+			
+			% save images
+			if do_images(this, DART)
+			
+				if ~iscell(this.save_images) && strcmp(this.save_images, 'yes')
+					output_image(this, DART, 'S')
+				elseif iscell(this.save_images)
+					for i = 1:numel(this.save_images)
+						output_image(this, DART, this.save_images{i});
+					end
+				end
+				
+			end					
+			
+		end
+		%------------------------------------------------------------------
+		
+	end
+	
+	%----------------------------------------------------------------------
+	methods (Access=private)
+		
+		function output_image(this, DART, data)
+			% 2D
+			if numel(size(DART.S)) == 2 
+				eval(['imwritesc(DART.' data ', sprintf(''%s%s' data '_%i.png'', this.directory, this.pre, DART.iterationcount))']);	
+			% 3D
+			elseif numel(size(DART.S)) == 3  
+				for slice = this.slices
+					eval(['imwritesc(DART.' data '(:,:,slice), sprintf(''%s%s' data '_%i_slice%i.png'', this.directory, this.pre, DART.iterationcount, slice))']);
+				end
+			end
+		end
+	
+		%------------------------------------------------------------------
+		function out = do_object(this, DART)
+			if strcmp(this.save_object,'no')
+				out = 0;
+				return
+			end			
+			if numel(this.save_object_interval) == 0 && mod(DART.iterationcount, this.save_interval) == 0 
+				out = 1;
+			elseif mod(DART.iterationcount, this.save_object_interval) == 0 
+				out = 1;
+			else
+				out = 0;
+			end
+		end				
+		%------------------------------------------------------------------
+		function out = do_results(this, DART)
+			if strcmp(this.save_results,'no')
+				out = 0;
+				return
+			end
+			if numel(this.save_results_interval) == 0 && mod(DART.iterationcount, this.save_interval) == 0 
+				out = 1;
+			elseif mod(DART.iterationcount, this.save_results_interval) == 0 
+				out = 1;
+			else
+				out = 0;
+			end
+		end		
+		
+		%------------------------------------------------------------------
+		function out = do_images(this, DART)
+			if numel(this.save_images_interval) == 0 && mod(DART.iterationcount, this.save_interval) == 0 
+				out = 1;
+			elseif mod(DART.iterationcount, this.save_images_interval) == 0 
+				out = 1;
+			else
+				out = 0;				
+			end
+		end
+		%------------------------------------------------------------------
+		
+	end
+	
+end
+
diff --git a/matlab/algorithms/DART/SegmentationDefault.m b/matlab/algorithms/DART/SegmentationDefault.m
new file mode 100644
index 0000000..c1d7d99
--- /dev/null
+++ b/matlab/algorithms/DART/SegmentationDefault.m
@@ -0,0 +1,55 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef SegmentationDefault < matlab.mixin.Copyable
+		
+	% Default policy class for segmentation for DART.
+	
+	%----------------------------------------------------------------------
+	properties (Access=public)
+		rho			= [];		% SETTING: Grey levels.
+		tau			= [];		% SETTING: Threshold values.
+	end
+	
+	%----------------------------------------------------------------------
+	methods (Access=public)
+		
+		%------------------------------------------------------------------
+		function settings = getsettings(this)
+			% Returns a structure containing all settings of this object.
+			% >> settings = DART.segmentation.getsettings();			
+			settings.rho	= this.rho;
+			settings.tau	= this.tau;
+		end
+	
+		%------------------------------------------------------------------
+		function this = estimate_grey_levels(this, ~, ~)
+			% Estimates grey levels
+			% >> DART.segmentation.estimate_grey_levels();				
+		end
+		
+		%------------------------------------------------------------------
+		function V_out = apply(this, ~, V_in)
+			% Applies segmentation.
+			% >> V_out = DART.segmentation.apply(DART, V_in);	
+			
+			V_out = ones(size(V_in)) * this.rho(1);
+			for n = 2:length(this.rho)
+				V_out(this.tau(n-1) < V_in) = this.rho(n);
+			end
+			
+		end
+		%------------------------------------------------------------------
+			
+	end
+	
+end
+
diff --git a/matlab/algorithms/DART/SmoothingDefault.m b/matlab/algorithms/DART/SmoothingDefault.m
new file mode 100644
index 0000000..58a8baa
--- /dev/null
+++ b/matlab/algorithms/DART/SmoothingDefault.m
@@ -0,0 +1,179 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef SmoothingDefault < matlab.mixin.Copyable
+	 
+	% Default policy class for smoothing for DART.
+	
+	%----------------------------------------------------------------------
+	properties (Access=public)
+		radius		= 1;			% SETTING: Radius of smoothing kernel.
+		b			= 0.1;			% SETTING: Intensity of smoothing.  Between 0 and 1.	
+		full3d		= 'yes';		% SETTING: smooth in 3D? {'yes','no'}
+		gpu			= 'yes';		% SETTING: Use gpu? {'yes', 'no'}
+		gpu_core	= 0;			% SETTING: Which gpu core to use, only when gpu='yes'.		
+	end
+	
+	
+	%----------------------------------------------------------------------
+	methods (Access=public)
+
+		%------------------------------------------------------------------
+		function settings = getsettings(this)
+			% Returns a structure containing all settings of this object.
+			% >> settings = DART.smoothing.getsettings();					
+			settings.radius				= this.radius;
+			settings.b					= this.b;
+			settings.full3d				= this.full3d;
+		end		
+		
+		%------------------------------------------------------------------
+		function V_out = apply(this, ~, V_in)
+			% Applies smoothing.
+			% >> V_out = DART.smoothing.apply(DART, V_in);	
+			
+			% 2D, one slice
+			if size(V_in,3) == 1
+				if strcmp(this.gpu,'yes')
+					V_out = this.apply_2D_gpu(V_in);
+				else
+					V_out = this.apply_2D(V_in);
+				end
+						
+			% 3D, slice by slice
+			elseif ~strcmp(this.full3d,'yes')
+				V_out = zeros(size(V_in));
+				for slice = 1:size(V_in,3)
+					if strcmp(this.gpu,'yes')
+						V_out(:,:,slice) = this.apply_2D_gpu(V_in(:,:,slice));
+					else
+						V_out(:,:,slice) = this.apply_2D(V_in(:,:,slice));						
+					end
+				end
+			
+			% 3D, full
+			else
+				if strcmp(this.gpu,'yes')
+					V_out = this.apply_3D_gpu(V_in);
+				else
+					V_out = this.apply_3D(V_in);
+				end
+			end
+			
+		end	
+			
+	end
+		
+	%----------------------------------------------------------------------
+	methods (Access=protected)
+		
+		%------------------------------------------------------------------
+		function V_out = apply_2D(this, V_in)
+			
+			r = this.radius;
+			w = 2 * r + 1;
+			
+			% Set Kernel
+			K = ones(w) * this.b / (w.^2-1); % edges
+			K(r+1,r+1) = 1 - this.b; % center
+			
+			% output window
+			V_out = zeros(size(V_in,1) + w-1, size(V_in,2) + w - 1);
+
+			% blur convolution
+			for s = -r:r 
+				for t = -r:r 
+					V_out(1+r+s:end-r+s, 1+r+t:end-r+t) = V_out(1+r+s:end-r+s, 1+r+t:end-r+t) + K(r+1+s, r+1+t) * V_in;
+				end
+			end
+			
+			% shrink output window
+			V_out = V_out(1+r:end-r, 1+r:end-r);
+			
+		end
+
+		%------------------------------------------------------------------
+		function V_out = apply_2D_gpu(this, V_in)
+			
+			vol_geom = astra_create_vol_geom(size(V_in));
+			in_id = astra_mex_data2d('create', '-vol', vol_geom, V_in);
+			out_id = astra_mex_data2d('create', '-vol', vol_geom, 0);
+
+			cfg = astra_struct('DARTSMOOTHING_CUDA');
+			cfg.InDataId = in_id;
+			cfg.OutDataId = out_id;
+			cfg.option.Intensity = this.b;
+			cfg.option.Radius = this.radius;			
+			cfg.option.GPUindex = this.gpu_core;
+			
+			alg_id = astra_mex_algorithm('create',cfg);	
+			astra_mex_algorithm('iterate',alg_id,1);
+			V_out = astra_mex_data2d('get', out_id);
+		
+			astra_mex_algorithm('delete', alg_id);
+			astra_mex_data2d('delete', in_id, out_id);
+			
+		end		
+		
+		%------------------------------------------------------------------
+		function I_out = apply_3D(this, I_in)
+			
+			r = this.radius;
+			w = 2 * r + 1;
+			
+			% Set Kernel
+			K = ones(w,w,w) * this.b / (w.^3-1); % edges
+			K(r+1,r+1,r+1) = 1 - this.b; % center
+			
+			% output window
+			I_out = zeros(size(I_in,1)+w-1, size(I_in,2)+w-1, size(I_in,3)+w-1);
+
+			% blur convolution
+			for s = -r:r 
+				for t = -r:r 
+					for u = -r:r 
+						I_out(1+r+s:end-r+s, 1+r+t:end-r+t, 1+r+u:end-r+u) = I_out(1+r+s:end-r+s, 1+r+t:end-r+t, 1+r+u:end-r+u) + K(r+1+s, r+1+t, r+1+u) * I_in;
+					end
+				end
+			end
+			
+			% shrink output window
+			I_out = I_out(1+r:end-r, 1+r:end-r, 1+r:end-r);
+			
+		end		
+		
+		%------------------------------------------------------------------
+		function V_out = apply_3D_gpu(this, V_in)
+
+			vol_geom = astra_create_vol_geom(size(V_in));
+			data_id = astra_mex_data3d('create', '-vol', vol_geom, V_in);
+
+			cfg = astra_struct('DARTSMOOTHING3D_CUDA');
+			cfg.InDataId = data_id;
+			cfg.OutDataId = data_id;
+			cfg.option.Intensity = this.b;
+			cfg.option.Radius = this.radius;
+			cfg.option.GPUindex = this.gpu_core;
+			
+			alg_id = astra_mex_algorithm('create',cfg);	
+			astra_mex_algorithm('iterate', alg_id, 1);
+			V_out = astra_mex_data3d('get', data_id);
+		
+			astra_mex_algorithm('delete', alg_id);
+			astra_mex_data3d('delete', data_id);
+			
+		end		
+		%------------------------------------------------------------------
+		
+	end
+	
+end
+
diff --git a/matlab/algorithms/DART/SmoothingGPU.m b/matlab/algorithms/DART/SmoothingGPU.m
new file mode 100644
index 0000000..857da37
--- /dev/null
+++ b/matlab/algorithms/DART/SmoothingGPU.m
@@ -0,0 +1,119 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef SmoothingGPU < matlab.mixin.Copyable
+	 
+	% Default policy class for smoothing for DART.
+	
+	%----------------------------------------------------------------------
+	properties (Access=public)
+		radius		= 1;			% SETTING: Radius of smoothing kernel.
+		b			= 0.1;			% SETTING: Intensity of smoothing.  Between 0 and 1.	
+		full3d		= 'yes';		% SETTING: smooth in 3D? {'yes','no'}
+		gpu_core	= 0;			% SETTING:		
+	end
+	
+	
+	%----------------------------------------------------------------------
+	methods (Access=public)
+
+		%------------------------------------------------------------------
+		function settings = getsettings(this)
+			% Returns a structure containing all settings of this object.
+			% >> settings = DART.smoothing.getsettings();					
+			settings.radius				= this.radius;
+			settings.b					= this.b;
+			settings.full3d				= this.full3d;
+		end		
+		
+		%------------------------------------------------------------------
+		function V_out = apply(this, ~, V_in)
+			% Applies smoothing.
+			% >> V_out = DART.smoothing.apply(DART, V_in);	
+			
+			% 2D, one slice
+			if size(V_in,3) == 1
+				V_out = this.apply_2D(V_in);
+						
+			% 3D, slice by slice
+			elseif ~strcmp(this.full3d,'yes')
+				V_out = zeros(size(V_in));
+				for slice = 1:size(V_in,3)
+					V_out(:,:,slice) = this.apply_2D(V_in(:,:,slice)); 
+				end
+			
+			% 3D, full
+			else
+				V_out = this.apply_3D(V_in);
+			end
+			
+		end	
+			
+	end
+		
+	%----------------------------------------------------------------------
+	methods (Access=protected)
+		
+		%------------------------------------------------------------------
+		function V_out = apply_2D(this, V_in)
+			
+			vol_geom = astra_create_vol_geom(size(V_in));
+			in_id = astra_mex_data2d('create', '-vol', vol_geom, V_in);
+			out_id = astra_mex_data2d('create', '-vol', vol_geom, 0);
+
+			cfg = astra_struct('DARTSMOOTHING_CUDA');
+			cfg.InDataId = in_id;
+			cfg.OutDataId = out_id;
+			cfg.Intensity = this.b;
+			cfg.option.GPUindex = this.gpu_core;
+			
+			alg_id = astra_mex_algorithm('create',cfg);	
+			astra_mex_algorithm('iterate',alg_id,1);
+			V_out = astra_mex_data2d('get', out_id);
+		
+			astra_mex_algorithm('delete', alg_id);
+			astra_mex_data2d('delete', in_id, out_id);
+			
+			
+		end
+		
+		%------------------------------------------------------------------
+		function I_out = apply_3D(this, I_in)
+			
+			r = this.radius;
+			w = 2 * r + 1;
+			
+			% Set Kernel
+			K = ones(w,w,w) * this.b / (w.^3-1); % edges
+			K(r+1,r+1,r+1) = 1 - this.b; % center
+			
+			% output window
+			I_out = zeros(size(I_in,1)+w-1, size(I_in,2)+w-1, size(I_in,3)+w-1);
+
+			% blur convolution
+			for s = -r:r 
+				for t = -r:r 
+					for u = -r:r 
+						I_out(1+r+s:end-r+s, 1+r+t:end-r+t, 1+r+u:end-r+u) = I_out(1+r+s:end-r+s, 1+r+t:end-r+t, 1+r+u:end-r+u) + K(r+1+s, r+1+t, r+1+u) * I_in;
+					end
+				end
+			end
+			
+			% shrink output window
+			I_out = I_out(1+r:end-r, 1+r:end-r, 1+r:end-r);
+			
+		end		
+		%------------------------------------------------------------------
+		
+	end
+	
+end
+
diff --git a/matlab/algorithms/DART/StatisticsDefault.m b/matlab/algorithms/DART/StatisticsDefault.m
new file mode 100644
index 0000000..7822c5f
--- /dev/null
+++ b/matlab/algorithms/DART/StatisticsDefault.m
@@ -0,0 +1,72 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef StatisticsDefault < matlab.mixin.Copyable
+	
+	% Default policy class for statistics for DART.	
+	
+	properties (Access=public)
+		pixel_error		= 'yes';		% SETTING: Store pixel error? {'yes','no'}
+		proj_diff		= 'yes';		% SETTING: Store projection difference? {'yes','no'}
+		timing			= 'yes';		% SETTING: Store timings? {'yes','no'}
+	end
+
+	
+	methods (Access=public)
+		
+		%------------------------------------------------------------------
+		function stats = apply(this, DART)
+			% Applies statistics.
+			% >> stats = DART.statistics.apply(DART);				
+			
+			stats = DART.stats;
+			
+			% timing
+			if strcmp(this.timing, 'yes')
+				stats.timing(DART.iterationcount) = toc(DART.start_tic);
+			end
+			
+			% pixel error
+			if strcmp(this.pixel_error, 'yes') && isfield(DART.base,'phantom')
+				[stats.rnmp, stats.nmp] = compute_rnmp(DART.base.phantom, DART.S);
+				stats.rnmp_hist(DART.iterationcount) = stats.rnmp;
+				stats.nmp_hist(DART.iterationcount) = stats.nmp;
+			end
+			
+			% projection difference
+			if strcmp(this.proj_diff, 'yes') 
+				new_sino = DART.tomography.createForwardProjection(DART, DART.S);
+				stats.proj_diff = sum((new_sino(:) - DART.base.sinogram(:)) .^2 ) ./ (sum(DART.base.sinogram(:)) );
+				stats.proj_diff_hist(DART.iterationcount) = stats.proj_diff;
+			end
+			
+		end
+		
+		%------------------------------------------------------------------
+		function s = tostring(~, stats)
+			% To string.
+			% >> stats = DART.statistics.apply(stats);				
+			
+			s = '';
+			if isfield(stats, 'nmp')
+				s = sprintf('%s [%d]', s, stats.nmp); 
+			end
+			if isfield(stats, 'proj_diff')
+				s = sprintf('%s {%0.2d}', s, stats.proj_diff); 
+			end			
+			
+		end
+		%------------------------------------------------------------------
+		
+	end
+	
+end
+
diff --git a/matlab/algorithms/DART/TomographyDefault.m b/matlab/algorithms/DART/TomographyDefault.m
new file mode 100644
index 0000000..4db3905
--- /dev/null
+++ b/matlab/algorithms/DART/TomographyDefault.m
@@ -0,0 +1,73 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef TomographyDefault < IterativeTomography
+
+	% Policy class for tomography for DART.
+
+	%----------------------------------------------------------------------
+	properties (Access=public)
+		t					= 5;		% SETTING: # ARMiterations, each DART iteration.
+		t0					= 100;		% SETTING: # ARM iterations at DART initialization.
+	end
+	%----------------------------------------------------------------------
+
+	methods
+
+		%------------------------------------------------------------------
+		function settings = getsettings(this)
+			% Returns a structure containing all settings of this object.
+			% >> settings = DART.tomography.getsettings();
+%			settings = getsettings@IterativeTomography();
+			settings.t					= this.t;
+			settings.t0					= this.t0;
+		end
+
+		%------------------------------------------------------------------
+		function initialize(this, DART)
+			% Initializes this object. 
+			% >> DART.tomography.initialize();
+			this.proj_geom = DART.base.proj_geom;
+			this.initialize@IterativeTomography();
+		end
+		
+		%------------------------------------------------------------------
+		function P = createForwardProjection(this, ~, volume)
+			% Compute forward projection. 
+			% >> DART.tomography.createForwardProjection(DART, volume);
+			P = this.project_c(volume);
+		end
+		
+		%------------------------------------------------------------------
+		function I = createReconstruction(this, ~, sinogram, V0, mask)
+			% Compute reconstruction (with mask).
+			% >> DART.tomography.createReconstruction(DART, sinogram, V0, mask);			
+			if strcmp(this.inner_circle,'yes')
+				mask = ROIselectfull(mask, size(mask,1));
+			end
+			I = this.reconstruct_c(sinogram, V0, mask, this.t);
+		end	
+		
+		%------------------------------------------------------------------
+		function I = createInitialReconstruction(this, ~, sinogram)
+			% Compute reconstruction (initial).
+			% >> DART.tomography.createInitialReconstruction(DART, sinogram);						
+			I = this.reconstruct_c(sinogram, [], [], this.t0);
+			if strcmp(this.inner_circle,'yes')
+				I = ROIselectfull(I, size(I,1));
+			end
+		end			
+		%------------------------------------------------------------------
+	
+	end		
+		
+end
+
diff --git a/matlab/algorithms/DART/TomographyDefault3D.m b/matlab/algorithms/DART/TomographyDefault3D.m
new file mode 100644
index 0000000..2be1b17
--- /dev/null
+++ b/matlab/algorithms/DART/TomographyDefault3D.m
@@ -0,0 +1,73 @@
+% This file is part of the
+% All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA-Toolbox")
+%
+% Copyright: iMinds-Vision Lab, University of Antwerp
+% License: Open Source under GPLv3
+% Contact: mailto:astra@ua.ac.be
+% Website: http://astra.ua.ac.be
+%
+% Author of this DART Algorithm: Wim van Aarle
+
+
+classdef TomographyDefault3D < IterativeTomography3D
+
+	% Policy class for 3D tomography for DART.
+
+	%----------------------------------------------------------------------
+	properties (Access=public)
+		t					= 5;		% SETTING: # ARMiterations, each DART iteration.
+		t0					= 100;		% SETTING: # ARM iterations at DART initialization.
+	end
+	%----------------------------------------------------------------------
+
+	methods
+
+		%------------------------------------------------------------------
+		function settings = getsettings(this)
+			% Returns a structure containing all settings of this object.
+			% >> settings = DART.tomography.getsettings();
+%			settings = getsettings@IterativeTomography();
+			settings.t					= this.t;
+			settings.t0					= this.t0;
+		end
+
+		%------------------------------------------------------------------
+		function initialize(this, DART)
+			% Initializes this object. 
+			% >> DART.tomography.initialize();
+			this.proj_geom = DART.base.proj_geom;
+			this.initialize@IterativeTomography3D();
+		end
+		
+		%------------------------------------------------------------------
+		function P = createForwardProjection(this, ~, volume)
+			% Compute forward projection. 
+			% >> DART.tomography.createForwardProjection(DART, volume);
+			P = this.project_c(volume);
+		end
+		
+		%------------------------------------------------------------------
+		function I = createReconstruction(this, ~, sinogram, V0, mask)
+			% Compute reconstruction (with mask).
+			% >> DART.tomography.createReconstruction(DART, sinogram, V0, mask);			
+			if strcmp(this.inner_circle,'yes')
+				mask = ROIselectfull(mask, size(mask,1));
+			end
+			I = this.reconstruct_c(sinogram, V0, mask, this.t);
+		end	
+		
+		%------------------------------------------------------------------
+		function I = createInitialReconstruction(this, ~, sinogram)
+			% Compute reconstruction (initial).
+			% >> DART.tomography.createInitialReconstruction(DART, sinogram);						
+			I = this.reconstruct_c(sinogram, [], [], this.t0);
+			if strcmp(this.inner_circle,'yes')
+				I = ROIselectfull(I, size(I,1));
+			end
+		end			
+		%------------------------------------------------------------------
+	
+	end		
+		
+end
+
diff --git a/matlab/algorithms/DART/examples/cylinders.png b/matlab/algorithms/DART/examples/cylinders.png
new file mode 100644
index 0000000..8d1c0b2
--- /dev/null
+++ b/matlab/algorithms/DART/examples/cylinders.png
diff --git a/matlab/algorithms/DART/examples/example1.m b/matlab/algorithms/DART/examples/example1.m
new file mode 100644
index 0000000..daa3ce8
--- /dev/null
+++ b/matlab/algorithms/DART/examples/example1.m
@@ -0,0 +1,79 @@
+clear all;
+
+addpath('..');
+
+%
+% Example 1: parallel beam, three slices.
+%
+
+% Configuration
+proj_count = 20;
+slice_count = 3;
+dart_iterations = 20;
+filename = 'cylinders.png';
+outdir = './';
+prefix = 'example1';
+rho = [0, 1];
+tau = 0.5;
+gpu_core = 0;
+
+% Load phantom.
+I = double(imread(filename)) / 255;
+
+% Create projection and volume geometries.
+det_count = size(I, 1);
+angles = linspace(0, pi - pi / proj_count, proj_count);
+proj_geom = astra_create_proj_geom('parallel3d', 1, 1, slice_count, det_count, angles);
+vol_geom = astra_create_vol_geom(det_count, det_count, 1);
+
+% Create sinogram.
+[sinogram_id, sinogram] = astra_create_sino3d_cuda(I, proj_geom, vol_geom);
+astra_mex_data3d('delete', sinogram_id);
+
+%
+% DART
+%
+
+base.sinogram = sinogram;
+base.proj_geom = proj_geom;
+
+D                      = DARTalgorithm(base);
+
+D.tomography           = TomographyDefault3D();
+D.tomography.t0        = 100;
+D.tomography.t         = 10;
+D.tomography.method    = 'SIRT3D_CUDA';
+D.tomography.gpu_core  = gpu_core;
+D.tomography.use_minc  = 'yes';
+% D.tomography.maxc      = 0.003; % Not a sensible value, just for demonstration.
+
+D.segmentation.rho     = rho;
+D.segmentation.tau     = tau;
+
+D.smoothing.b          = 0.1;
+D.smoothing.full3d     = 'yes';
+D.smoothing.gpu_core   = gpu_core;
+ 
+D.masking.random       = 0.1;
+D.masking.conn         = 6;
+D.masking.gpu_core     = gpu_core;
+
+D.output.directory     = outdir;
+D.output.pre           = [prefix '_'];
+D.output.save_images   = 'no';
+D.output.save_results  = {'stats', 'settings', 'S', 'V'};
+D.output.save_interval = dart_iterations;
+D.output.verbose       = 'yes';
+
+D.statistics.proj_diff = 'no';
+
+D.initialize();
+
+disp([D.output.directory D.output.pre]);
+
+D.iterate(dart_iterations);
+
+% Convert middle slice of final iteration to png.
+load([outdir '/' prefix '_results_' num2str(dart_iterations) '.mat']);
+imwritesc(D.S(:, :, round(slice_count / 2)), [outdir '/' prefix '_slice_2_S.png']);
+imwritesc(D.V(:, :, round(slice_count / 2)), [outdir '/' prefix '_slice_2_V.png']);
diff --git a/matlab/algorithms/DART/examples/example2.m b/matlab/algorithms/DART/examples/example2.m
new file mode 100644
index 0000000..8ee8cba
--- /dev/null
+++ b/matlab/algorithms/DART/examples/example2.m
@@ -0,0 +1,80 @@
+clear all;
+
+addpath('..');
+
+%
+% Example 2: cone beam, full cube.
+%
+
+% Configuration
+det_count = 128;
+proj_count = 45;
+slice_count = det_count;
+dart_iterations = 20;
+outdir = './';
+prefix = 'example2';
+rho = [0 0.5 1];
+tau = [0.25 0.75];
+gpu_core = 0;
+
+% Create phantom.
+% I = phantom3d([1 0.9 0.9 0.9 0 0 0 0 0 0; -0.5 0.8 0.8 0.8 0 0 0 0 0 0; -0.5 0.3 0.3 0.3 0 0 0 0 0 0], det_count);
+% save('phantom3d', 'I');
+load('phantom3d'); % Loads I.
+
+% Create projection and volume geometries.
+angles = linspace(0, pi - pi / proj_count, proj_count);
+proj_geom = astra_create_proj_geom('cone', 1, 1, slice_count, det_count, angles, 500, 0);
+vol_geom = astra_create_vol_geom(det_count, det_count, slice_count);
+
+% Create sinogram.
+[sinogram_id, sinogram] = astra_create_sino3d_cuda(I, proj_geom, vol_geom);
+astra_mex_data3d('delete', sinogram_id);
+
+%
+% DART
+%
+
+base.sinogram = sinogram;
+base.proj_geom = proj_geom;
+
+D                      = DARTalgorithm(base);
+
+D.tomography           = TomographyDefault3D();
+D.tomography.t0        = 100;
+D.tomography.t         = 10;
+D.tomography.method    = 'SIRT3D_CUDA';
+D.tomography.gpu_core  = gpu_core;
+D.tomography.use_minc  = 'yes';
+% D.tomography.maxc      = 0.003; % Not a sensible value, just for demonstration.
+
+D.segmentation.rho     = rho;
+D.segmentation.tau     = tau;
+
+D.smoothing.b          = 0.1;
+D.smoothing.full3d     = 'yes';
+D.smoothing.gpu_core   = gpu_core;
+ 
+D.masking.random       = 0.1;
+D.masking.conn         = 6;
+D.masking.gpu_core     = gpu_core;
+
+D.output.directory     = outdir;
+D.output.pre           = [prefix '_'];
+D.output.save_images   = 'no';
+D.output.save_results  = {'stats', 'settings', 'S', 'V'};
+D.output.save_interval = dart_iterations;
+D.output.verbose       = 'yes';
+
+D.statistics.proj_diff = 'no';
+
+D.initialize();
+
+disp([D.output.directory D.output.pre]);
+
+D.iterate(dart_iterations);
+
+% Convert middle slice of final iteration to png.
+load([outdir '/' prefix '_results_' num2str(dart_iterations) '.mat']);
+imwritesc(D.S(:, :, round(slice_count / 2)), [outdir '/' prefix '_slice_2_S.png']);
+imwritesc(D.V(:, :, round(slice_count / 2)), [outdir '/' prefix '_slice_2_V.png']);
diff --git a/matlab/algorithms/DART/examples/example3.m b/matlab/algorithms/DART/examples/example3.m
new file mode 100644
index 0000000..f6e360e
--- /dev/null
+++ b/matlab/algorithms/DART/examples/example3.m
@@ -0,0 +1,79 @@
+clear all;
+
+addpath('..');
+
+%
+% Example 3: parallel beam, 2D
+%
+
+% Configuration
+proj_count = 30;
+dart_iterations = 20;
+filename = 'cylinders.png';
+outdir = './';
+prefix = 'example3';
+rho = [0, 1];
+tau = 0.5;
+gpu_core = 0;
+
+% Load phantom.
+I = double(imread(filename)) / 255;
+
+% Create projection and volume geometries.
+det_count = size(I, 1);
+angles = linspace(0, pi - pi / proj_count, proj_count);
+proj_geom = astra_create_proj_geom('parallel', 1, det_count, angles);
+vol_geom = astra_create_vol_geom(det_count, det_count);
+
+% Create sinogram.
+[sinogram_id, sinogram] = astra_create_sino_cuda(I, proj_geom, vol_geom);
+astra_mex_data2d('delete', sinogram_id);
+
+%
+% DART
+%
+
+base.sinogram = sinogram;
+base.proj_geom = proj_geom;
+
+D                      = DARTalgorithm(base);
+
+%D.tomography           = TomographyDefault3D();
+D.tomography.t0        = 100;
+D.tomography.t         = 10;
+D.tomography.method    = 'SIRT_CUDA';
+D.tomography.proj_type = 'strip';
+D.tomography.gpu_core  = gpu_core;
+D.tomography.use_minc  = 'yes';
+% D.tomography.maxc      = 0.003; % Not a sensible value, just for demonstration.
+
+D.segmentation.rho     = rho;
+D.segmentation.tau     = tau;
+
+D.smoothing.b          = 0.1;
+D.smoothing.full3d     = 'yes';
+D.smoothing.gpu_core   = gpu_core;
+ 
+D.masking.random       = 0.1;
+D.masking.conn         = 6;
+D.masking.gpu_core     = gpu_core;
+
+D.output.directory     = outdir;
+D.output.pre           = [prefix '_'];
+D.output.save_images   = 'no';
+D.output.save_results  = {'stats', 'settings', 'S', 'V'};
+D.output.save_interval = dart_iterations;
+D.output.verbose       = 'yes';
+
+D.statistics.proj_diff = 'no';
+
+D.initialize();
+
+disp([D.output.directory D.output.pre]);
+
+D.iterate(dart_iterations);
+
+% Convert output of final iteration to png.
+load([outdir '/' prefix '_results_' num2str(dart_iterations) '.mat']);
+imwritesc(D.S, [outdir '/' prefix '_S.png']);
+imwritesc(D.V, [outdir '/' prefix '_V.png']);
diff --git a/matlab/algorithms/DART/examples/phantom3d.mat b/matlab/algorithms/DART/examples/phantom3d.mat
new file mode 100644
index 0000000..6d70c16
--- /dev/null
+++ b/matlab/algorithms/DART/examples/phantom3d.mat