Merge branch 'master' into volgeom3d

Conflicts:
	src/CudaBackProjectionAlgorithm3D.cpp

5 files changed, 145 insertions, 22 deletions

diff --git a/samples/matlab/s010_supersampling.m b/samples/matlab/s010_supersampling.m
index 80f6f56..148f6ad 100644
--- a/samples/matlab/s010_supersampling.m
+++ b/samples/matlab/s010_supersampling.m
@@ -12,23 +12,15 @@ vol_geom = astra_create_vol_geom(256, 256);
 proj_geom = astra_create_proj_geom('parallel', 3.0, 128, linspace2(0,pi,180));
 P = phantom(256);
 
-% Because the astra_create_sino_gpu wrapper does not have support for
-% all possible algorithm options, we manually create a sinogram
-phantom_id = astra_mex_data2d('create', '-vol', vol_geom, P);
-sinogram_id = astra_mex_data2d('create', '-sino', proj_geom);
-cfg = astra_struct('FP_CUDA');
-cfg.VolumeDataId = phantom_id;
-cfg.ProjectionDataId = sinogram_id;
+% We create a projector set up to use 3 rays per detector element
+cfg_proj = astra_struct('cuda');
+cfg_proj.option.DetectorSuperSampling = 3;
+cfg_proj.ProjectionGeometry = proj_geom;
+cfg_proj.VolumeGeometry = vol_geom;
+proj_id = astra_mex_projector('create', cfg_proj);
 
-% Set up 3 rays per detector element
-cfg.option.DetectorSuperSampling = 3;
 
-alg_id = astra_mex_algorithm('create', cfg);
-astra_mex_algorithm('run', alg_id);
-astra_mex_algorithm('delete', alg_id);
-astra_mex_data2d('delete', phantom_id);
-
-sinogram3 = astra_mex_data2d('get', sinogram_id);
+[sinogram3 sinogram_id] = astra_create_sino(P, proj_id);
 
 figure(1); imshow(P, []);
 figure(2); imshow(sinogram3, []);
@@ -39,14 +31,14 @@ rec_id = astra_mex_data2d('create', '-vol', vol_geom);
 cfg = astra_struct('SIRT_CUDA');
 cfg.ReconstructionDataId = rec_id;
 cfg.ProjectionDataId = sinogram_id;
-% Set up 3 rays per detector element
-cfg.option.DetectorSuperSampling = 3;
+cfg.ProjectorId = proj_id;
+
 
 % There is also an option for supersampling during the backprojection step.
 % This should be used if your detector pixels are smaller than the voxels.
 
 % Set up 2 rays per image pixel dimension, for 4 rays total per image pixel.
-% cfg.option.PixelSuperSampling = 2;
+% cfg_proj.option.PixelSuperSampling = 2;
 
 
 alg_id = astra_mex_algorithm('create', cfg);
diff --git a/samples/matlab/s017_opTomo.m b/samples/matlab/s017_opTomo.m
new file mode 100644
index 0000000..891a93d
--- /dev/null
+++ b/samples/matlab/s017_opTomo.m
@@ -0,0 +1,62 @@
+% -----------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+% 
+% Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+%            2014-2015, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Contact: astra@uantwerpen.be
+% Website: http://sf.net/projects/astra-toolbox
+% -----------------------------------------------------------------------
+
+% This sample illustrates the use of opTomo.
+%
+% opTomo is a wrapper around the FP and BP operations of the ASTRA Toolbox,
+% to allow you to use them as you would a matrix.
+%
+% This class requires the Spot Linear-Operator Toolbox to be installed.
+% You can download this at http://www.cs.ubc.ca/labs/scl/spot/
+
+% load a phantom image
+im = phantom(256);
+% and flatten it to a vector
+x = im(:);
+
+%% Setting up the geometry
+% projection geometry
+proj_geom = astra_create_proj_geom('parallel', 1, 256, linspace2(0,pi,180));
+% object dimensions
+vol_geom  = astra_create_vol_geom(256,256);
+
+%% Generate projection data
+% Create the Spot operator for ASTRA using the GPU.
+W = opTomo('cuda', proj_geom, vol_geom);
+
+p = W*x;
+
+% reshape the vector into a sinogram
+sinogram = reshape(p, W.proj_size);
+imshow(sinogram, []);
+
+
+%% Reconstruction
+% We use a least squares solver lsqr from Matlab to solve the 
+% equation W*x = p.
+% Max number of iterations is 100, convergence tolerance of 1e-6.
+y = lsqr(W, p, 1e-6, 100);
+
+% the output is a vector, so we reshape it into an image
+reconstruction = reshape(y, W.vol_size);
+
+subplot(1,3,1);
+imshow(reconstruction, []);
+title('Reconstruction');
+
+subplot(1,3,2);
+imshow(im, []);
+title('Ground truth');
+
+% The transpose of the operator corresponds to the backprojection.
+backProjection = W'*p;
+subplot(1,3,3);
+imshow(reshape(backProjection, W.vol_size), []);
+title('Backprojection');
diff --git a/samples/python/s005_3d_geometry.py b/samples/python/s005_3d_geometry.py
index f43fc7e..a7f7a3d 100644
--- a/samples/python/s005_3d_geometry.py
+++ b/samples/python/s005_3d_geometry.py
@@ -24,7 +24,11 @@
 #
 #-----------------------------------------------------------------------
 
-from six.moves import range
+try:
+    from six.moves import range
+except ImportError:
+    # six 1.3.0
+    from six.moves import xrange as range
 import astra
 import numpy as np
 
diff --git a/samples/python/s016_plots.py b/samples/python/s016_plots.py
index cd4d98c..8a8ba64 100644
--- a/samples/python/s016_plots.py
+++ b/samples/python/s016_plots.py
@@ -24,7 +24,11 @@
 #
 #-----------------------------------------------------------------------
 
-from six.moves import range
+try:
+    from six.moves import range
+except ImportError:
+    # six 1.3.0
+    from six.moves import xrange as range
 import astra
 import numpy as np
 
@@ -35,8 +39,8 @@ proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180
 # As before, create a sinogram from a phantom
 import scipy.io
 P = scipy.io.loadmat('phantom.mat')['phantom256']
-proj_id = astra.create_projector('line',proj_geom,vol_geom)
-sinogram_id, sinogram = astra.create_sino(P, proj_id,useCUDA=True)
+proj_id = astra.create_projector('cuda',proj_geom,vol_geom)
+sinogram_id, sinogram = astra.create_sino(P, proj_id)
 
 import pylab
 pylab.gray()
diff --git a/samples/python/s017_OpTomo.py b/samples/python/s017_OpTomo.py
new file mode 100644
index 0000000..967fa64
--- /dev/null
+++ b/samples/python/s017_OpTomo.py
@@ -0,0 +1,61 @@
+#-----------------------------------------------------------------------
+#Copyright 2013 Centrum Wiskunde & Informatica, Amsterdam
+#
+#Author: Daniel M. Pelt
+#Contact: D.M.Pelt@cwi.nl
+#Website: http://dmpelt.github.io/pyastratoolbox/
+#
+#
+#This file is part of the Python interface to the
+#All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox").
+#
+#The Python interface to 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
+#the Free Software Foundation, either version 3 of the License, or
+#(at your option) any later version.
+#
+#The Python interface to the ASTRA Toolbox is distributed in the hope that it will be useful,
+#but WITHOUT ANY WARRANTY; without even the implied warranty of
+#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#GNU General Public License for more details.
+#
+#You should have received a copy of the GNU General Public License
+#along with the Python interface to the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+#
+#-----------------------------------------------------------------------
+
+import astra
+import numpy as np
+import scipy.sparse.linalg
+
+vol_geom = astra.create_vol_geom(256, 256)
+proj_geom = astra.create_proj_geom('parallel', 1.0, 384, np.linspace(0,np.pi,180,False))
+
+# As before, create a sinogram from a phantom
+import scipy.io
+P = scipy.io.loadmat('phantom.mat')['phantom256']
+proj_id = astra.create_projector('cuda',proj_geom,vol_geom)
+
+# construct the OpTomo object
+W = astra.OpTomo(proj_id)
+
+sinogram = W * P
+sinogram = sinogram.reshape([180, 384])
+
+import pylab
+pylab.gray()
+pylab.figure(1)
+pylab.imshow(P)
+pylab.figure(2)
+pylab.imshow(sinogram)
+
+# Run the lsqr linear solver
+output = scipy.sparse.linalg.lsqr(W, sinogram.flatten(), iter_lim=150)
+rec = output[0].reshape([256, 256])
+
+pylab.figure(3)
+pylab.imshow(rec)
+pylab.show()
+
+# Clean up.
+astra.projector.delete(proj_id)