summaryrefslogtreecommitdiffstats
path: root/matlab
diff options
context:
space:
mode:
Diffstat (limited to 'matlab')
-rw-r--r--matlab/algorithms/DART/IterativeTomography.m2
-rw-r--r--matlab/algorithms/DART/examples/example1.m53
-rw-r--r--matlab/mex/astra_mex_data2d_c.cpp5
-rw-r--r--matlab/tools/astra_create_proj_geom.m13
-rw-r--r--matlab/tools/astra_geom_2vec.m76
-rw-r--r--matlab/tools/astra_geom_postalignment.m33
-rw-r--r--matlab/tools/astra_geom_size.m36
-rw-r--r--matlab/tools/astra_geom_visualize.m216
8 files changed, 379 insertions, 55 deletions
diff --git a/matlab/algorithms/DART/IterativeTomography.m b/matlab/algorithms/DART/IterativeTomography.m
index 3875e6b..c3e737f 100644
--- a/matlab/algorithms/DART/IterativeTomography.m
+++ b/matlab/algorithms/DART/IterativeTomography.m
@@ -77,7 +77,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 cb02e0f..791e440 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 909d229..6863b7a 100644
--- a/matlab/mex/astra_mex_data2d_c.cpp
+++ b/matlab/mex/astra_mex_data2d_c.cpp
@@ -45,6 +45,7 @@ $Id$
#include "astra/SparseMatrixProjectionGeometry2D.h"
#include "astra/FanFlatProjectionGeometry2D.h"
#include "astra/FanFlatVecProjectionGeometry2D.h"
+#include "astra/ParallelVecProjectionGeometry2D.h"
using namespace std;
using namespace astra;
@@ -160,6 +161,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();
}
@@ -449,6 +452,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 862e410..ff7d74d 100644
--- a/matlab/tools/astra_create_proj_geom.m
+++ b/matlab/tools/astra_create_proj_geom.m
@@ -107,6 +107,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
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 16:45:20 +0000