mirror of
https://github.com/c-sooyoung/fold_slice.git
synced 2026-09-17 19:39:08 +09:00
539 lines
19 KiB
Matlab
539 lines
19 KiB
Matlab
function [img,H] = iradonfast_v3(varargin)
|
|
%IRADON Inverse Radon transform.
|
|
% I = iradon(R,THETA) reconstructs the image I from projection data in the
|
|
% 2-D array R. The columns of R are parallel beam projection data.
|
|
% IRADON assumes that the center of rotation is the center point of the
|
|
% projections, which is defined as ceil(size(R,1)/2).
|
|
%
|
|
% THETA describes the angles (in degrees) at which the projections were
|
|
% taken. It can be either a vector containing the angles or a scalar
|
|
% specifying D_theta, the incremental angle between projections. If THETA
|
|
% is a vector, it must contain angles with equal spacing between them. If
|
|
% THETA is a scalar specifying D_theta, the projections are taken at
|
|
% angles THETA = m * D_theta; m = 0,1,2,...,size(R,2)-1. If the input is
|
|
% the empty matrix ([]), D_theta defaults to 180/size(R,2).
|
|
%
|
|
% IRADON uses the filtered backprojection algorithm to perform the inverse
|
|
% Radon transform. The filter is designed directly in the frequency
|
|
% domain and then multiplied by the FFT of the projections. The
|
|
% projections are zero-padded to a power of 2 before filtering to prevent
|
|
% spatial domain aliasing and to speed up the FFT.
|
|
%
|
|
% I = IRADON(R,THETA,INTERPOLATION,FILTER,FREQUENCY_SCALING,OUTPUT_SIZE)
|
|
% specifies parameters to use in the inverse Radon transform. You can
|
|
% specify any combination of the last four arguments. IRADON uses default
|
|
% values for any of these arguments that you omit.
|
|
%
|
|
% INTERPOLATION specifies the type of interpolation to use in the
|
|
% backprojection. The default is linear interpolation. Available methods
|
|
% are:
|
|
%
|
|
% 'nearest' - nearest neighbor interpolation
|
|
% 'linear' - linear interpolation (default)
|
|
% 'spline' - spline interpolation
|
|
% 'pchip' - shape-preserving piecewise cubic interpolation
|
|
% 'cubic' - same as 'pchip'
|
|
% 'v5cubic' - the cubic interpolation from MATLAB 5, which does not
|
|
% extrapolate and uses 'spline' if X is not equally spaced.
|
|
%
|
|
% FILTER specifies the filter to use for frequency domain filtering.
|
|
% FILTER is a string that specifies any of the following standard filters:
|
|
%
|
|
% 'Ram-Lak' The cropped Ram-Lak or ramp filter (default). The
|
|
% frequency response of this filter is |f|. Because this
|
|
% filter is sensitive to noise in the projections, one of
|
|
% the filters listed below may be preferable.
|
|
% 'Shepp-Logan' The Shepp-Logan filter multiplies the Ram-Lak filter by
|
|
% a sinc function.
|
|
% 'Cosine' The cosine filter multiplies the Ram-Lak filter by a
|
|
% cosine function.
|
|
% 'Hamming' The Hamming filter multiplies the Ram-Lak filter by a
|
|
% Hamming window.
|
|
% 'Hann' The Hann filter multiplies the Ram-Lak filter by a
|
|
% Hann window.
|
|
% 'parzen' The parzen filter multiplies the Ram-Lak filter by a
|
|
% Parzen window. Guizar - Nov 30 2010
|
|
%
|
|
% FREQUENCY_SCALING is a scalar in the range (0,1] that modifies the
|
|
% filter by rescaling its frequency axis. The default is 1. If
|
|
% FREQUENCY_SCALING is less than 1, the filter is compressed to fit into
|
|
% the frequency range [0,FREQUENCY_SCALING], in normalized frequencies;
|
|
% all frequencies above FREQUENCY_SCALING are set to 0.
|
|
%
|
|
% OUTPUT_SIZE is a scalar that specifies the number of rows and columns in
|
|
% the reconstructed image. If OUTPUT_SIZE is not specified, the size is
|
|
% determined from the length of the projections:
|
|
%
|
|
% OUTPUT_SIZE = 2*floor(size(R,1)/(2*sqrt(2)))
|
|
%
|
|
% If you specify OUTPUT_SIZE, IRADON reconstructs a smaller or larger
|
|
% portion of the image, but does not change the scaling of the data.
|
|
%
|
|
% If the projections were calculated with the RADON function, the
|
|
% reconstructed image may not be the same size as the original image.
|
|
%
|
|
% [I,H] = iradon(...) returns the frequency response of the filter in the
|
|
% vector H.
|
|
%
|
|
% Class Support
|
|
% -------------
|
|
% R can be double or single. All other numeric input arguments must be double.
|
|
% I has the same class as R. H is double.
|
|
%
|
|
% Example
|
|
% -------
|
|
% P = phantom(128);
|
|
% R = radon(P,0:179);
|
|
% I = iradon(R,0:179,'nearest','Hann');
|
|
% figure, imshow(P), figure, imshow(I);
|
|
%
|
|
% See also FAN2PARA, FANBEAM, IFANBEAM, PARA2FAN, PHANTOM, RADON.
|
|
%
|
|
% Copyright 1993-2004 The MathWorks, Inc.
|
|
% $Revision: 1.1 $ $Date: 2013/10/31 13:29:18 $
|
|
%
|
|
% References:
|
|
% A. C. Kak, Malcolm Slaney, "Principles of Computerized Tomographic
|
|
% Imaging", IEEE Press 1988.
|
|
%
|
|
% 'derivative' - Optional input argument to use the filter for input
|
|
% derivative of projections - Manuel Guizar - Nov 30 2010
|
|
%
|
|
% Weights for uneven angles
|
|
% If the angles are between 0 and 180 degrees it computes custom weights to
|
|
% the projections in order to allow for not equal angular sampling. Other
|
|
% angles are not considered because it would need significantly more
|
|
% checks.
|
|
% If you don't want to use this functionality, add 360 to your angles
|
|
% (theta+360)
|
|
% Manuel Guizar - Oct 10 2015
|
|
|
|
import utils.*
|
|
|
|
[p,theta,filter,d,interp,N,derivative] = parse_inputs(varargin{:});
|
|
|
|
% use Matlab or C-code for the linear interpolation
|
|
use_original_matlab_code = 0;
|
|
determine_weights = 1;
|
|
|
|
%%% Determine weights for uneven angular sampling %%%
|
|
if any(theta<0)
|
|
warning('There are some theta < 0 angles. Using constant angular sampling code.')
|
|
determine_weights = 0;
|
|
end
|
|
if any(theta>=pi)
|
|
warning('There are some theta >= 180 angles. Using constant angular sampling code.')
|
|
determine_weights = 0;
|
|
end
|
|
if any(diff(theta)==0)
|
|
warning('There are some repeated angles. Using constant angular sampling code.')
|
|
determine_weights = 0;
|
|
end
|
|
|
|
if abs(max(theta)-min(theta)-pi) > 5*mean(diff(sort(theta)))
|
|
warning('Missing wedge is to large for weighting')
|
|
determine_weights = 0;
|
|
end
|
|
if verLessThan('matlab', '9.1')
|
|
error('Only matlab versions >= 2016b are supported')
|
|
end
|
|
if determine_weights
|
|
[theta,ind_sort] = sort(theta);
|
|
weights = theta*0;
|
|
weights(2:end-1) = - theta(1:end-2)/2 + theta(3:end)/2;
|
|
weights(1) = - (-pi + theta(end))/2 + theta(2)/2;
|
|
weights(end) = - theta(end-1)/2 + (pi+theta(1))/2;
|
|
p = p(:,ind_sort,:);
|
|
p = bsxfun(@times, p, weights/2);
|
|
else
|
|
p = p*pi/(2*length(theta));
|
|
end
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
p = single(p); % single will make it ~2x faster
|
|
|
|
|
|
% Design the filter
|
|
len=size(p,1);
|
|
H = designFilter(filter, len, d, derivative);
|
|
p(length(H),1)=0; % Zero pad projections
|
|
|
|
% In the code below, I continuously reuse the array p so as to
|
|
% save memory. This makes it harder to read, but the comments
|
|
% explain what is going on.
|
|
|
|
p = fft(p); % p holds fft of projections
|
|
|
|
p = p.*H; % frequency domain filtering
|
|
|
|
p = real(ifft(p)); % p is the filtered projections
|
|
p(len+1:end,:,:) = []; % Truncate the filtered projections
|
|
|
|
% Define the x & y axes for the reconstructed image so that the origin
|
|
% (center) is in the spot which RADON would choose.
|
|
center = floor((N + 1)/2);
|
|
xleft = -center + 1;
|
|
x = (1:N) - 1 + xleft;
|
|
x = repmat(x, N, 1);
|
|
|
|
ytop = center - 1;
|
|
y = (N:-1:1).' - N + ytop;
|
|
y = repmat(y, 1, N);
|
|
|
|
if ((~strcmp(interp, 'linear')) || (use_original_matlab_code))
|
|
costheta = cos(theta);
|
|
sintheta = sin(theta);
|
|
img = zeros(N,class(p)); % Allocate memory for the image.
|
|
end
|
|
|
|
ctrIdx = ceil(len/2); % index of the center of the projections
|
|
|
|
% Zero pad the projections to size 1+2*ceil(N/sqrt(2)) if this
|
|
% quantity is greater than the length of the projections
|
|
imgDiag = 2*ceil(N/sqrt(2))+1; % largest distance through image.
|
|
|
|
if size(p,1) < imgDiag
|
|
rz = imgDiag - size(p,1); % how many rows of zeros
|
|
p = [zeros(ceil(rz/2),size(p,2),size(p,3),class(p)); p; zeros(floor(rz/2),size(p,2),size(p,3), class(p))];
|
|
ctrIdx = ctrIdx+ceil(rz/2);
|
|
end
|
|
|
|
Nlayers = size(p,3);
|
|
|
|
|
|
|
|
% Backprojection - vectorized in (x,y), looping over theta
|
|
switch interp
|
|
case 'nearest neighbor'
|
|
|
|
for i=1:length(theta)
|
|
proj = p(:,i);
|
|
t = round(x*costheta(i) + y*sintheta(i));
|
|
img = img + proj(t+ctrIdx);
|
|
end
|
|
|
|
case 'linear'
|
|
if (use_original_matlab_code)
|
|
for i=1:length(theta)
|
|
proj = p(:,i);
|
|
t = x.*costheta(i) + y.*sintheta(i);
|
|
a = floor(t);
|
|
img = img + (t-a).*proj(a+1+ctrIdx) + (a+1-t).*proj(a+ctrIdx);
|
|
end
|
|
else
|
|
if gpuDeviceCount
|
|
try
|
|
img = iradon_cuda( p, theta, x, y );
|
|
catch err
|
|
warning('GPU failed, running slower CPU version\n Error message: %s', err.message)
|
|
end
|
|
end
|
|
if ~exist('img', 'var')
|
|
img = zeros(N,N,Nlayers, 'single');
|
|
p = double(p);
|
|
for i = 1:size(p,3)
|
|
progressbar(i,Nlayers)
|
|
img(:,:,i) = iradon_c( p(:,:,i), theta, x, y );
|
|
end
|
|
end
|
|
end
|
|
|
|
|
|
case {'spline','pchip','cubic','v5cubic'}
|
|
|
|
interp_method = sprintf('*%s',interp); % Add asterisk to assert
|
|
% even-spacing of taxis
|
|
|
|
for i=1:length(theta)
|
|
proj = p(:,i);
|
|
taxis = (1:size(p,1)) - ctrIdx;
|
|
t = x.*costheta(i) + y.*sintheta(i);
|
|
projContrib = interp1(taxis,proj,t(:),interp_method);
|
|
img = img + reshape(projContrib,N,N);
|
|
end
|
|
|
|
end
|
|
|
|
% img = img*pi/(2*length(theta));
|
|
end
|
|
|
|
%%%
|
|
%%% Sub-Function: designFilter
|
|
%%%
|
|
|
|
function filt = designFilter(filter, len, d, derivative)
|
|
% Returns the Fourier Transform of the filter which will be
|
|
% used to filter the projections
|
|
%
|
|
% INPUT ARGS: filter - either the string specifying the filter
|
|
% len - the length of the projections
|
|
% d - the fraction of frequencies below the nyquist
|
|
% which we want to pass
|
|
%
|
|
% OUTPUT ARGS: filt - the filter to use on the projections
|
|
|
|
|
|
order = max(64,2^nextpow2(2*len));
|
|
|
|
% First create a ramp filter - go up to the next highest
|
|
% power of 2.
|
|
if derivative
|
|
filt = 0*( 0:(order/2) )+1;
|
|
else
|
|
filt = 2*( 0:(order/2) )./order;
|
|
end
|
|
w = 2*pi*(0:size(filt,2)-1)/order; % frequency axis up to Nyquist
|
|
|
|
switch filter
|
|
case 'ram-lak'
|
|
% Do nothing
|
|
case 'shepp-logan'
|
|
% be careful not to divide by 0:
|
|
filt(2:end) = filt(2:end) .* (sin(w(2:end)/(2*d))./(w(2:end)/(2*d)));
|
|
case 'cosine'
|
|
filt(2:end) = filt(2:end) .* cos(w(2:end)/(2*d));
|
|
case 'hamming'
|
|
filt(2:end) = filt(2:end) .* (.54 + .46 * cos(w(2:end)/d));
|
|
case 'hann'
|
|
filt(2:end) = filt(2:end) .*(1+cos(w(2:end)./d)) / 2;
|
|
case 'parzen'
|
|
aux = parzenwin(round(2*size(filt,2)*d)-1)';
|
|
aux = aux(round(size(aux,2)/2):round(size(aux,2)));
|
|
filt(1:size(aux,2)) = filt(1:size(aux,2)).*aux;
|
|
filt(size(aux,2)+1:end) = 0;
|
|
otherwise
|
|
eid = sprintf('Images:%s:invalidFilter',mfilename);
|
|
msg = 'Invalid filter selected.';
|
|
error(eid,'%s',msg);
|
|
end
|
|
|
|
filt(w>pi*d) = 0; % Crop the frequency response
|
|
if derivative
|
|
filt = [filt' ; -filt(end-1:-1:2)']/(1i*pi); % Symmetry of the filter
|
|
else
|
|
filt = [filt' ; filt(end-1:-1:2)']; % Symmetry of the filter
|
|
end
|
|
end
|
|
|
|
%%%
|
|
%%% Sub-Function: parse_inputs
|
|
%%%
|
|
|
|
function [p,theta,filter,d,interp,N,derivative] = parse_inputs(varargin)
|
|
% Parse the input arguments and retun things
|
|
%
|
|
% Inputs: varargin - Cell array containing all of the actual inputs
|
|
%
|
|
% Outputs: p - Projection data
|
|
% theta - the angles at which the projections were taken
|
|
% filter - string specifying filter or the actual filter
|
|
% d - a scalar specifying normalized freq. at which to crop
|
|
% the frequency response of the filter
|
|
% interp - the type of interpolation to use
|
|
% N - The size of the reconstructed image
|
|
|
|
if nargin<2
|
|
eid = sprintf('Images:%s:tooFewInputs',mfilename);
|
|
msg = 'Invalid input arguments.';
|
|
error(eid,'%s',msg);
|
|
end
|
|
|
|
p = varargin{1};
|
|
theta = pi*varargin{2}/180;
|
|
|
|
% Default values
|
|
N = 0; % Size of the reconstructed image
|
|
d = 1; % Defaults to no cropping of filters frequency response
|
|
filter = 'ram-lak'; % The ramp filter is the default
|
|
interp = 'linear'; % default interpolation is linear
|
|
string_args = {'nearest neighbor', 'linear', 'spline', 'pchip', 'cubic', 'v5cubic', ...
|
|
'ram-lak','shepp-logan','cosine','hamming', 'hann','parzen','derivative'};
|
|
|
|
for i=3:nargin
|
|
arg = varargin{i};
|
|
if ischar(arg)
|
|
idx = strmatch(lower(arg),string_args);
|
|
if isempty(idx)
|
|
eid = sprintf('Images:%s:unknownInputString',mfilename);
|
|
msg = sprintf('Unknown input string: %s.', arg);
|
|
error(eid,'%s',msg);
|
|
elseif numel(idx) > 1
|
|
eid = sprintf('Images:%s:ambiguousInputString',mfilename);
|
|
msg = sprintf('Ambiguous input string: %s.', arg);
|
|
error(eid,'%s',msg);
|
|
elseif numel(idx) == 1
|
|
if idx <= 6 % It is the interpolation
|
|
interp = string_args{idx};
|
|
elseif (idx > 6) && (idx <= 12)
|
|
filter = string_args{idx};
|
|
elseif idx == 13
|
|
derivative = true; % Input is a derivative of sinogram
|
|
end
|
|
end
|
|
elseif numel(arg)==1
|
|
if arg <=1
|
|
d = arg;
|
|
else
|
|
N = arg;
|
|
end
|
|
else
|
|
eid = sprintf('Images:%s:invalidInputParameters',mfilename);
|
|
msg = 'Invalid input parameters';
|
|
error(eid,'%s',msg);
|
|
end
|
|
end
|
|
|
|
% If the user didn't specify the size of the reconstruction, so
|
|
% deduce it from the length of projections
|
|
if N==0
|
|
N = 2*floor( size(p,1)/(2*sqrt(2)) ); % This doesn't always jive with RADON
|
|
end
|
|
|
|
% for empty theta, choose an intelligent default delta-theta
|
|
if isempty(theta)
|
|
theta = pi / size(p,2);
|
|
end
|
|
|
|
% If the user passed in delta-theta, build the vector of theta values
|
|
if numel(theta)==1
|
|
theta = (0:(size(p,2)-1))* theta;
|
|
end
|
|
|
|
if length(theta) ~= size(p,2)
|
|
eid = sprintf('Images:%s:thetaNotMatchingProjectionNumber',mfilename);
|
|
msg = 'THETA does not match the number of projections.';
|
|
error(eid,'%s',msg);
|
|
end
|
|
|
|
if ~exist('derivative')
|
|
derivative = false;
|
|
end
|
|
end
|
|
|
|
function img = iradon_cuda( sinogram_full, theta, x, y )
|
|
%% preprocess data for the ASTRA toolbox wrapper
|
|
import utils.*
|
|
|
|
[Wsin, Nproj, Nlayers] = size(sinogram_full);
|
|
[Nx, Ny] = size(x);
|
|
|
|
if any(size(sinogram_full) > 4096)
|
|
error('Maximal size of CUDA 3D texture is 4096x4096x4096, size of sinogram is %ix%ix%i', size(sinogram_full))
|
|
end
|
|
|
|
%% create data and geometry
|
|
Nangles = length(theta);
|
|
assert(Nproj == Nangles, 'Wrong input size')
|
|
|
|
vectors_all = zeros(Nangles, 12);
|
|
for i = 1:Nangles
|
|
% ray direction
|
|
vectors_all(i,1) = sin(theta(i));
|
|
vectors_all(i,2) = -cos(theta(i));
|
|
vectors_all(i,3) = 0;
|
|
vectors_all(i,1:3) = vectors_all(i,1:3);
|
|
% center of detector
|
|
vectors_all(i,4:6) = 0;
|
|
% vector from detector pixel (0,0) to (0,1)
|
|
vectors_all(i,7) = cos(theta(i));
|
|
vectors_all(i,8) = sin(theta(i));
|
|
vectors_all(i,9) = 0;
|
|
% vector from detector pixel (0,0) to (1,0)
|
|
vectors_all(i,10) = 0;
|
|
vectors_all(i,11) = 0;
|
|
vectors_all(i,12) = 1;
|
|
end
|
|
%% astra settings
|
|
cfg.iVolX = Nx;
|
|
cfg.iVolY = Ny;
|
|
cfg.iVolZ = Nlayers;
|
|
cfg.iProjAngles = Nangles;
|
|
cfg.iProjU = Wsin;
|
|
cfg.iProjV = Nlayers;
|
|
cfg.iRaysPerDet = 1;
|
|
cfg.iRaysPerDetDim = 1;
|
|
cfg.iRaysPerVoxelDim = 1;
|
|
|
|
|
|
|
|
|
|
|
|
%% call ASTRA wrapper
|
|
gpu = gpuDevice;
|
|
AvailableMemory = gpu.AvailableMemory;
|
|
% required memory is 2x dataset size, max data size allowed by GPU is 1024MB
|
|
% max number of GPU projections is 1024 (limit in ASTRA code, can be fixed ... )
|
|
Nangular_slices = ceil(Nangles/1024); % limitation in the ASTRA code
|
|
Nslices = ceil(numel(sinogram_full)*8 / min(1024e6, AvailableMemory));
|
|
|
|
if Nangular_slices > 1 || Nslices > 1
|
|
if Nslices > 1
|
|
fprintf('Dataset does not fit to GPU memory => autosplitting\nFree memory: %iMB\tDataset size: %iMB\n\n', ceil(AvailableMemory/1e6), ceil(numel(sinogram_full)*4/1e6))
|
|
else
|
|
fprintf('Number of angles axceeded limit 1024 => autosplitting\nFree memory: %iMB\tDataset size: %iMB\n\n', ceil(AvailableMemory/1e6), ceil(numel(sinogram_full)*4/1e6))
|
|
end
|
|
img = zeros([Nx, Ny, Nlayers], 'single');
|
|
end
|
|
for j = 1:Nangular_slices % split in the angular space
|
|
|
|
|
|
ind_angle = 1+(j-1)*ceil(Nangles/Nangular_slices):min(Nangles, j*ceil(Nangles/Nangular_slices));
|
|
vectors = vectors_all(ind_angle,:);
|
|
|
|
%% subpixel shift correction to make it consistent with the matlab iradon code
|
|
|
|
cfg.iProjAngles = length(ind_angle);
|
|
% split the sinogram (avoid copying if possible)
|
|
if Nangular_slices > 1
|
|
sinogram = sinogram_full(:,ind_angle,:);
|
|
else
|
|
sinogram = sinogram_full;
|
|
end
|
|
if Nslices == 1
|
|
%%%% apply geometry correction to shift reconstruction into center %%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
vectors(:,4:6) = vectors(:,4:6) -(vectors(:,10:12)*cfg.iProjV/2+vectors(:,7:9)*cfg.iProjU/2);
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% reshape the inputs for ASTRA
|
|
sinogram = reshape(sinogram,[cfg.iProjU,cfg.iProjV,cfg.iProjAngles]);
|
|
sinogram = gpuArray(single(sinogram));
|
|
%% ASTRA WRAPPER
|
|
vol = astra.iradon_gpu_wrapper(sinogram, cfg, vectors);
|
|
vol = gather(vol); % backprojection
|
|
if Nangular_slices > 1
|
|
img = img + vol;
|
|
else
|
|
img = vol;
|
|
end
|
|
else % low memory case (a bit slower) => split to several chunks rotation axis
|
|
for i = 1:Nslices
|
|
progressbar(i,Nslices+1)
|
|
ind = 1+(i-1)*ceil(Nlayers/Nslices):min(Nlayers, i*ceil(Nlayers/Nslices));
|
|
p_tmp = sinogram(:,:,ind);
|
|
cfg.iVolZ = length(ind);
|
|
cfg.iProjV = length(ind);
|
|
p_tmp = reshape(p_tmp,[cfg.iProjU,cfg.iProjV,cfg.iProjAngles]);
|
|
p_tmp = gpuArray(single(p_tmp));
|
|
vectors_all_tmp = vectors;
|
|
|
|
%%%% apply geometry correction to shift reconstruction into center
|
|
vectors_all_tmp(:,4:6) = vectors_all_tmp(:,4:6) -(vectors_all_tmp(:,10:12)*cfg.iProjV/2+vectors_all_tmp(:,7:9)*cfg.iProjU/2);
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
|
|
%% ASTRA WRAPPER
|
|
vol = astra.iradon_gpu_wrapper(p_tmp, cfg, vectors_all_tmp);
|
|
vol = gather(vol); % gather from GPU takes 30% of time !!!
|
|
if Nangular_slices > 1
|
|
img(:,:,ind) = img(:,:,ind) + vol; % adding up to the array is also pretty slow
|
|
else
|
|
img(:,:,ind) = vol;
|
|
end
|
|
end
|
|
progressbar(i,Nslices)
|
|
end
|
|
end
|
|
% rotate as the output for Matlab
|
|
img = rot90(img, 1);
|
|
|
|
end
|
|
|