diff options
| -rw-r--r-- | README.md | 82 | ||||
| -rw-r--r-- | astra_vc11.vcxproj | 4 | ||||
| -rw-r--r-- | build/linux/Makefile.in | 3 | ||||
| -rw-r--r-- | build/linux/acinclude.m4 | 15 | ||||
| -rwxr-xr-x | build/linux/autogen.sh | 7 | ||||
| -rw-r--r-- | build/linux/configure.ac | 61 | ||||
| -rw-r--r-- | include/astra/AsyncAlgorithm.h | 8 | ||||
| -rw-r--r-- | include/astra/Globals.h | 7 | ||||
| -rw-r--r-- | matlab/algorithms/DART/DARTalgorithm.m | 11 | ||||
| -rw-r--r-- | matlab/algorithms/DART/IterativeTomography.m | 435 | ||||
| -rw-r--r-- | matlab/algorithms/DART/IterativeTomography3D.m | 405 | ||||
| -rw-r--r-- | matlab/algorithms/DART/examples/example1.m | 40 | ||||
| -rw-r--r-- | matlab/algorithms/DART/examples/example2.m | 55 | ||||
| -rw-r--r-- | matlab/algorithms/DART/examples/example3.m | 21 | ||||
| -rw-r--r-- | matlab/mex/astra_mex_algorithm_c.cpp | 4 | ||||
| -rw-r--r-- | matlab/mex/astra_mex_data3d_c.cpp | 2 | ||||
| -rw-r--r-- | matlab/mex/mexCopyDataHelpFunctions.cpp | 42 | ||||
| -rw-r--r-- | matlab/mex/mexDataManagerHelpFunctions.cpp | 12 | ||||
| -rw-r--r-- | matlab/mex/mexDataManagerHelpFunctions.h | 10 | ||||
| -rw-r--r-- | matlab/mex/mexHelpFunctions.h | 10 | ||||
| -rw-r--r-- | src/AsyncAlgorithm.cpp | 26 |
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) |