mirror of
https://github.com/c-sooyoung/fold_slice.git
synced 2026-09-17 19:39:08 +09:00
171 lines
8.4 KiB
Matlab
171 lines
8.4 KiB
Matlab
% FBP_ZSPLIT filtered back projection with simple splitting along vertical direction -> works only for classical tomography, not for laminography
|
||
%
|
||
% [rec] = FBP_zsplit(sinogram, cfg, vectors, varargin)
|
||
%
|
||
% Inputs:
|
||
% **sino - sinogram (Nlayers x width x Nangles)
|
||
% **cfg - config struct from astra.ASTRA_initialize
|
||
% **vectors - vectors of projection rotation generated by astra.ASTRA_initialize
|
||
% *optional*
|
||
% ** split =[1,1,1] - split the solved volume, split(3 is used to split in separated blocks, split(1:2) is used inside Atx_partial to du subplitting for ASTRA
|
||
% ** valid_angles = [] - list of valid angles, []==all are valid
|
||
% ** filter = 'ram-lak' - name of the FBP filter
|
||
% ** filter_value = 1 - fitlering value for the FBP filter
|
||
% ** deformation_fields = {} - cell 3x1 of deformation arrays
|
||
% ** GPU = [] - list of GPUs to be used in reconstruction
|
||
% ** split_sub = [1,1,1] - splitting of the sub block on smaller tasks in the Atx_partial method , 1 == no splitting
|
||
% ** verbose = 1 - verbose = 0 : quiet, verbose : standard info , verbose = 2: debug
|
||
% ** use_derivative = false - calculate reconstruction from the phase derivative
|
||
% ** extra_padding = false - surround the projection by void space to enforce zero around tomogram
|
||
% ** keep_on_GPU - if false, move the reconstruction back from GPU before returning
|
||
% ** determine_weights = true - reweight projections if the angles are not equidistant
|
||
% ** mask = [] - apply mask on reconstruction , inputs is 2D or 3D array
|
||
% ** padding = 0 - zero padding is improving standard tomography. 'symmetric' is good for lamino / local tomo
|
||
% ** only_filter_sinogram = false - return filtered sinogram, do not backproject
|
||
%
|
||
% *returns*
|
||
% ++rec - reconstructed volume
|
||
|
||
|
||
%*-----------------------------------------------------------------------*
|
||
%| |
|
||
%| Except where otherwise noted, this work is licensed under a |
|
||
%| Creative Commons Attribution-NonCommercial-ShareAlike 4.0 |
|
||
%| International (CC BY-NC-SA 4.0) license. |
|
||
%| |
|
||
%| Copyright (c) 2017 by Paul Scherrer Institute (http://www.psi.ch) |
|
||
%| |
|
||
%| Author: CXS group, PSI |
|
||
%*-----------------------------------------------------------------------*
|
||
% You may use this code with the following provisions:
|
||
%
|
||
% If the code is fully or partially redistributed, or rewritten in another
|
||
% computing language this notice should be included in the redistribution.
|
||
%
|
||
% If this code, or subfunctions or parts of it, is used for research in a
|
||
% publication or if it is fully or partially rewritten for another
|
||
% computing language the authors and institution should be acknowledged
|
||
% in written form in the publication: “Data processing was carried out
|
||
% using the “cSAXS matlab package” developed by the CXS group,
|
||
% Paul Scherrer Institut, Switzerland.”
|
||
% Variations on the latter text can be incorporated upon discussion with
|
||
% the CXS group if needed to more specifically reflect the use of the package
|
||
% for the published work.
|
||
%
|
||
% A publication that focuses on describing features, or parameters, that
|
||
% are already existing in the code should be first discussed with the
|
||
% authors.
|
||
%
|
||
% This code and subroutines are part of a continuous development, they
|
||
% are provided “as they are” without guarantees or liability on part
|
||
% of PSI or the authors. It is the user responsibility to ensure its
|
||
% proper use and the correctness of the results.
|
||
|
||
function [rec] = FBP_zsplit(sinogram, cfg0, vectors0, varargin)
|
||
|
||
par = inputParser;
|
||
par.KeepUnmatched = true;
|
||
par.addOptional('split', [1,1,1,1]) % split the solved volume, split(3) is used to split in separated blocks, split(1:2) is used inside Atx_partial to du subplitting for ASTRA
|
||
par.addOptional('keep_on_GPU', false) % keep results in GPU
|
||
par.addOptional('verbose', true) % verbosity level
|
||
par.addOptional('mask', []) % apply reconstruction mask, inputs is 2D or 3D array
|
||
par.addOptional('GPU', []) % ids of the used GPUs
|
||
par.addOptional('use_GPU', true, @islogical) % if false, use CPU based reconstruction
|
||
|
||
par.parse(varargin{:})
|
||
r = par.Results;
|
||
|
||
|
||
[Nlayers,Nw,Nproj] = size(sinogram);
|
||
cfg0.iProjAngles = Nproj;
|
||
assert(cfg0.iProjU == Nw, 'Wrong sinogram width')
|
||
assert(cfg0.iProjV == Nlayers, 'Wrong sinogram height')
|
||
assert(mod(Nw,2)==0, 'Only even width of sinogram is supported')
|
||
assert(length(r.GPU)<=1, 'For multiGPU use tomo.FBP function')
|
||
|
||
if r.verbose > 0; fprintf('====== FBP split ==========\n'); end
|
||
|
||
|
||
Nelements = numel(sinogram) ;
|
||
|
||
if gpuDeviceCount
|
||
Nblocks = r.split(3);
|
||
r.split(3) = 1;
|
||
gpu = gpuDevice;
|
||
% empirical condition, may be too pesimistic
|
||
Nblocks = max(Nblocks,ceil( (2*8*4* Nelements) / gpu.AvailableMemory)) ;
|
||
Nblocks = max(Nblocks, Nelements/ double(intmax('int32')));
|
||
Nblocks = max(Nblocks, length(r.GPU));
|
||
else
|
||
% CPU processing
|
||
max_block_size = min(utils.check_available_memory*1e6, 20e9); %% work with 10GB blocks
|
||
Nblocks = ceil( (6*8* Nelements) / max_block_size) ;
|
||
end
|
||
|
||
if ~r.use_GPU
|
||
rec = tomo.FBP_CPU(sinogram, cfg0, vectors0, varargin{:}, 'verbose', 0);
|
||
return
|
||
end
|
||
|
||
if Nblocks == 1 || cfg0.iProjU>4096 || cfg0.iProjV>4096
|
||
% for small datatsets process everthing in a single block
|
||
if Nblocks == 1
|
||
sinogram = utils.Garray(sinogram);
|
||
end
|
||
rec = tomo.FBP(sinogram, cfg0, vectors0, varargin{:}, 'verbose', 0);
|
||
if ~isempty(r.mask)
|
||
rec = rec .* r.mask;
|
||
end
|
||
if ~r.keep_on_GPU
|
||
rec = gather(rec);
|
||
end
|
||
else
|
||
%% split tomogram into vertical blocks, it is the most effecient way how to calculate it
|
||
% -> upload each block to GPU and keep it there for maximal speed
|
||
Nl_small = ceil(Nlayers / Nblocks);
|
||
% keep on GPU only if the volume is small enough !!
|
||
r.keep_on_GPU = r.keep_on_GPU && (4*cfg0.iVolX*cfg0.iVolY*cfg0.iVolZ) < gpu.AvailableMemory / 10;
|
||
|
||
if r.keep_on_GPU && gpuDeviceCount
|
||
rec = gpuArray.zeros(cfg0.iVolX, cfg0.iVolY, cfg0.iVolZ, 'single');
|
||
else
|
||
rec = zeros(cfg0.iVolX, cfg0.iVolY, cfg0.iVolZ, 'single');
|
||
end
|
||
for ii = 1:Nblocks
|
||
if r.verbose>0; utils.progressbar(ii, Nblocks); end
|
||
ind = 1+Nl_small*(ii-1) : min(Nlayers, Nl_small*ii);
|
||
if isempty(ind); continue; end
|
||
if isa(sinogram, 'gpuArray')
|
||
% if on GPU, use matlab memcpy
|
||
sinogram_small = sinogram(ind,:,:);
|
||
else
|
||
% fast MEX memory copy
|
||
pos = [Nl_small*(ii-1),0];
|
||
sinogram_small = zeros(length(ind),Nw,Nproj, 'like', sinogram);
|
||
sinogram_small = utils.get_from_3D_projection(sinogram_small, sinogram, repmat(pos,Nproj,1), 1:Nproj);
|
||
end
|
||
cfg = cfg0;
|
||
vectors = vectors0;
|
||
cfg.iProjV = size(sinogram_small,1);
|
||
cfg.iVolZ = cfg.iProjV;
|
||
vectors(:,4:6) = vectors(:,4:6) - vectors(:,10:12)*(cfg.iProjV - cfg0.iProjV)/2;
|
||
sinogram_small = utils.Garray(fp16.get(sinogram_small));
|
||
|
||
%% call standard FBP on the data that are on GPU
|
||
try
|
||
rec_tmp = tomo.FBP(sinogram_small, cfg, vectors,varargin{:}, 'split', r.split, 'verbose', 0, 'keep_on_GPU', r.keep_on_GPU);
|
||
catch err
|
||
if strcmp(err.identifier, 'parallel:gpu:array:OOM')
|
||
warning(err.message)
|
||
end
|
||
keyboard
|
||
end
|
||
|
||
rec(:,:,ind) = rec_tmp;
|
||
|
||
end
|
||
end
|
||
|
||
|
||
end
|