Merge branch 'parallel_vec'

8 files changed, 379 insertions, 55 deletions

diff --git a/matlab/algorithms/DART/IterativeTomography.m b/matlab/algorithms/DART/IterativeTomography.m
index e94dd36..df414a0 100644
--- a/matlab/algorithms/DART/IterativeTomography.m
+++ b/matlab/algorithms/DART/IterativeTomography.m
@@ -87,7 +87,7 @@ classdef IterativeTomography < matlab.mixin.Copyable
 		function ok = initialize(this)
 			% Initialize this object.  Returns 1 if succesful.
 			% >> tomography.initialize();
-			
+			disp('sdfqnlmkqdsfmlkjdfqsjklm');
 			% create projection geometry with super-resolution
 			if this.superresolution > 1
 				this.proj_geom_sr = astra_geom_superresolution(this.proj_geom, this.superresolution);
diff --git a/matlab/algorithms/DART/examples/example1.m b/matlab/algorithms/DART/examples/example1.m
index 6c86473..25d5215 100644
--- a/matlab/algorithms/DART/examples/example1.m
+++ b/matlab/algorithms/DART/examples/example1.m
@@ -35,36 +35,41 @@ vol_geom = astra_create_vol_geom(det_count, det_count);
 [sinogram_id, sinogram] = astra_create_sino_cuda(I, proj_geom, vol_geom);
 astra_mex_data2d('delete', sinogram_id);
 
-% DART
-D						= DARTalgorithm(sinogram, proj_geom);
-D.t0					= 100;
-D.t						= 10;
+	% DART
+	D						= DARTalgorithm(sinogram, proj_geom);
+	D.t0					= 100;
+	D.t						= 10;
 
-D.tomography.method		= 'SIRT_CUDA';
-D.tomography.gpu_core	= gpu_core;
-D.tomography.use_minc	= 'yes';
+	D.tomography.method		= 'SIRT';
+	D.tomography.gpu_core	= gpu_core;
+	D.tomography.use_minc	= 'yes';
+	D.tomography.gpu        = 'no';
 
-D.segmentation.rho		= rho;
-D.segmentation.tau		= tau;
+	D.segmentation          = SegmentationPDM();
+	D.segmentation.rho		= rho*1.8;
+	D.segmentation.tau		= tau*1.5;
+	D.segmentation.interval = 5;
 
-D.smoothing.b			= 0.1;
-D.smoothing.gpu_core	= gpu_core;
- 
-D.masking.random		= 0.1;
-D.masking.gpu_core		= gpu_core;
+	D.smoothing.b			= 0.1;
+	D.smoothing.gpu_core	= gpu_core;
+	D.smoothing.gpu        = 'no';
+	
+	D.masking.random		= 0.1;
+	D.masking.gpu_core		= gpu_core;
+	D.masking.gpu        = 'no';
+	
+	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.initialize();
-
-D.iterate(dart_iterations);
+	D.iterate(dart_iterations);
 
 % save the reconstruction and the segmentation to file
 imwritesc(D.S, [outdir '/' prefix '_S.png']);
diff --git a/matlab/mex/astra_mex_data2d_c.cpp b/matlab/mex/astra_mex_data2d_c.cpp
index ee558e0..3d3a253 100644
--- a/matlab/mex/astra_mex_data2d_c.cpp
+++ b/matlab/mex/astra_mex_data2d_c.cpp
@@ -44,6 +44,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #include "astra/SparseMatrixProjectionGeometry2D.h"
 #include "astra/FanFlatProjectionGeometry2D.h"
 #include "astra/FanFlatVecProjectionGeometry2D.h"
+#include "astra/ParallelVecProjectionGeometry2D.h"
 
 using namespace std;
 using namespace astra;
@@ -159,6 +160,8 @@ void astra_mex_data2d_create(int& nlhs, mxArray* plhs[], int& nrhs, const mxArra
 			pGeometry = new CFanFlatProjectionGeometry2D();	
 		} else if (type == "fanflat_vec") {
 			pGeometry = new CFanFlatVecProjectionGeometry2D();	
+		} else if (type == "parallel_vec") {
+			pGeometry = new CParallelVecProjectionGeometry2D();	
 		} else {
 			pGeometry = new CParallelProjectionGeometry2D();	
 		}
@@ -448,6 +451,8 @@ void astra_mex_data2d_change_geometry(int nlhs, mxArray* plhs[], int nrhs, const
 			pGeometry = new CFanFlatProjectionGeometry2D();	
 		} else if (type == "fanflat_vec") {
 			pGeometry = new CFanFlatVecProjectionGeometry2D();	
+		} else if (type == "parallel_vec") {
+			pGeometry = new CParallelVecProjectionGeometry2D();	
 		} else {
 			pGeometry = new CParallelProjectionGeometry2D();	
 		}
diff --git a/matlab/tools/astra_create_proj_geom.m b/matlab/tools/astra_create_proj_geom.m
index 0a30f23..01c9177 100644
--- a/matlab/tools/astra_create_proj_geom.m
+++ b/matlab/tools/astra_create_proj_geom.m
@@ -106,6 +106,19 @@ if strcmp(type,'parallel')
 		'ProjectionAngles',		varargin{3}  ...
 	);
 
+elseif strcmp(type,'parallel_vec')
+	if numel(varargin) < 2
+		error('not enough variables: astra_create_proj_geom(parallel_vec, det_count, V')
+	end
+	if size(varargin{2}, 2) ~= 6
+		error('V should be a Nx6 matrix, with N the number of projections')
+	end
+	proj_geom = struct( ...
+		'type',					'parallel_vec',  ...
+		'DetectorCount',		varargin{1}, ...
+		'Vectors',				varargin{2}  ...
+	);
+
 elseif strcmp(type,'fanflat')
 	if numel(varargin) < 5
 		error('not enough variables: astra_create_proj_geom(fanflat, det_width, det_count, angles, source_origin, source_det)');
diff --git a/matlab/tools/astra_geom_2vec.m b/matlab/tools/astra_geom_2vec.m
index 0abd07c..e563f47 100644
--- a/matlab/tools/astra_geom_2vec.m
+++ b/matlab/tools/astra_geom_2vec.m
@@ -1,25 +1,65 @@
-function proj_geom_out = astra_geom_2vec(proj_geom)
+function proj_geom_vec = astra_geom_2vec(proj_geom)
+%--------------------------------------------------------------------------
+% proj_geom_vec = astra_geom_2vec(proj_geom)
+%
+% Convert a conventional projection geometry to a corresponding vector-base
+% projection geometry
+%
+% proj_geom: input projection geometry (parallel, fanflat, parallel3d, cone)
+% proj_geom_vec: output vector-base projection geometry
+%--------------------------------------------------------------------------
+% 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
+%--------------------------------------------------------------------------
+% $Id$
+
+	% PARALLEL 2D
+	if strcmp(proj_geom.type,'parallel')
+
+		vectors = zeros(numel(proj_geom.ProjectionAngles), 6);
+		for i = 1:numel(proj_geom.ProjectionAngles)
+
+			% ray direction
+			vectors(i,1) = sin(proj_geom.ProjectionAngles(i));
+			vectors(i,2) = -cos(proj_geom.ProjectionAngles(i));
+
+			% center of detector
+			vectors(i,3) = 0;
+			vectors(i,4) = 0;
+
+			% vector from detector pixel 0 to 1
+			vectors(i,5) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DetectorWidth;
+			vectors(i,6) = sin(proj_geom.ProjectionAngles(i)) * proj_geom.DetectorWidth;
+
+		end
+
+		proj_geom_vec = astra_create_proj_geom('parallel_vec', proj_geom.DetectorCount, vectors);
 
 	% FANFLAT
-	if strcmp(proj_geom.type,'fanflat')
+	elseif strcmp(proj_geom.type,'fanflat')
 
 		vectors = zeros(numel(proj_geom.ProjectionAngles), 6);
 		for i = 1:numel(proj_geom.ProjectionAngles)
 
 			% source
 			vectors(i,1) = sin(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginSource;
-			vectors(i,2) = -cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginSource;	
+			vectors(i,2) = -cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginSource;
 
 			% center of detector
 			vectors(i,3) = -sin(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginDetector;
-			vectors(i,4) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginDetector;	
+			vectors(i,4) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginDetector;
 
 			% vector from detector pixel 0 to 1
 			vectors(i,5) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DetectorWidth;
 			vectors(i,6) = sin(proj_geom.ProjectionAngles(i)) * proj_geom.DetectorWidth;
 		end
 
-		proj_geom_out = astra_create_proj_geom('fanflat_vec', proj_geom.DetectorCount, vectors);
+		proj_geom_vec = astra_create_proj_geom('fanflat_vec', proj_geom.DetectorCount, vectors);
 
 	% CONE
 	elseif strcmp(proj_geom.type,'cone')
@@ -29,13 +69,13 @@ function proj_geom_out = astra_geom_2vec(proj_geom)
 
 			% source
 			vectors(i,1) = sin(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginSource;
-			vectors(i,2) = -cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginSource;	
-			vectors(i,3) = 0;	
+			vectors(i,2) = -cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginSource;
+			vectors(i,3) = 0;
 
 			% center of detector
 			vectors(i,4) = -sin(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginDetector;
-			vectors(i,5) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginDetector;	
-			vectors(i,6) = 0;	
+			vectors(i,5) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DistanceOriginDetector;
+			vectors(i,6) = 0;
 
 			% vector from detector pixel (0,0) to (0,1)
 			vectors(i,7) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DetectorSpacingX;
@@ -45,13 +85,13 @@ function proj_geom_out = astra_geom_2vec(proj_geom)
 			% vector from detector pixel (0,0) to (1,0)
 			vectors(i,10) = 0;
 			vectors(i,11) = 0;
-			vectors(i,12) = proj_geom.DetectorSpacingY;		
+			vectors(i,12) = proj_geom.DetectorSpacingY;
 		end
 
-		proj_geom_out = astra_create_proj_geom('cone_vec', proj_geom.DetectorRowCount, proj_geom.DetectorColCount, vectors);	
+		proj_geom_vec = astra_create_proj_geom('cone_vec', proj_geom.DetectorRowCount, proj_geom.DetectorColCount, vectors);
 
-	% PARALLEL
-	elseif strcmp(proj_geom.type,'parallel3d')	
+	% PARALLEL 3D
+	elseif strcmp(proj_geom.type,'parallel3d')
 
 		vectors = zeros(numel(proj_geom.ProjectionAngles), 12);
 		for i = 1:numel(proj_geom.ProjectionAngles)
@@ -59,12 +99,12 @@ function proj_geom_out = astra_geom_2vec(proj_geom)
 			% ray direction
 			vectors(i,1) = sin(proj_geom.ProjectionAngles(i));
 			vectors(i,2) = -cos(proj_geom.ProjectionAngles(i));
-			vectors(i,3) = 0;	
+			vectors(i,3) = 0;
 
 			% center of detector
 			vectors(i,4) = 0;
 			vectors(i,5) = 0;
-			vectors(i,6) = 0;	
+			vectors(i,6) = 0;
 
 			% vector from detector pixel (0,0) to (0,1)
 			vectors(i,7) = cos(proj_geom.ProjectionAngles(i)) * proj_geom.DetectorSpacingX;
@@ -74,11 +114,13 @@ function proj_geom_out = astra_geom_2vec(proj_geom)
 			% vector from detector pixel (0,0) to (1,0)
 			vectors(i,10) = 0;
 			vectors(i,11) = 0;
-			vectors(i,12) = proj_geom.DetectorSpacingY;		
+			vectors(i,12) = proj_geom.DetectorSpacingY;
 		end
 
-		proj_geom_out = astra_create_proj_geom('parallel3d_vec', proj_geom.DetectorRowCount, proj_geom.DetectorColCount, vectors);		
+		proj_geom_vec = astra_create_proj_geom('parallel3d_vec', proj_geom.DetectorRowCount, proj_geom.DetectorColCount, vectors);
 
 	else
 		error(['No suitable vector geometry found for type: ' proj_geom.type])
 	end
+
+end
diff --git a/matlab/tools/astra_geom_postalignment.m b/matlab/tools/astra_geom_postalignment.m
index 4115af2..c16f8ed 100644
--- a/matlab/tools/astra_geom_postalignment.m
+++ b/matlab/tools/astra_geom_postalignment.m
@@ -1,11 +1,36 @@
 function proj_geom = astra_geom_postalignment(proj_geom, factor)
+%--------------------------------------------------------------------------
+% proj_geom = astra_geom_postalignment(proj_geom, factorU)
+% proj_geom = astra_geom_postalignment(proj_geom, [factorU factorV])
+%
+% Apply a postalignment to a vector-based projection geometry.  Can be used to model the rotation axis offset.
+%
+% proj_geom: input projection geometry (vector-based only, use astra_geom_2vec to convert conventional projection geometries)
+% dim (optional): which dimension
+% s: output
+%--------------------------------------------------------------------------
+% 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
+%--------------------------------------------------------------------------
+% $Id$
+
+	if strcmp(proj_geom.type,'fanflat_vec') || strcmp(proj_geom.type,'parallel_vec')
+		proj_geom.Vectors(:,3:4) = proj_geom.Vectors(:,3:4) + factor(1) * proj_geom.Vectors(:,5:6);
 
-	if strcmp(proj_geom.type,'fanflat_vec')
-		proj_geom.Vectors(:,3:4) = proj_geom.Vectors(:,3:4) + factor * proj_geom.Vectors(:,5:6);
-		
 	elseif strcmp(proj_geom.type,'cone_vec') || strcmp(proj_geom.type,'parallel3d_vec')
-		proj_geom.Vectors(:,4:6) = proj_geom.Vectors(:,4:6) + factor * proj_geom.Vectors(:,7:9);
+		if numel(factor) == 1
+			proj_geom.Vectors(:,4:6) = proj_geom.Vectors(:,4:6) + factor * proj_geom.Vectors(:,7:9);
+		elseif numel(factor) > 1
+			proj_geom.Vectors(:,4:6) = proj_geom.Vectors(:,4:6) + factor(1) * proj_geom.Vectors(:,7:9) + factor(2) * proj_geom.Vectors(:,10:12);
+		end
 
 	else
 		error('Projection geometry not suited for postalignment correction.')
 	end
+
+end
\ No newline at end of file
diff --git a/matlab/tools/astra_geom_size.m b/matlab/tools/astra_geom_size.m
index 7044515..b3c1522 100644
--- a/matlab/tools/astra_geom_size.m
+++ b/matlab/tools/astra_geom_size.m
@@ -1,4 +1,22 @@
 function s = astra_geom_size(geom, dim)
+%--------------------------------------------------------------------------
+% s = astra_geom_size(geom, dim)
+%
+% Get the size of a volume or projection geometry.
+%
+% geom: volume or projection geometry
+% dim (optional): which dimension
+% s: output
+%--------------------------------------------------------------------------
+% 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
+%--------------------------------------------------------------------------
+% $Id$
 
 	if isfield(geom, 'GridSliceCount')
 		% 3D Volume geometry?
@@ -6,23 +24,23 @@ function s = astra_geom_size(geom, dim)
 	elseif isfield(geom, 'GridColCount')
 		% 2D Volume geometry?
 		s = [ geom.GridRowCount, geom.GridColCount ];
-	elseif strcmp(geom.type,'parallel') || strcmp(geom.type,'fanflat') 
+	elseif strcmp(geom.type,'parallel') || strcmp(geom.type,'fanflat')
 		s = [numel(geom.ProjectionAngles), geom.DetectorCount];
-		
-	elseif strcmp(geom.type,'parallel3d') || strcmp(geom.type,'cone') 
+
+	elseif strcmp(geom.type,'parallel3d') || strcmp(geom.type,'cone')
 		s =  [geom.DetectorColCount, numel(geom.ProjectionAngles), geom.DetectorRowCount];
-		
-	elseif strcmp(geom.type,'fanflat_vec')
+
+	elseif strcmp(geom.type,'fanflat_vec') || strcmp(geom.type,'parallel_vec')
 		s = [size(geom.Vectors,1), geom.DetectorCount];
-		
-	elseif strcmp(geom.type,'parallel3d_vec') || strcmp(geom.type,'cone_vec') 		
+
+	elseif strcmp(geom.type,'parallel3d_vec') || strcmp(geom.type,'cone_vec')
 		s = [geom.DetectorColCount, size(geom.Vectors,1), geom.DetectorRowCount];
-		
+
 	end
 
 	if nargin == 2
 		s = s(dim);
 	end
-	
+
 end
 
diff --git a/matlab/tools/astra_geom_visualize.m b/matlab/tools/astra_geom_visualize.m
new file mode 100644
index 0000000..0044844
--- /dev/null
+++ b/matlab/tools/astra_geom_visualize.m
@@ -0,0 +1,216 @@
+function astra_geom_visualize(proj_geom, vol_geom)
+
+	if strcmp(proj_geom.type,'parallel') ||  strcmp(proj_geom.type,'fanflat') ||strcmp(proj_geom.type,'parallel3d') ||  strcmp(proj_geom.type,'cone')
+		proj_geom = astra_geom_2vec(proj_geom);
+	end
+
+	% open window
+	f = figure('Visible','off');
+	hold on
+
+	% display projection 1
+	displayProjection(1);
+
+	% label
+	txt = uicontrol('Style','text', 'Position',[10 10 70 20], 'String','Projection');
+
+	% slider
+	anglecount = size(proj_geom.Vectors,1);
+	sld = uicontrol('Style', 'slider', ...
+		'Min', 1, 'Max', anglecount, 'SliderStep', [1 1]/anglecount, 'Value', 1, ...
+		'Position', [80 10 200 20], ...
+		'Callback', @updateProjection); 
+
+	f.Visible = 'on';
+
+    function updateProjection(source, callbackdata)
+		displayProjection(floor(source.Value));
+    end
+
+    function displayProjection(a)
+
+		colours = get(gca,'ColorOrder');
+
+    	
+    	% set title
+    	title(['projection ' num2str(a)]);
+
+		if strcmp(proj_geom.type,'parallel_vec')
+
+			v = proj_geom.Vectors;
+			d = proj_geom.DetectorCount;
+
+			if ~isfield(vol_geom, 'option')
+				minx = -vol_geom.GridRowCount/2;
+				miny = -vol_geom.GridColCount/2;
+				minz = -vol_geom.GridSliceCount/2;
+				maxx = vol_geom.GridRowCount/2;
+			else
+				minx = vol_geom.option.WindowMinX;
+				miny = vol_geom.option.WindowMinY;
+				maxx = vol_geom.option.WindowMaxX;
+				maxy = vol_geom.option.WindowMaxY;
+			end
+		
+			% axis
+			cla
+			axis([minx maxx miny maxy]*2.25)
+			axis square
+
+			% volume
+			plot([minx minx maxx maxx minx], [miny maxy maxy miny miny], 'LineWidth', 1, 'Color', colours(1,:))
+
+			% ray
+			s = maxx - minx;
+			plot([0 v(a,1)]*s*0.33, [0 v(a,2)]*s*0.33, 'LineWidth', 2, 'Color', colours(3,:))
+
+			% detector
+			s2 = s*0.75;
+			plot([-d/2 d/2]*v(a,5) + v(a,3) + s2*v(a,1), [-d/2 d/2]*v(a,6) + v(a,4) + s2*v(a,2), 'LineWidth', 2, 'Color', colours(5,:))
+
+		elseif strcmp(proj_geom.type,'fanflat_vec')
+
+			v = proj_geom.Vectors;
+			d = proj_geom.DetectorCount;
+
+			if ~isfield(vol_geom, 'option')
+				minx = -vol_geom.GridRowCount/2;
+				miny = -vol_geom.GridColCount/2;
+				minz = -vol_geom.GridSliceCount/2;
+				maxx = vol_geom.GridRowCount/2;
+			else
+				minx = vol_geom.option.WindowMinX;
+				miny = vol_geom.option.WindowMinY;
+				maxx = vol_geom.option.WindowMaxX;
+				maxy = vol_geom.option.WindowMaxY;
+			end
+
+			% axis
+			cla
+			axis([minx maxx miny maxy]*2.25)
+			axis square
+
+			% volume
+			plot([minx minx maxx maxx minx], [miny maxy maxy miny miny], 'LineWidth', 1, 'Color', colours(1,:))
+
+			% detector
+			D1 = v(a,3:4) - d/2*v(a,5:6);
+			D2 = v(a,3:4) + d/2*v(a,5:6);
+			plot([D1(1) D2(1)], [D1(2) D2(2)], 'LineWidth', 2, 'Color', colours(5,:))
+
+			% beam
+			plot([v(a,1) D1(1)], [v(a,2) D1(2)], 'LineWidth', 1, 'Color', colours(3,:))
+			plot([v(a,1) D2(1)], [v(a,2) D2(2)], 'LineWidth', 1, 'Color', colours(3,:))
+
+		elseif strcmp(proj_geom.type,'parallel3d_vec')
+
+			v = proj_geom.Vectors;
+			d1 = proj_geom.DetectorRowCount;
+			d2 = proj_geom.DetectorColCount;
+
+			if ~isfield(vol_geom, 'option')
+				minx = -vol_geom.GridRowCount/2;
+				miny = -vol_geom.GridColCount/2;
+				minz = -vol_geom.GridSliceCount/2;
+				maxx = vol_geom.GridRowCount/2;
+				maxy = vol_geom.GridColCount/2;
+				maxz = vol_geom.GridSliceCount/2;
+			else
+				minx = vol_geom.option.WindowMinX;
+				miny = vol_geom.option.WindowMinY;
+				minz = vol_geom.option.WindowMinZ;
+				maxx = vol_geom.option.WindowMaxX;
+				maxy = vol_geom.option.WindowMaxY;
+				maxz = vol_geom.option.WindowMaxZ;
+			end
+
+			% axis
+			windowminx = min(v(:,4));
+			windowminy = min(v(:,5));
+			windowminz = max(v(:,6));
+			windowmaxx = max(v(:,4));
+			windowmaxy = max(v(:,5));
+			windowmaxz = max(v(:,6));
+			cla
+			axis([minx maxx miny maxy minz maxz]*5.10)
+
+			% volume
+			plot3([minx minx maxx maxx minx], [miny maxy maxy miny miny], [minz minz minz minz minz], 'LineWidth', 1, 'Color', colours(1,:))
+			plot3([minx minx maxx maxx minx], [miny maxy maxy miny miny], [maxz maxz maxz maxz maxz], 'LineWidth', 1, 'Color', colours(1,:))
+			plot3([minx minx], [miny miny], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:))
+			plot3([maxx maxx], [miny miny], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:))
+			plot3([minx minx], [maxy maxy], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:))
+			plot3([maxx maxx], [maxy maxy], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:))
+
+			% detector
+			D1 = v(a,4:6) - d1/2*v(a,7:9) - d2/2*v(a,10:12);
+			D2 = v(a,4:6) + d1/2*v(a,7:9) - d2/2*v(a,10:12);
+			D3 = v(a,4:6) + d1/2*v(a,7:9) + d2/2*v(a,10:12);
+			D4 = v(a,4:6) - d1/2*v(a,7:9) + d2/2*v(a,10:12);
+			plot3([D1(1) D2(1) D3(1) D4(1) D1(1)], [D1(2) D2(2) D3(2) D4(2) D1(2)], [D1(3) D2(3) D3(3) D4(3) D1(3)], 'LineWidth', 2, 'Color', colours(5,:))
+
+			% ray
+			s = maxx - minx;
+			plot3([0 v(a,1)]*s*0.30, [0 v(a,2)]*s*0.30, [0 v(a,3)]*s*0.30, 'LineWidth', 2, 'Color', colours(3,:))
+
+		elseif strcmp(proj_geom.type,'cone_vec')
+
+			v = proj_geom.Vectors;
+			d1 = proj_geom.DetectorRowCount;
+			d2 = proj_geom.DetectorColCount;
+
+			if ~isfield(vol_geom, 'option')
+				minx = -vol_geom.GridRowCount/2;
+				miny = -vol_geom.GridColCount/2;
+				minz = -vol_geom.GridSliceCount/2;
+				maxx = vol_geom.GridRowCount/2;
+				maxy = vol_geom.GridColCount/2;
+				maxz = vol_geom.GridSliceCount/2;
+			else
+				minx = vol_geom.option.WindowMinX;
+				miny = vol_geom.option.WindowMinY;
+				minz = vol_geom.option.WindowMinZ;
+				maxx = vol_geom.option.WindowMaxX;
+				maxy = vol_geom.option.WindowMaxY;
+				maxz = vol_geom.option.WindowMaxZ;
+			end
+
+			% axis
+			windowminx = min(v(:,4));
+			windowminy = min(v(:,5));
+			windowminz = max(v(:,6));
+			windowmaxx = max(v(:,4));
+			windowmaxy = max(v(:,5));
+			windowmaxz = max(v(:,6));
+			cla
+			axis([minx maxx miny maxy minz maxz]*5.10)
+
+			% volume
+			plot3([minx minx maxx maxx minx], [miny maxy maxy miny miny], [minz minz minz minz minz], 'LineWidth', 1, 'Color', colours(1,:))
+			plot3([minx minx maxx maxx minx], [miny maxy maxy miny miny], [maxz maxz maxz maxz maxz], 'LineWidth', 1, 'Color', colours(1,:))
+			plot3([minx minx], [miny miny], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:))
+			plot3([maxx maxx], [miny miny], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:))
+			plot3([minx minx], [maxy maxy], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:))
+			plot3([maxx maxx], [maxy maxy], [minz maxz], 'LineWidth', 1, 'Color', colours(1,:))
+
+			% detector
+			D1 = v(a,4:6) - d1/2*v(a,7:9) - d2/2*v(a,10:12);
+			D2 = v(a,4:6) + d1/2*v(a,7:9) - d2/2*v(a,10:12);
+			D3 = v(a,4:6) + d1/2*v(a,7:9) + d2/2*v(a,10:12);
+			D4 = v(a,4:6) - d1/2*v(a,7:9) + d2/2*v(a,10:12);
+			plot3([D1(1) D2(1) D3(1) D4(1) D1(1)], [D1(2) D2(2) D3(2) D4(2) D1(2)], [D1(3) D2(3) D3(3) D4(3) D1(3)], 'LineWidth', 2, 'Color', colours(5,:))
+
+			% beam
+			plot3([v(a,1) D1(1)], [v(a,2) D1(2)], [v(a,3) D1(3)], 'LineWidth', 1, 'Color', colours(3,:))
+			plot3([v(a,1) D2(1)], [v(a,2) D2(2)], [v(a,3) D2(3)], 'LineWidth', 1, 'Color', colours(3,:))
+			plot3([v(a,1) D3(1)], [v(a,2) D3(2)], [v(a,3) D3(3)], 'LineWidth', 1, 'Color', colours(3,:))
+			plot3([v(a,1) D4(1)], [v(a,2) D4(2)], [v(a,3) D4(3)], 'LineWidth', 1, 'Color', colours(3,:))
+
+
+		else
+			error('invalid projector type')
+
+		end
+    end
+
+end