Add sample script using parallel_vec

2 files changed, 167 insertions, 0 deletions

diff --git a/samples/matlab/s022_fbp_cor.m b/samples/matlab/s022_fbp_cor.m
new file mode 100644
index 0000000..b33d5ce
--- /dev/null
+++ b/samples/matlab/s022_fbp_cor.m
@@ -0,0 +1,69 @@
+% -----------------------------------------------------------------------
+% This file is part of the ASTRA Toolbox
+% 
+% Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp
+%            2014-2016, CWI, Amsterdam
+% License: Open Source under GPLv3
+% Contact: astra@uantwerpen.be
+% Website: http://www.astra-toolbox.com/
+% -----------------------------------------------------------------------
+
+cor_shift = 3.6;
+
+vol_geom = astra_create_vol_geom(256, 256);
+proj_geom = astra_create_proj_geom('parallel', 1.0, 256, linspace2(0,pi,180));
+
+% Projection geometry with shifted center of rotation
+% We create this by shifting the detector center V(:,3:4) by a multiple of the detector
+% orientation V(:,5:6).
+proj_geom_cor = astra_geom_2vec(proj_geom);
+proj_geom_cor.Vectors(:,3:4) = proj_geom_cor.Vectors(:,3:4) + cor_shift * proj_geom_cor.Vectors(:,5:6);
+
+% As before, create a sinogram from a phantom, using the shifted center of rotation
+P = phantom(256);
+[sinogram_id, sinogram] = astra_create_sino_gpu(P, proj_geom_cor, vol_geom);
+figure(1); imshow(P, []);
+figure(2); imshow(sinogram, []);
+
+astra_mex_data2d('delete', sinogram_id);
+
+% We now re-create the sinogram data object as we would do when loading
+% an external sinogram, using a standard geometry, and try to do a reconstruction,
+% to show the misalignment artifacts caused by the shifted center of rotation
+sinogram_id = astra_mex_data2d('create', '-sino', proj_geom, sinogram);
+
+% Create a data object for the reconstruction
+rec_id = astra_mex_data2d('create', '-vol', vol_geom);
+
+% Set up the parameters for a reconstruction algorithm using the GPU
+cfg = astra_struct('FBP_CUDA');
+cfg.ReconstructionDataId = rec_id;
+cfg.ProjectionDataId = sinogram_id;
+alg_id = astra_mex_algorithm('create', cfg);
+astra_mex_algorithm('run', alg_id);
+
+% Get the result
+rec = astra_mex_data2d('get', rec_id);
+figure(3); imshow(rec, []);
+
+astra_mex_algorithm('delete', alg_id);
+astra_mex_data2d('delete', rec_id);
+
+% Now change back to the proper, shifted geometry, and do another reconstruction
+astra_mex_data2d('change_geometry', sinogram_id, proj_geom_cor);
+rec_id = astra_mex_data2d('create', '-vol', vol_geom);
+
+cfg = astra_struct('FBP_CUDA');
+cfg.ReconstructionDataId = rec_id;
+cfg.ProjectionDataId = sinogram_id;
+alg_id = astra_mex_algorithm('create', cfg);
+astra_mex_algorithm('run', alg_id);
+
+% Get the result
+rec = astra_mex_data2d('get', rec_id);
+figure(4); imshow(rec, []);
+
+
+astra_mex_algorithm('delete', alg_id);
+astra_mex_data2d('delete', rec_id);
+astra_mex_data2d('delete', sinogram_id);
diff --git a/samples/python/s022_fbp_cor.py b/samples/python/s022_fbp_cor.py
new file mode 100644
index 0000000..eb2f56b
--- /dev/null
+++ b/samples/python/s022_fbp_cor.py
@@ -0,0 +1,98 @@
+# -----------------------------------------------------------------------
+# Copyright: 2010-2016, iMinds-Vision Lab, University of Antwerp
+#            2013-2016, CWI, Amsterdam
+#
+# Contact: astra@uantwerpen.be
+# Website: http://www.astra-toolbox.com/
+#
+# This file is part of the 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
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# 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 ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
+#
+# -----------------------------------------------------------------------
+
+import astra
+import numpy as np
+
+cor_shift = 3.6
+
+vol_geom = astra.create_vol_geom(256, 256)
+proj_geom = astra.create_proj_geom('parallel', 1.0, 256, np.linspace(0,np.pi,180,False))
+
+# Projection geometry with shifted center of rotation
+# We create this by shifting the detector center V[:,2:4] by a multiple of the detector
+# orientation V[:,4:6].
+proj_geom_cor = astra.functions.geom_2vec(proj_geom)
+V = proj_geom_cor['Vectors']
+V[:,2:4] = V[:,2:4] + cor_shift * V[:,4:6]
+
+# As before, create a sinogram from a phantom, using the shifted center of rotation
+import scipy.io
+P = scipy.io.loadmat('phantom.mat')['phantom256']
+
+proj_id_cor = astra.create_projector('cuda',proj_geom_cor,vol_geom)
+sinogram_id, sinogram = astra.create_sino(P, proj_id_cor)
+
+# Change the projection geometry metadata attached to the sinogram to standard geometry,
+# and try to do a reconstruction, to show the misalignment artifacts caused by
+# the shifted center of rotation
+astra.data2d.change_geometry(sinogram_id, proj_geom)
+
+import pylab
+pylab.gray()
+pylab.figure(1)
+pylab.imshow(P)
+pylab.figure(2)
+pylab.imshow(sinogram)
+
+# Create a data object for the reconstruction
+rec_id = astra.data2d.create('-vol', vol_geom)
+
+# Set up the parameters for a reconstruction algorithm using the GPU
+cfg = astra.astra_dict('FBP_CUDA')
+cfg['ReconstructionDataId'] = rec_id
+cfg['ProjectionDataId'] = sinogram_id
+alg_id = astra.algorithm.create(cfg)
+astra.algorithm.run(alg_id)
+
+# Get the result
+rec = astra.data2d.get(rec_id)
+pylab.figure(3)
+pylab.imshow(rec)
+
+astra.algorithm.delete(alg_id)
+astra.data2d.delete(rec_id)
+
+# Now change back to the proper, shifted geometry, and do another reconstruction
+astra.data2d.change_geometry(sinogram_id, proj_geom_cor)
+rec_id = astra.data2d.create('-vol', vol_geom)
+cfg = astra.astra_dict('FBP_CUDA')
+cfg['ReconstructionDataId'] = rec_id
+cfg['ProjectionDataId'] = sinogram_id
+alg_id = astra.algorithm.create(cfg)
+astra.algorithm.run(alg_id)
+
+# Get the result
+rec = astra.data2d.get(rec_id)
+pylab.figure(4)
+pylab.imshow(rec)
+pylab.show()
+
+
+
+astra.data2d.delete(sinogram_id)
+astra.projector.delete(proj_id_cor)
+astra.algorithm.delete(alg_id)
+astra.data2d.delete(rec_id)