Merge branch 'master' into python-interface

21 files changed, 1081 insertions, 179 deletions

diff --git a/README.md b/README.md
new file mode 100644
index 0000000..ec86136
--- /dev/null
+++ b/README.md
@@ -0,0 +1,82 @@
+# The ASTRA Toolbox
+
+The ASTRA Toolbox is a MATLAB toolbox of high-performance GPU primitives for 2D and 3D tomography.
+
+We support 2D parallel and fan beam geometries, and 3D parallel and cone beam.  All of them have highly flexible source/detector positioning.
+
+A large number of 2D and 3D algorithms are available, including FBP, SIRT, SART, CGLS.
+
+The basic forward and backward projection operations are GPU-accelerated, and directly callable from MATLAB to enable building new algorithms.
+
+
+
+## Documentation / samples
+
+See the matlab code samples in samples/ and on http://sf.net/projects/astra-toolbox .
+
+
+## Installation instructions
+
+### Windows, binary
+
+Add the mex and tools subdirectories to your matlab path.
+
+### Linux, from source
+
+Requirements: g++, boost, CUDA (driver+toolkit), matlab
+
+```
+cd build/linux
+./configure --with-cuda=/usr/local/cuda \
+            --with-matlab=/usr/local/MATLAB/R2012a \
+            --prefix=/usr/local/astra
+make
+make install
+```
+Add /usr/local/astra/lib to your LD_LIBRARY_PATH.
+Add /usr/local/astra/matlab and its subdirectories (tools, mex)
+  to your matlab path.
+
+
+NB: Each matlab version only supports a specific range of g++ versions.
+Despite this, if you have a newer g++ and if you get errors related to missing
+GLIBCXX_3.4.xx symbols, it is often possible to work around this requirement
+by deleting the version of libstdc++ supplied by matlab in
+MATLAB_PATH/bin/glnx86 or MATLAB_PATH/bin/glnxa64 (at your own risk).
+
+
+### Windows, from source using Visual Studio 2008
+
+Requirements: Visual Studio 2008, boost, CUDA (driver+toolkit), matlab.
+Note that a .zip with all required (and precompiled) boost files is
+  available from our website.
+
+Set the environment variable MATLAB_ROOT to your matlab install location.
+Open astra_vc08.sln in Visual Studio.
+Select the appropriate solution configuration.
+  (typically Release_CUDA|win32 or Release_CUDA|x64)
+Build the solution.
+Install by copying AstraCuda32.dll or AstraCuda64.dll from bin/ and
+  all .mexw32 or .mexw64 files from bin/Release_CUDA or bin/Debug_CUDA
+  and the entire matlab/tools directory to a directory to be added to
+  your matlab path.
+
+
+## References
+
+If you use parallel beam GPU code for your research, we would appreciate it if you would refer to the following paper:
+
+W. J. Palenstijn, K J. Batenburg, and J. Sijbers, "Performance improvements for iterative electron tomography reconstruction using graphics processing units (GPUs)", Journal of Structural Biology, vol. 176, issue 2, pp. 250-253, 2011.
+
+## License
+
+The ASTRA Toolbox is open source under the GPLv3 license.
+
+## Contact
+
+email: astra@uantwerpen.be
+website: http://sf.net/projects/astra-toolbox
+
+Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+           2014-2015, CWI, Amsterdam
+           http://visielab.uantwerpen.be/ and http://www.cwi.nl/
\ No newline at end of file
diff --git a/astra_vc11.vcxproj b/astra_vc11.vcxproj
index 1288d35..482bb1b 100644
--- a/astra_vc11.vcxproj
+++ b/astra_vc11.vcxproj
@@ -194,7 +194,7 @@
       <OptimizeReferences>true</OptimizeReferences>
       <OutputFile>bin\x64\Release_CUDA\AstraCuda64.dll</OutputFile>
       <AdditionalLibraryDirectories>lib\x64;$(CUDA_LIB_PATH);$(NVSDKCUDA_ROOT)\common\lib;$(NVSDKCOMPUTE_ROOT)/C/common/lib;%(AdditionalLibraryDirectories)</AdditionalLibraryDirectories>
-      <AdditionalDependencies>cudart.lib;FreeImage.lib;cufft.lib;%(AdditionalDependencies)</AdditionalDependencies>
+      <AdditionalDependencies>cudart.lib;cufft.lib;%(AdditionalDependencies)</AdditionalDependencies>
       <IgnoreAllDefaultLibraries>false</IgnoreAllDefaultLibraries>
     </Link>
     <CudaCompile>
@@ -436,4 +436,4 @@
   <ImportGroup Label="ExtensionTargets">
     <Import Project="$(VCTargetsPath)\BuildCustomizations\CUDA 5.5.targets" />
   </ImportGroup>
-</Project>
+</Project>
\ No newline at end of file
diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in
index fe1ba91..6620446 100644
--- a/build/linux/Makefile.in
+++ b/build/linux/Makefile.in
@@ -25,10 +25,11 @@ VPATH=../..
 CPPFLAGS=@SAVED_CPPFLAGS@
 CXXFLAGS=@SAVED_CXXFLAGS@
 LDFLAGS=@SAVED_LDFLAGS@
+LIBS=@SAVED_LIBS@
 
 CPPFLAGS+=-I../.. -I../../include -I../../lib/include/rapidxml
 CXXFLAGS+=-g -O3 -Wall -Wshadow
-LIBS=-lpthread -lrt
+LIBS+=-lpthread
 LDFLAGS+=-g
 
 ifeq ($(cuda),yes)
diff --git a/build/linux/acinclude.m4 b/build/linux/acinclude.m4
index 4ff9e4b..713b5d3 100644
--- a/build/linux/acinclude.m4
+++ b/build/linux/acinclude.m4
@@ -60,8 +60,14 @@ AC_DEFUN([ASTRA_RUN_STOREOUTPUT],[{
   test $ac_status = 0;
  }])
 
-dnl ASTRA_RUN(command)
-AC_DEFUN([ASTRA_RUN],[ASTRA_RUN_STOREOUTPUT($1,/dev/null)])
+dnl ASTRA_RUN_LOGOUTPUT(command)
+AC_DEFUN([ASTRA_RUN_LOGOUTPUT],[{
+  AS_ECHO(["$as_me:${as_lineno-$LINENO}: $1"]) >&AS_MESSAGE_LOG_FD
+  ( $1 ) >&AS_MESSAGE_LOG_FD 2>&1
+  ac_status=$?
+  AS_ECHO(["$as_me:${as_lineno-$LINENO}: \$? = $ac_status"]) >&AS_MESSAGE_LOG_FD
+  test $ac_status = 0;
+ }])
 
 
 
@@ -79,9 +85,9 @@ ASTRA_RUN_STOREOUTPUT([$NVCC -c -o conftest.o conftest.cu $$2],conftest.nvcc.out
   $1="no"
   # Check if hack for gcc 4.4 helps
   if grep -q __builtin_stdarg_start conftest.nvcc.out; then
+    AS_ECHO(["$as_me:${as_lineno-$LINENO}: Trying CUDA hack for gcc 4.4"]) >&AS_MESSAGE_LOG_FD
     NVCC_OPT="-Xcompiler -D__builtin_stdarg_start=__builtin_va_start"
-
-    ASTRA_RUN([$NVCC -c -o conftest.o conftest.cu $$2 $NVCC_OPT]) && {
+    ASTRA_RUN_LOGOUTPUT([$NVCC -c -o conftest.o conftest.cu $$2 $NVCC_OPT]) && {
       $1="yes"
       $2="$$2 $NVCC_OPT"
     }
@@ -94,6 +100,7 @@ fi
 rm -f conftest.cu conftest.o conftest.nvcc.out
 ])
 
+
 dnl ASTRA_FIND_NVCC_ARCHS(archs-to-try,cppflags-to-extend,output-list)
 dnl Architectures should be of the form 10,20,30,35,
 dnl and should be in order. The last accepted one will be used for PTX output.
diff --git a/build/linux/autogen.sh b/build/linux/autogen.sh
index c856793..544fdeb 100755
--- a/build/linux/autogen.sh
+++ b/build/linux/autogen.sh
@@ -12,9 +12,12 @@ if test $? -ne 0; then
   exit 1
 fi
 
-libtoolize --install --force > /dev/null 2>&1
+case `uname` in Darwin*) LIBTOOLIZEBIN=glibtoolize ;;
+  *) LIBTOOLIZEBIN=libtoolize ;; esac
+
+$LIBTOOLIZEBIN --install --force > /dev/null 2>&1
 if test $? -ne 0; then
-  libtoolize --force
+  $LIBTOOLIZEBIN --force
   if test $? -ne 0; then
     echo "Error running libtoolize"
     exit 1
diff --git a/build/linux/configure.ac b/build/linux/configure.ac
index 7ad03c3..f4cc82e 100644
--- a/build/linux/configure.ac
+++ b/build/linux/configure.ac
@@ -31,6 +31,7 @@ LT_INIT([disable-static])
 SAVED_CPPFLAGS="$CPPFLAGS"
 SAVED_CXXFLAGS="$CXXFLAGS"
 SAVED_LDFLAGS="$LDFLAGS"
+SAVED_LIBS="$LIBS"
 
 AC_CANONICAL_BUILD
 AC_CANONICAL_HOST
@@ -54,7 +55,7 @@ AC_MSG_CHECKING([for boost-unit-test-framework])
 ASTRA_CHECK_BOOST_UNIT_TEST_FRAMEWORK(-lboost_unit_test_framework-mt, BOOSTUTF=yes_mt, BOOSTUTF=no)
 if test x$BOOSTUTF = xno; then
   ASTRA_CHECK_BOOST_UNIT_TEST_FRAMEWORK(-lboost_unit_test_framework, BOOSTUTF=yes, BOOSTUTF=no)
-  if test x$BOOSTTHREAD = xno; then
+  if test x$BOOSTUTF = xno; then
     AC_MSG_RESULT(no)
     AC_MSG_ERROR([No boost-unit-test-framework library found])
   else
@@ -65,48 +66,41 @@ else
   AC_MSG_RESULT([yes, libboost_unit_test_framework-mt])
   LIBS_BOOSTUTF="-lboost_unit_test_framework-mt"
 fi
-# TODO: do something with the result
 
 
 # nvcc, cuda
 
 AC_ARG_WITH(cuda, [[  --with-cuda=path        path of CUDA SDK (optional)]],,)
 
-NVCC_PATH=$PATH
-if test x"$with_cuda" != x; then
-  NVCC_PATH="$with_cuda/bin:$NVCC_PATH"
+if test x"$with_cuda" != xno; then
+  NVCC_PATH=$PATH
+  if test x"$with_cuda" != x -a x"$with_cuda" != xyes; then
+    NVCC_PATH="$with_cuda/bin:$NVCC_PATH"
+  fi
+  AC_PATH_PROG([NVCC], [nvcc], [no], [$NVCC_PATH])
+else
+  NVCC=no
 fi
-AC_PATH_PROG([NVCC], [nvcc], [no], [$NVCC_PATH])
-# TODO: do something with the result
 
 HAVECUDA=no
 if test x"$NVCC" != xno; then
   HAVECUDA=yes
   BACKUP_CUDA_LDFLAGS="$LDFLAGS"
-  if test x"$with_cuda" != x; then
-    LDFLAGS_CUDA="-L$with_cuda/lib"
+  if test x"$with_cuda" != x -a x"$with_cuda" != xyes; then
+    case $host_cpu in
+    x86_64)
+      LDFLAGS_CUDA="-L$with_cuda/lib64"
+      ;;
+    *)
+      LDFLAGS_CUDA="-L$with_cuda/lib"
+      ;;
+    esac
     CPPFLAGS_CUDA="-I$with_cuda/include"
     LDFLAGS="$LDFLAGS $LDFLAGS_CUDA"
   fi
   AC_CHECK_LIB(cudart,cudaMalloc, ,HAVECUDA=no)
   AC_CHECK_LIB(cufft,cufftPlan1d, ,HAVECUDA=no)
 
-  if test x"$HAVECUDA" = xno; then
-    # try lib64 instead of lib
-
-    HAVECUDA=yes
-    LDFLAGS="$BACKUP_CUDA_LDFLAGS"
-
-    # prevent cached values from being used
-    unset ac_cv_lib_cudart_cudaMalloc
-    unset ac_cv_lib_cufft_cufftPlan1d
-
-    LDFLAGS_CUDA="-L$with_cuda/lib64"
-    LDFLAGS="$LDFLAGS $LDFLAGS_CUDA"
-    AC_CHECK_LIB(cudart,cudaMalloc, ,HAVECUDA=no)
-    AC_CHECK_LIB(cufft,cufftPlan1d, ,HAVECUDA=no)
-  fi
-
   LDFLAGS="$BACKUP_CUDA_LDFLAGS"
   unset BACKUP_CUDA_LDFLAGS
   # TODO: check for cuda headers?
@@ -115,11 +109,11 @@ if test x"$NVCC" != xno; then
 fi
 
 NVCCFLAGS=""
-AC_MSG_CHECKING([if nvcc works])
 if test x"$HAVECUDA" = xyes; then
+  AC_MSG_CHECKING([if nvcc works])
   ASTRA_CHECK_NVCC(HAVECUDA,NVCCFLAGS)
+  AC_MSG_RESULT($HAVECUDA)
 fi
-AC_MSG_RESULT($HAVECUDA)
 
 AC_ARG_WITH(cuda_compute, [[  --with-cuda-compute=archs  comma separated list of CUDA compute models (optional)]],,)
 if test x"$HAVECUDA" = xyes; then
@@ -154,7 +148,7 @@ if test x"$with_matlab" != x; then
     AC_SUBST(MEX)
     MATLAB_ROOT="$with_matlab"
     AC_SUBST(MATLAB_ROOT)
-
+    # TODO: maybe catch mex warnings
     ASTRA_CHECK_MEX_SUFFIX([mexa64 mexglx mexmaci64 mexmaci],[MEXSUFFIX])
     if test x$MEXSUFFIX = x; then
       AC_MSG_FAILURE([Unable to determine matlab mex suffix])
@@ -196,7 +190,14 @@ AC_SUBST(HAVEPYTHON)
 AC_SUBST(SAVED_CPPFLAGS)
 AC_SUBST(SAVED_CXXFLAGS)
 AC_SUBST(SAVED_LDFLAGS)
-
-
+AC_SUBST(SAVED_LIBS)
 AC_CONFIG_FILES([Makefile])
 AC_OUTPUT
+
+echo
+echo "Summary of ASTRA Toolbox build options:"
+echo "  CUDA   : $HAVECUDA"
+echo "  Matlab : $HAVEMATLAB"
+echo
+echo "  prefix : $prefix"
+echo
diff --git a/include/astra/AsyncAlgorithm.h b/include/astra/AsyncAlgorithm.h
index 2d5d31e..a3157fc 100644
--- a/include/astra/AsyncAlgorithm.h
+++ b/include/astra/AsyncAlgorithm.h
@@ -32,14 +32,12 @@ $Id$
 #include "Config.h"
 #include "Algorithm.h"
 
-#ifdef __linux__
-#define USE_PTHREADS
+#ifdef USE_PTHREADS
 #include <pthread.h>
 #else
 #include <boost/thread.hpp>
 #endif
 
-
 namespace astra {
 	
 /**
@@ -75,10 +73,6 @@ public:
 	 */
 	virtual void run(int _iNrIterations = 0);
 
-	/** Wait for thread to complete and delete thread. 
-	 */
-	virtual void timedJoin(int _milliseconds);
-
 	/** Return pointer to the wrapped algorithm. 
 	 */
 	CAlgorithm* getWrappedAlgorithm() { return m_pAlg; }
diff --git a/include/astra/Globals.h b/include/astra/Globals.h
index fdeaa23..9c8ddfb 100644
--- a/include/astra/Globals.h
+++ b/include/astra/Globals.h
@@ -306,4 +306,11 @@ _AstraExport inline bool cudaEnabled() { return false; }
 
 #endif
 
+//----------------------------------------------------------------------------------------
+// use pthreads on Linux and OSX
+#if defined(__linux__) || defined(__MACH__)
+#define USE_PTHREADS
+#endif
+
+
 #endif
diff --git a/matlab/algorithms/DART/DARTalgorithm.m b/matlab/algorithms/DART/DARTalgorithm.m
index 85a3ca0..fc707dd 100644
--- a/matlab/algorithms/DART/DARTalgorithm.m
+++ b/matlab/algorithms/DART/DARTalgorithm.m
@@ -9,13 +9,8 @@
 %--------------------------------------------------------------------------
 
 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)
 	
@@ -78,7 +73,6 @@ classdef DARTalgorithm < matlab.mixin.Copyable
 				error('invalid arguments')
 			end
 		end	
-	
 		
 		%------------------------------------------------------------------
 		function D = deepcopy(this)
@@ -100,7 +94,6 @@ classdef DARTalgorithm < matlab.mixin.Copyable
 			
 			% Initialize tomography part
 			if ~this.tomography.initialized
-				this.tomography.sinogram = this.base.sinogram;
 				this.tomography.proj_geom = this.base.proj_geom;	
 				this.tomography.initialize();
 			end
@@ -110,7 +103,7 @@ classdef DARTalgorithm < matlab.mixin.Copyable
 				this.V0 = this.base.V0;
 			else
 				this.output.pre_initial_iteration(this);
-				this.V0 = this.tomography.reconstruct2(this.base.sinogram, [], this.t0);
+				this.V0 = this.tomography.reconstruct(this.base.sinogram, this.t0);
 				this.output.post_initial_iteration(this);
 			end
 			this.V = this.V0;
@@ -163,7 +156,7 @@ classdef DARTalgorithm < matlab.mixin.Copyable
 				this.R = this.base.sinogram - this.tomography.project(F);
 				
 				% ART update part
-				this.V = this.tomography.reconstruct2_mask(this.R, this.V, this.Mask, this.t);
+				this.V = this.tomography.reconstruct_mask(this.R, this.V, this.Mask, this.t);
 
 				% blur
 				this.V = this.smoothing.apply(this, this.V);
diff --git a/matlab/algorithms/DART/IterativeTomography.m b/matlab/algorithms/DART/IterativeTomography.m
new file mode 100644
index 0000000..3875e6b
--- /dev/null
+++ b/matlab/algorithms/DART/IterativeTomography.m
@@ -0,0 +1,435 @@
+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)
+		proj_geom			= [];				% DATA: Projection geometry.
+		vol_geom			= [];				% DATA: Volume geometry.
+	end
+	%----------------------------------------------------------------------
+	properties (SetAccess=private, GetAccess=public)
+		initialized			= 0;				% Is this object initialized?
+	end
+	%----------------------------------------------------------------------
+	properties (SetAccess=protected, GetAccess=protected)
+		cfg_base			= struct();			% PROTECTED: base configuration structure for the reconstruction algorithm.
+		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')	
+    end
+	%----------------------------------------------------------------------
+
+	methods (Access=public)
+		
+		%------------------------------------------------------------------
+		function this = IterativeTomography(varargin)
+			% Constructor
+			% >> tomography = IterativeTomography(proj_geom);
+			% >> tomography = IterativeTomography(proj_geom, vol_geom);
+			
+			% Input: IterativeTomography(proj_geom)
+			if nargin == 1
+				this.proj_geom = varargin{1};
+			
+			% Input: IterativeTomography(proj_geom, vol_geom)
+			elseif nargin == 2
+				this.proj_geom = varargin{1};
+				this.vol_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
+			if this.superresolution > 1
+				this.proj_geom_sr = astra_geom_superresolution(this.proj_geom, this.superresolution);
+			else
+				this.proj_geom_sr = this.proj_geom;
+			end
+				
+			% 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;
+			end
+			this.cfg_base.option.DetectorSuperSampling = this.superresolution;
+			if strcmp(this.gpu,'yes')
+				this.cfg_base.option.GPUindex = this.gpu_core;
+			end
+			if this.use_minc
+				this.cfg_base.option.MinConstraint = 0;
+			end
+		
+			this.initialized = 1;
+			ok = this.initialized;
+		end
+		
+		%------------------------------------------------------------------
+		function projections = project(this, volume)
+			% Compute forward projection.  
+			% >> projections = tomography.project(volume);
+			
+			if ~this.initialized 
+				this.initialize();
+			end
+
+			% project			
+			projections = this.project_c(volume);
+		end
+		
+		%------------------------------------------------------------------
+		function reconstruction = reconstruct(this, varargin)
+			% Compute reconstruction.
+			% Uses tomography.sinogram
+			% Initial solution (if available) should be stored in tomography.V  
+			% >> reconstruction = tomography.reconstruct(projections, iterations);
+			% >> reconstruction = tomography.reconstruct(projections, volume0, iterations);
+			
+			if ~this.initialized 
+				this.initialize();
+			end			
+			
+			if numel(varargin) == 2
+				reconstruction = this.reconstruct_c(varargin{1}, [], [], varargin{2});
+			elseif numel(varargin) == 3
+				reconstruction = this.reconstruct_c(varargin{1}, varargin{2}, [], varargin{3});
+			else
+				error('invalid parameter list')
+			end
+				
+			if strcmp(this.inner_circle,'yes')
+				reconstruction = this.selectROI(reconstruction);
+			end
+		end
+
+		%------------------------------------------------------------------
+		function reconstruction = reconstruct_mask(this, varargin)
+			% Compute reconstruction with mask.
+			% Uses tomography.sinogram  
+			% Initial solution (if available) should be stored in tomography.V  
+			% >> reconstruction = tomography.reconstructMask(projections, mask, iterations);
+			% >> reconstruction = tomography.reconstructMask(projections, volume0, mask, iterations);
+			
+			if ~this.initialized 
+				this.initialize();
+			end
+			
+			if numel(varargin) == 3
+				reconstruction = this.reconstruct_c(varargin{1}, [], varargin{2}, varargin{3});
+			elseif numel(varargin) == 4
+				reconstruction = this.reconstruct_c(varargin{1}, varargin{2}, varargin{3}, varargin{4});
+			else
+				error('invalid parameter list')
+			end
+		
+			if strcmp(this.inner_circle,'yes')
+				reconstruction = this.selectROI(reconstruction);
+			end
+
+		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);	
+					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)
+				this.cfg_base.option.ProjectionOrder = 'random';
+			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..29b963f
--- /dev/null
+++ b/matlab/algorithms/DART/IterativeTomography3D.m
@@ -0,0 +1,405 @@
+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 = IterativeTomography3D(varargin)
+			% Constructor
+			% >> tomography = IterativeTomography3D(proj_geom);
+			% >> tomography = IterativeTomography3D(proj_geom, vol_geom);
+			
+			% Input: IterativeTomography(proj_geom)
+			if nargin == 1
+				this.proj_geom = varargin{1};
+			
+			% Input: IterativeTomography(proj_geom, vol_geom)
+			elseif nargin == 2
+				this.proj_geom = varargin{1};
+				this.vol_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;
+			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 projections = project(this, volume)
+			% Compute forward projection.  
+			% >> projections = tomography.project(volume);
+			
+			if ~this.initialized 
+				this.initialize();
+			end
+
+			% project			
+			projections = this.project_c(volume);
+		end
+
+		%------------------------------------------------------------------
+		function reconstruction = reconstruct(this, varargin)
+			% Compute reconstruction.
+			% Uses tomography.sinogram
+			% Initial solution (if available) should be stored in tomography.V  
+			% >> reconstruction = tomography.reconstruct(projections, iterations);
+			% >> reconstruction = tomography.reconstruct(projections, volume0, iterations);
+			
+			if ~this.initialized 
+				this.initialize();
+			end			
+			
+			if numel(varargin) == 2
+				reconstruction = this.reconstruct_c(varargin{1}, [], [], varargin{2});
+			elseif numel(varargin) == 3
+				reconstruction = this.reconstruct_c(varargin{1}, varargin{2}, [], varargin{3});
+			else
+				error('invalid parameter list')
+			end
+				
+			if strcmp(this.inner_circle,'yes')
+				reconstruction = this.selectROI(reconstruction);
+			end
+		end
+		
+		%------------------------------------------------------------------
+		function reconstruction = reconstruct_mask(this, varargin)
+			% Compute reconstruction with mask.
+			% Uses tomography.sinogram  
+			% Initial solution (if available) should be stored in tomography.V  
+			% >> reconstruction = tomography.reconstructMask(projections, mask, iterations);
+			% >> reconstruction = tomography.reconstructMask(projections, volume0, mask, iterations);
+			
+			if ~this.initialized 
+				this.initialize();
+			end
+			
+			if numel(varargin) == 3
+				reconstruction = this.reconstruct_c(varargin{1}, [], varargin{2}, varargin{3});
+			elseif numel(varargin) == 4
+				reconstruction = this.reconstruct_c(varargin{1}, varargin{2}, varargin{3}, varargin{4});
+			else
+				error('invalid parameter list')
+			end
+		
+			if strcmp(this.inner_circle,'yes')
+				reconstruction = this.selectROI(reconstruction);
+			end
+
+		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);	
+					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/examples/example1.m b/matlab/algorithms/DART/examples/example1.m
index 9a836f8..cb02e0f 100644
--- a/matlab/algorithms/DART/examples/example1.m
+++ b/matlab/algorithms/DART/examples/example1.m
@@ -8,15 +8,11 @@
 % Website: http://sf.net/projects/astra-toolbox
 %--------------------------------------------------------------------------
 
-clear all;
-
 addpath('../');
 
-%%
-% Example 1: 2D parallel beam, cuda
-%
+%% Example 1: 2D parallel beam, cuda
 
-% Configuration
+% configuration
 proj_count		= 20;
 dart_iterations = 20;
 filename		= 'cylinders.png';
@@ -26,26 +22,20 @@ rho				= [0, 255];
 tau				= 128;
 gpu_core		= 0;
 
-% Load phantom
+% load phantom
 I = imreadgs(filename);
 
-% Create projection and volume geometries
+% create projection and volume geometries
 det_count = size(I, 1);
 angles = linspace2(0, pi, proj_count);
 proj_geom = astra_create_proj_geom('parallel', 1, det_count, angles);
-vol_geom = astra_create_vol_geom(det_count, det_count, 1);
+vol_geom = astra_create_vol_geom(det_count, det_count);
 
-% Create sinogram.
+% 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(sinogram, proj_geom);
 D.t0					= 100;
 D.t						= 10;
@@ -63,21 +53,19 @@ D.smoothing.gpu_core	= gpu_core;
 D.masking.random		= 0.1;
 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.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.statistics.proj_diff	= 'no';
 
 D.initialize();
 
 D.iterate(dart_iterations);
 
-%%
-% Convert middle slice of final iteration to png.
-%
+% save the reconstruction and the segmentation to file
 imwritesc(D.S, [outdir '/' prefix '_S.png']);
 imwritesc(D.V, [outdir '/' prefix '_V.png']);
diff --git a/matlab/algorithms/DART/examples/example2.m b/matlab/algorithms/DART/examples/example2.m
index 7fe6988..89660a5 100644
--- a/matlab/algorithms/DART/examples/example2.m
+++ b/matlab/algorithms/DART/examples/example2.m
@@ -8,15 +8,11 @@
 % Website: http://sf.net/projects/astra-toolbox
 %--------------------------------------------------------------------------
 
-clear all;
-
 addpath('../');
 
-%%
-% Example Z: 3D parallel beam, cuda
-%
+%% Example 2: 3D parallel beam, cuda
 
-% Configuration
+% configuration
 proj_count		= 20;
 dart_iterations = 20;
 outdir			= './';
@@ -25,50 +21,47 @@ rho				= [0, 0.5, 1];
 tau				= [0.25, 0.75];
 gpu_core		= 0;
 
-% Load phantom
+% load phantom
 load('phantom3d');
 
-% Create projection and volume geometries
+% create projection and volume geometries
 det_count = size(I, 1);
 slice_count = size(I,3);
 angles = linspace2(0, pi, proj_count);
 proj_geom = astra_create_proj_geom('parallel3d', 1, 1, slice_count, det_count, angles);
 vol_geom = astra_create_vol_geom(size(I));
 
-% Create sinogram
+% create sinogram
 [sinogram_id, sinogram] = astra_create_sino3d_cuda(I, proj_geom, vol_geom);
 astra_mex_data3d('delete', sinogram_id);
 
-%%
 % DART
-%
-
 D						= DARTalgorithm(sinogram, proj_geom);
 D.t0					= 100;
 D.t						= 10;
 
 D.tomography			= IterativeTomography3D();
-D.tomography.method 	= 'SIRT3D_CUDA';
-D.tomography.gpu_core 	= gpu_core;
-D.tomography.use_minc 	= 'yes';
+D.tomography.method		= 'SIRT3D_CUDA';
+D.tomography.gpu_core	= gpu_core;
+D.tomography.use_minc	= 'yes';
 
-D.segmentation.rho     = rho;
-D.segmentation.tau     = tau;
+D.segmentation.rho		= rho;
+D.segmentation.tau		= tau;
 
-D.smoothing.b          = 0.1;
-D.smoothing.full3d     = 'yes';
-D.smoothing.gpu_core   = gpu_core;
+D.smoothing.b			= 0.1;
+D.smoothing.full3d		= 'yes';
+D.smoothing.gpu_core	= gpu_core;
  
-D.masking.random       = 0.1;
-D.masking.conn         = 4;
-D.masking.gpu_core     = gpu_core;
+D.masking.random		= 0.1;
+D.masking.conn			= 4;
+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.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';
 
@@ -76,8 +69,6 @@ D.initialize();
 
 D.iterate(dart_iterations);
 
-%%
-% Convert output of final iteration to png.
-%
+% save the central slice of the reconstruction and the segmentation to file
 imwritesc(D.S(:,:,64), [outdir '/' prefix '_S_slice_64.png']);
 imwritesc(D.V(:,:,64), [outdir '/' prefix '_V_slice_64.png']);
diff --git a/matlab/algorithms/DART/examples/example3.m b/matlab/algorithms/DART/examples/example3.m
index 895630b..cc80b0f 100644
--- a/matlab/algorithms/DART/examples/example3.m
+++ b/matlab/algorithms/DART/examples/example3.m
@@ -8,15 +8,11 @@
 % Website: http://sf.net/projects/astra-toolbox
 %--------------------------------------------------------------------------
 
-clear all;
-
 addpath('../');
 
-%%
-% Example 3: 3D cone beam, cuda
-%
+%% Example 3: 3D cone beam, cuda
 
-% Configuration
+% configuration
 proj_count		= 20;
 dart_iterations = 20;
 outdir			= './';
@@ -25,24 +21,21 @@ rho				= [0, 0.5, 1];
 tau				= [0.25, 0.75];
 gpu_core		= 0;
 
-% Load phantom
+% load phantom
 load('phantom3d');
 
-% Create projection and volume geometries
+% create projection and volume geometries
 det_count = size(I, 1);
 slice_count = size(I,3);
 angles = linspace2(0, pi, proj_count);
 proj_geom = astra_create_proj_geom('cone', 1, 1, slice_count, det_count, angles, 500, 0);
 vol_geom = astra_create_vol_geom(size(I));
 
-% Create sinogram
+% create sinogram
 [sinogram_id, sinogram] = astra_create_sino3d_cuda(I, proj_geom, vol_geom);
 astra_mex_data3d('delete', sinogram_id);
 
-%%
 % DART
-%
-
 D						= DARTalgorithm(sinogram, proj_geom);
 D.t0					= 100;
 D.t						= 10;
@@ -76,8 +69,6 @@ D.initialize();
 
 D.iterate(dart_iterations);
 
-%%
-% Convert output of final iteration to png.
-%
+% save the central slice of the reconstruction and the segmentation to file
 imwritesc(D.S(:,:,64), [outdir '/' prefix '_S_slice_64.png']);
 imwritesc(D.V(:,:,64), [outdir '/' prefix '_V_slice_64.png']);
diff --git a/matlab/mex/astra_mex_algorithm_c.cpp b/matlab/mex/astra_mex_algorithm_c.cpp
index f719a6b..be1c89f 100644
--- a/matlab/mex/astra_mex_algorithm_c.cpp
+++ b/matlab/mex/astra_mex_algorithm_c.cpp
@@ -32,13 +32,14 @@ $Id$
  */
 #include <mex.h>
 #include "mexHelpFunctions.h"
+#include "astra/Globals.h"
 
 #define USE_MATLAB_UNDOCUMENTED
 
 #ifdef USE_MATLAB_UNDOCUMENTED
 extern "C" { bool utIsInterruptPending(); }
 
-#ifdef __linux__
+#ifdef USE_PTHREADS
 #define USE_PTHREADS_CTRLC
 #include <pthread.h>
 #else
@@ -49,7 +50,6 @@ extern "C" { bool utIsInterruptPending(); }
 
 
 
-#include "astra/Globals.h"
 
 #include "astra/AstraObjectManager.h"
 #include "astra/AstraObjectFactory.h"
diff --git a/matlab/mex/astra_mex_data3d_c.cpp b/matlab/mex/astra_mex_data3d_c.cpp
index db81519..35a7512 100644
--- a/matlab/mex/astra_mex_data3d_c.cpp
+++ b/matlab/mex/astra_mex_data3d_c.cpp
@@ -79,7 +79,7 @@ void astra_mex_data3d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArra
 	}
 
 	if (data && !checkDataType(data)) {
-		mexErrMsgTxt("Data must be single or double.");
+		mexErrMsgTxt("Data must be single, double or logical.");
 		return;
 	}
 
diff --git a/matlab/mex/mexCopyDataHelpFunctions.cpp b/matlab/mex/mexCopyDataHelpFunctions.cpp
index bbcad30..6dfd4a6 100644
--- a/matlab/mex/mexCopyDataHelpFunctions.cpp
+++ b/matlab/mex/mexCopyDataHelpFunctions.cpp
@@ -1,13 +1,13 @@
 /*
 -----------------------------------------------------------------------
-Copyright 2013 iMinds-Vision Lab, University of Antwerp
+Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+           2014-2015, CWI, Amsterdam
 
-Contact: astra@ua.ac.be
-Website: http://astra.ua.ac.be
+Contact: astra@uantwerpen.be
+Website: http://sf.net/projects/astra-toolbox
 
+This file is part of the ASTRA Toolbox.
 
-This file is part of the
-All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox").
 
 The ASTRA Toolbox is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -40,7 +40,6 @@ $Id$
 
 #if defined(__SSE2__)
 # include <emmintrin.h>
-
 # define STORE_32F_64F_CORE8(in, out, count, base) \
 		{\
 			const __m128 inV0 = *((const __m128 *) &in[count + 0 + base]);\
@@ -108,6 +107,17 @@ $Id$
 			out[count + 7 + base] = in[count + 7 + base];\
 		}
 #endif
+#define STORE_8F_32F_CORE8(in, out, count, base) \
+		{\
+			out[count + 0 + base] = (float)in[count + 0 + base];\
+			out[count + 1 + base] = (float)in[count + 1 + base];\
+			out[count + 2 + base] = (float)in[count + 2 + base];\
+			out[count + 3 + base] = (float)in[count + 3 + base];\
+			out[count + 4 + base] = (float)in[count + 4 + base];\
+			out[count + 5 + base] = (float)in[count + 5 + base];\
+			out[count + 6 + base] = (float)in[count + 6 + base];\
+			out[count + 7 + base] = (float)in[count + 7 + base];\
+		}
 
 const char * warnDataTypeNotSupported = "Data type not supported: nothing was copied";
 
@@ -118,6 +128,7 @@ _copyMexToCFloat32Array(const mxArray * const inArray, astra::float32 * const ou
 	const long long block = 32;
 	const long long totRoundedPixels = ROUND_DOWN(tot_size, block);
 
+	// Array of doubles
 	if (mxIsDouble(inArray)) {
 		const double * const pdMatlabData = mxGetPr(inArray);
 
@@ -133,6 +144,8 @@ _copyMexToCFloat32Array(const mxArray * const inArray, astra::float32 * const ou
 			out[count] = pdMatlabData[count];
 		}
 	}
+
+	// Array of floats
 	else if (mxIsSingle(inArray)) {
 		const float * const pfMatlabData = (const float *)mxGetData(inArray);
 
@@ -148,6 +161,23 @@ _copyMexToCFloat32Array(const mxArray * const inArray, astra::float32 * const ou
 			out[count] = pfMatlabData[count];
 		}
 	}
+
+	// Array of logicals
+	else if (mxIsLogical(inArray)) {
+		const mxLogical * const pfMatlabData = (const mxLogical *)mxGetLogicals(inArray);
+
+#pragma omp for nowait
+		for (long long count = 0; count < totRoundedPixels; count += block) {
+			STORE_8F_32F_CORE8(pfMatlabData, out, count,  0);
+			STORE_8F_32F_CORE8(pfMatlabData, out, count,  8);
+			STORE_8F_32F_CORE8(pfMatlabData, out, count, 16);
+			STORE_8F_32F_CORE8(pfMatlabData, out, count, 24);
+		}
+#pragma omp for nowait
+		for (long long count = totRoundedPixels; count < tot_size; count++) {
+			out[count] = pfMatlabData[count];
+		}
+	}
 	else {
 #pragma omp single nowait
 		mexWarnMsgIdAndTxt("ASTRA_MEX:wrong_datatype", warnDataTypeNotSupported);
diff --git a/matlab/mex/mexDataManagerHelpFunctions.cpp b/matlab/mex/mexDataManagerHelpFunctions.cpp
index 5e6020b..f9d971c 100644
--- a/matlab/mex/mexDataManagerHelpFunctions.cpp
+++ b/matlab/mex/mexDataManagerHelpFunctions.cpp
@@ -1,13 +1,13 @@
 /*
 -----------------------------------------------------------------------
-Copyright 2013 iMinds-Vision Lab, University of Antwerp
+Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+           2014-2015, CWI, Amsterdam
 
-Contact: astra@ua.ac.be
-Website: http://astra.ua.ac.be
+Contact: astra@uantwerpen.be
+Website: http://sf.net/projects/astra-toolbox
 
+This file is part of the ASTRA Toolbox.
 
-This file is part of the
-All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox").
 
 The ASTRA Toolbox is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -105,7 +105,7 @@ checkID(const astra::int32 & id, astra::CFloat32Data3DMemory *& pDataObj)
 bool
 checkDataType(const mxArray * const in)
 {
-	return (mex_is_scalar(in) || mxIsDouble(in) || mxIsSingle(in));
+	return (mex_is_scalar(in) || mxIsDouble(in) || mxIsSingle(in) || mxIsLogical(in));
 }
 
 //-----------------------------------------------------------------------------------------
diff --git a/matlab/mex/mexDataManagerHelpFunctions.h b/matlab/mex/mexDataManagerHelpFunctions.h
index 36b74d8..0614e05 100644
--- a/matlab/mex/mexDataManagerHelpFunctions.h
+++ b/matlab/mex/mexDataManagerHelpFunctions.h
@@ -1,13 +1,13 @@
 /*
 -----------------------------------------------------------------------
-Copyright 2013 iMinds-Vision Lab, University of Antwerp
+Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+           2014-2015, CWI, Amsterdam
 
-Contact: astra@ua.ac.be
-Website: http://astra.ua.ac.be
+Contact: astra@uantwerpen.be
+Website: http://sf.net/projects/astra-toolbox
 
+This file is part of the ASTRA Toolbox.
 
-This file is part of the
-All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox").
 
 The ASTRA Toolbox is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
diff --git a/matlab/mex/mexHelpFunctions.h b/matlab/mex/mexHelpFunctions.h
index 425b4ef..84372ba 100644
--- a/matlab/mex/mexHelpFunctions.h
+++ b/matlab/mex/mexHelpFunctions.h
@@ -1,13 +1,13 @@
 /*
 -----------------------------------------------------------------------
-Copyright 2012 iMinds-Vision Lab, University of Antwerp
+Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+           2014-2015, CWI, Amsterdam
 
-Contact: astra@ua.ac.be
-Website: http://astra.ua.ac.be
+Contact: astra@uantwerpen.be
+Website: http://sf.net/projects/astra-toolbox
 
+This file is part of the ASTRA Toolbox.
 
-This file is part of the
-All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox").
 
 The ASTRA Toolbox is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
diff --git a/src/AsyncAlgorithm.cpp b/src/AsyncAlgorithm.cpp
index fcc4dcb..b265f59 100644
--- a/src/AsyncAlgorithm.cpp
+++ b/src/AsyncAlgorithm.cpp
@@ -160,32 +160,6 @@ void CAsyncAlgorithm::runWrapped(int _iNrIterations)
 	m_bDone = true;
 }
 
-void CAsyncAlgorithm::timedJoin(int _milliseconds)
-{
-#ifndef USE_PTHREADS
-	if (m_pThread) {
-		boost::posix_time::milliseconds rel(_milliseconds);
-		bool res = m_pThread->timed_join(rel);
-		if (res) {
-			delete m_pThread;
-			m_pThread = 0;
-			m_bThreadStarted = false;
-		}
-	}
-#else
-	if (m_bThreadStarted) {
-		struct timespec abstime;
-		clock_gettime(CLOCK_REALTIME, &abstime);
-		abstime.tv_sec += _milliseconds / 1000;
-		abstime.tv_nsec += (_milliseconds % 1000) * 1000000L;
-		int err = pthread_timedjoin_np(m_thread, 0, &abstime);
-		if (err == 0) {
-			m_bThreadStarted = false;
-		}
-	}
-#endif
-}
-
 void CAsyncAlgorithm::signalAbort()
 {
 	if (m_pAlg)