mirror of
https://github.com/c-sooyoung/fold_slice.git
synced 2026-09-17 19:39:08 +09:00
627 lines
25 KiB
Matlab
627 lines
25 KiB
Matlab
% ATX_SUP_PARTIAL distributed (multiGPU) backprojector that allows to split the full volume into smaller pieces
|
||
% this allows to solve datasets much larger than memory of used GPU or
|
||
% spread calculations over several GPUs
|
||
%
|
||
% volData = Atx_sup_partial(projData, cfg, vectors, split, varargin)
|
||
%
|
||
% Inputs:
|
||
% **projData - array Nlayers x width_sinogram x Nangles of back-projected data
|
||
% **cfg - config structure generated by ASTRA_initialize
|
||
% **vectors - orientation of projections generated by ASTRA_initialize
|
||
% **split 3 or 4 elements vector, [split X, split Y, split Z, split angle ]
|
||
% **deformation_fields: 3x1 cell contaning 3D arrays of local deformation of the object
|
||
% **verbose - verbose <= 0 : quiet, verbose : standard info , verbose = 2: debug
|
||
% **use_shared_memory - true - share data between processed by shared memory, false = use matlab parfor distribution
|
||
% **max_memory_blocks - maximal size of used share memory memory
|
||
% **varargin - for additional parameters see the code and als the astra.Axt_partial function
|
||
% *returns*
|
||
% ++volData - backprojected volume
|
||
%
|
||
% recompile commands
|
||
% (Linux, GCC 4.8.5) mexcuda -outdir private +astra/ASTRA_GPU_wrapper/ASTRA_GPU_wrapper.cu +astra/ASTRA_GPU_wrapper/util3d.cu +astra/ASTRA_GPU_wrapper/par3d_fp.cu +astra/ASTRA_GPU_wrapper/par3d_bp.cu
|
||
% (Windows) mexcuda -outdir private ASTRA_GPU_wrapper\ASTRA_GPU_wrapper.cu ASTRA_GPU_wrapper\util3d.cu ASTRA_GPU_wrapper\par3d_fp.cu ASTRA_GPU_wrapper\par3d_bp.cu
|
||
|
||
|
||
%*-----------------------------------------------------------------------*
|
||
%| |
|
||
%| 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 volData = Atx_sup_partial(projData, cfg, vectors, split, varargin)
|
||
|
||
|
||
import utils.*
|
||
import math.*
|
||
|
||
par = inputParser;
|
||
par.addOptional('deformation_fields', {}) % deformation_fields: 3x1 cell contaning 3D arrays of local deformation of the object
|
||
par.addOptional('GPU', []) % list of GPUs to be used in reconstruction
|
||
par.addOptional('split_sub', 1) % splitting of the sub block on smaller tasks in the Atx_partial method , 1 == no splitting , otherwise [split_x,split_y,split_z,split_angles]
|
||
par.addOptional('verbose', 1) % verbose = 0 : quiet, verbose : standard info , verbose = 2: debug
|
||
par.addOptional('use_shared_memory', []) % true - share data between processed by shared memory, false = use matlab parfor distribution
|
||
par.addOptional('max_memory_blocks', min(50e9, utils.check_available_memory*1e6/4)) % maximal size of used share memory memory
|
||
|
||
|
||
par.parse(varargin{:})
|
||
r = par.Results;
|
||
|
||
if isscalar(split)
|
||
split = split .* ones(1,3);
|
||
end
|
||
|
||
split_orig = split;
|
||
|
||
if isempty(r.use_shared_memory)
|
||
r.use_shared_memory = length(r.GPU) > 1;
|
||
end
|
||
if isscalar(r.split_sub)
|
||
r.split_sub = r.split_sub .* ones(1,3);
|
||
end
|
||
|
||
gpu = gpuDevice;
|
||
if isempty(r.GPU)
|
||
r.GPU = gpu.Index;
|
||
end
|
||
N_GPU = length(r.GPU);
|
||
|
||
if ~isempty(r.deformation_fields)
|
||
% deformation field splitting not implemneted for split_sub;
|
||
split(1:3) = split(1:3) .* r.split_sub(1:3);
|
||
r.split_sub = 1;
|
||
end
|
||
%% if not splitting on this level is requirested, continue to tomo.Atx_partial
|
||
if all(split(1:3) == 1) && N_GPU == 1 && gpu.AvailableMemory > 4*(numel(projData)/10+cfg.iVolX*cfg.iVolY*cfg.iVolZ)
|
||
volData = astra.Atx_partial(projData, cfg, vectors, r.split_sub, ...
|
||
'GPU', r.GPU, 'deformation_fields', r.deformation_fields, 'verbose', r.verbose);
|
||
return
|
||
elseif all(split(1:2) == 1) && isempty(r.deformation_fields) && ...
|
||
utils.check_available_memory*1e6 > 4*(cfg.iVolX*cfg.iVolY*cfg.iVolZ*N_GPU+numel(projData)) && ...
|
||
cfg.iVolX*cfg.iVolY*cfg.iVolZ*4 < min(4*double(intmax('int32')),gpu.TotalMemory/2)
|
||
|
||
%% if no splitting is needed and the volume is small enough then at least split the data on multiple GPUs by angles
|
||
volData = Atx_angle_split(projData, cfg, vectors, r.split_sub, r.GPU, r.verbose);
|
||
return
|
||
end
|
||
|
||
|
||
%% otherwise do proper checking of all inputs before splitting
|
||
if ~(isa(projData, 'gpuArray') && strcmp(classUnderlying(projData), 'single')) && ...
|
||
~isa(projData, 'single')
|
||
error('Only single precision input array supported')
|
||
end
|
||
projData = gather(projData);
|
||
|
||
if any(cfg.pixel_scale ~= 1)
|
||
error('Non integer pixel size is not implemented, try to use Atx_partial')
|
||
end
|
||
if cfg.skewness_angle ~= 0
|
||
error('Skewness is not working well with sup-split')
|
||
end
|
||
if length(cfg.pixel_scale) == 2 && cfg.pixel_scale(1) ~= cfg.pixel_scale(2)
|
||
error('Variable pixel size for each axis is not implemented')
|
||
end
|
||
Nvol_orig = [cfg.iVolX,cfg.iVolY,cfg.iVolZ];
|
||
|
||
|
||
split(3) = max(split(3), ceil(split(3) * ( 4*(numel(projData)/5+cfg.iVolX*cfg.iVolY*cfg.iVolZ*(1+any(r.split_sub>1)) ) / gpu.AvailableMemory / prod(split) )));
|
||
split(3) = max(split(3), ceil(split(3) * ( ( cfg.iVolX*cfg.iVolY*cfg.iVolZ)/prod(split(1:3)) / double(intmax('int32')) )));
|
||
% make it equaly splitable among the GPUs
|
||
split(3) = ceil(max(split(3), ceil(prod(split(1:3)) / N_GPU) * N_GPU) / prod(split(1:2)));
|
||
split(1:3) = min(split(1:3), Nvol_orig(1:3));
|
||
|
||
% adjust splitting to make get equal blocks of the volume
|
||
for i = 1:100
|
||
split(1:3) = ceil(Nvol_orig ./ floor(Nvol_orig ./ split(1:3)));
|
||
end
|
||
if prod(r.split_sub) > 1
|
||
% reduce the split_sub if possible
|
||
r.split_sub(3) = ceil(r.split_sub(3) / 2^nextpow2(prod(split) / prod(split_orig)) );
|
||
end
|
||
if length(split) < 4
|
||
split(4) = 1;
|
||
end
|
||
|
||
% initial parameters check
|
||
Nvol_sub = Nvol_orig./split(1:3);
|
||
Nproj_full = size(projData);
|
||
|
||
assert(all(mod(Nvol_sub,1)==0), sprintf('Volume array cannot be divided to %i %i %i cubes', split(1:3)))
|
||
assert(all(Nproj_full==[cfg.iProjV,cfg.iProjU,cfg.iProjAngles]), 'Wrong inputs size')
|
||
assert(all(size(vectors)==[cfg.iProjAngles,12]), 'Wrong vectors size')
|
||
|
||
% size of the subprojection of single subvolume
|
||
Nproj_sub = [ (Nvol_sub(3)* sind(cfg.lamino_angle) + sqrt(sum(Nvol_sub(1:2).^2))*cosd(cfg.lamino_angle)), ...
|
||
sqrt(sum(Nvol_sub(1:2).^2))];
|
||
|
||
% adjust sub projection size to account for inplane rotation of the geometry
|
||
if cfg.tilt_angle ~= 0
|
||
Nproj_rot = [cosd(cfg.tilt_angle), -sind(cfg.tilt_angle); +sind(cfg.tilt_angle), cosd(cfg.tilt_angle)] * [0,0; Nproj_sub(1:2)];
|
||
% calculate projection window size after rotation
|
||
Nproj_rot = max(Nproj_rot) - min(Nproj_rot);
|
||
% add some extra padding
|
||
Nproj_sub(1:2) = 2*(Nproj_rot - Nproj_sub(1:2)) + Nproj_sub(1:2);
|
||
end
|
||
|
||
% provide extra space for subpixel (linear) interpolation at the borders of the split volumes, needed only for noninteger CoR_offset
|
||
Nproj_sub = Nproj_sub + 2*([split(3), max(split([1,2]))] - 1);
|
||
|
||
|
||
Nproj_sub = ceil(Nproj_sub/16)*16; % make it easier splitable for ASTRA
|
||
if cfg.lamino_angle == 90
|
||
Nproj_sub(1) = min(Nproj_sub(1), Nproj_full(1));
|
||
end
|
||
% only if there is not split in the horizontal dimension
|
||
if all(split(1:2) == 1)
|
||
% do not take larger than size of the inputs
|
||
Nproj_sub = min(Nproj_sub, Nproj_full(1:2));
|
||
end
|
||
|
||
% avoid low RAM issues due to paralelization
|
||
split(4) = max(split(4), ceil(4*prod(Nproj_sub)*cfg.iProjAngles/split(4)/(r.max_memory_blocks/N_GPU)));
|
||
|
||
Nproj_sub(3) = ceil(cfg.iProjAngles/split(4));
|
||
|
||
|
||
cfg_small = cfg;
|
||
cfg_small.iVolX = Nvol_sub(1);
|
||
cfg_small.iVolY = Nvol_sub(2);
|
||
cfg_small.iVolZ = Nvol_sub(3);
|
||
% get new size of projections
|
||
cfg_small.iProjU = Nproj_sub(2);
|
||
cfg_small.iProjV = Nproj_sub(1);
|
||
|
||
% calculate and store offset of the center of rotation, it will be used later
|
||
offset = vectors(:,4:6) +(vectors(:,10:12).*cfg.iProjV/2+vectors(:,7:9).*cfg.iProjU/2 );
|
||
CoR_offset = -[dot(offset', vectors(:,10:12)') ./ dot(vectors(:,10:12)', vectors(:,10:12)');
|
||
dot(offset', vectors(:,7:9)') ./ dot(vectors(:,7:9)', vectors(:,7:9)')]' ;
|
||
|
||
% remove centering offset && apply new offset
|
||
shift_vec = vectors(:,10:12)*(cfg.iProjV/2-cfg_small.iProjV/2)+vectors(:,7:9)*(cfg.iProjU/2-cfg_small.iProjU/2);
|
||
vectors(:,4:6) = vectors(:,4:6) + shift_vec;
|
||
|
||
%% keep the volume in RAM , transfer only sub-blocks
|
||
Nblocks = prod(split);
|
||
%% write back to the preallocated shared array
|
||
volData = zeros(Nvol_orig, 'single');
|
||
|
||
if N_GPU > 1
|
||
poolobj = gcp('nocreate');
|
||
if isempty(poolobj) || poolobj.NumWorkers < N_GPU
|
||
delete(poolobj);
|
||
poolobj = parpool(N_GPU);
|
||
end
|
||
poolobj.IdleTimeout = 600; % set idle timeout to 10 hours
|
||
end
|
||
|
||
% when the function is finished, make sure to execute following code
|
||
global status
|
||
status = true;
|
||
if r.use_shared_memory
|
||
out = onCleanup(@()myCleanupFun());
|
||
end
|
||
% run blocks in series
|
||
% run sub-blocks on each GPU in parallel
|
||
t_total = tic();
|
||
clear output
|
||
|
||
%% START OF OUTER GPU LOOP
|
||
outputs_blocks = [];
|
||
%% unitialize one solver per GPU
|
||
for thread_id = 1:N_GPU
|
||
% parse inputs and try to split them if possible
|
||
[outputs_blocks, inputs_block{thread_id},cfg_all{thread_id}] = ...
|
||
submit_block(thread_id, thread_id, outputs_blocks, projData, cfg, cfg_small, vectors,CoR_offset, split, Nproj_sub,Nvol_sub, r, varargin{:} );
|
||
end
|
||
|
||
unprocessed_blocks = N_GPU+1:Nblocks;
|
||
|
||
%% merge blocks back from GPUs and write to the shared array volData
|
||
for ii = 1:Nblocks
|
||
if r.verbose; utils.progressbar(ii, Nblocks); end
|
||
% set values from the small blocks to the final output arrays
|
||
[thread_id, timing, id] = gather_block( outputs_blocks,volData, cfg_all, split(4) > 1);
|
||
|
||
if ~isempty(unprocessed_blocks)
|
||
block_id = unprocessed_blocks(1);
|
||
unprocessed_blocks(1) = [];
|
||
|
||
if isa(outputs_blocks, 'parallel.FevalFuture') && sum([outputs_blocks.Read]) ~= 1
|
||
outputs_blocks
|
||
keyboard
|
||
end
|
||
|
||
% submit a new job once the previous is finished
|
||
[outputs_blocks, inputs_block{thread_id},cfg_all{block_id}] = ...
|
||
submit_block(block_id, thread_id, outputs_blocks, projData, cfg, cfg_small, vectors, CoR_offset, split, Nproj_sub,Nvol_sub, r, varargin{:} );
|
||
|
||
if isa(outputs_blocks, 'parallel.FevalFuture') && any(cat(1,[outputs_blocks.Read]))
|
||
outputs_blocks
|
||
keyboard
|
||
end
|
||
end
|
||
end
|
||
|
||
if r.verbose > 1
|
||
fprintf('Timing system: GPU init %3.2gs shared_mem down %3.2gs upload on GPU %3.2gs tomo projection %3.2gs download from GPU %3.2gs shared_mem %3.2gs \n',sum(timing,2) )
|
||
if length(r.GPU) > 1
|
||
fprintf('Timing local: GPU init %3.2gs shared_mem down %3.2gs upload on GPU %3.2gs tomo projection %3.2gs download from GPU %3.2gs shared_mem %3.2gs \n ',sum(timing,2)/max(1,length(r.GPU)) )
|
||
fprintf('Total time %3.2fs, parfor overhead %3.2fs \n', t_total, t_total - sum(sum(timing,2)/max(1,length(r.GPU))) )
|
||
end
|
||
end
|
||
|
||
|
||
|
||
% finish GPU blocks
|
||
% everything was fine -> no cleaning needed
|
||
status = false;
|
||
|
||
end
|
||
|
||
function [outputs_blocks,inputs_block, cfg_out] = submit_block(block_id, thread_id, outputs_blocks, projData, cfg, cfg_small, vectors, CoR_offset,split, Nproj_sub,Nvol_sub, r, varargin )
|
||
% prepare blocks for asynchonous processing
|
||
try
|
||
inputs_block = prepare_block(block_id, projData, cfg, cfg_small, vectors,CoR_offset, split, Nproj_sub,Nvol_sub, r, varargin{:});
|
||
catch err
|
||
disp(getReport(err))
|
||
keyboard
|
||
end
|
||
N_GPU = length(r.GPU);
|
||
if isempty(outputs_blocks); clear outputs_blocks; end
|
||
|
||
cfg_out = inputs_block{2};
|
||
|
||
|
||
try
|
||
%% process preloaded data
|
||
% no paralel toolbox
|
||
if N_GPU <= 1
|
||
[outputs_blocks{thread_id}.volData_small,outputs_blocks{thread_id}.timing, outputs_blocks{thread_id}.id]=...
|
||
run_partial_projector(inputs_block, block_id, 1,r.GPU,0);
|
||
else
|
||
% run it asynchronously
|
||
if r.verbose > 3
|
||
ticBytes(gcp);
|
||
end
|
||
outputs_blocks(thread_id) = parfeval(@run_partial_projector, 3, inputs_block,block_id,thread_id, r.GPU,0);
|
||
if r.verbose > 3
|
||
try; tocBytes(gcp); end
|
||
end
|
||
end
|
||
catch err
|
||
disp(getReport(err))
|
||
utils.check_available_memory
|
||
keyboard
|
||
end
|
||
end
|
||
|
||
function myCleanupFun()
|
||
% destroy all shared memory that could have been left behind
|
||
!ipcs -m | cut -d' ' -f2 | grep '^[0-9]' | while read x; do ipcrm -m $x; done
|
||
end
|
||
|
||
function [thread_id, timing, id] = gather_block(output_package,volData, cfg_all, add_values )
|
||
|
||
if isa(output_package, 'parallel.FevalFuture')
|
||
% gather results from cluster , WAIT FOR CALCULATIONS TO BE FINISHED
|
||
% [~, outputs_block, id] = fetchNext(output_package);
|
||
%% my version of the fetchNext function, it seems faster
|
||
id = [];
|
||
assert(any(~[output_package.Read]), 'All blocks are already read')
|
||
while true
|
||
for thread_id =1:length(output_package)
|
||
if strcmpi(output_package(thread_id).State, 'finished') && output_package(thread_id).Read == 0
|
||
try
|
||
[volData_small,timing,id] = output_package(thread_id).fetchOutputs;
|
||
catch err
|
||
if strcmpi(err.identifier, 'parallel:fevalqueue:InvalidExecutionResult')
|
||
warning('Unknown error, trying to restart parpool')
|
||
delete(gcp('nocreate'));
|
||
end
|
||
if strcmpi(err.identifier, 'parallel:fevalqueue:InvalidExecutionResult')
|
||
warning('Unknown error, trying to restart parpool')
|
||
delete(gcp('nocreate'));
|
||
end
|
||
if ~isempty(output_package(thread_id).Diary)
|
||
fprintf('============ THREAD %i FAILED, OUTPUT: ============= \n', thread_id)
|
||
disp(output_package(thread_id).Diary)
|
||
end
|
||
fprintf('============ THREAD %i FAILED, ERROR: ============= \n', thread_id)
|
||
disp(getReport(output_package(thread_id).Error))
|
||
|
||
keyboard
|
||
|
||
output_package.cancel
|
||
rethrow(err)
|
||
end
|
||
break
|
||
end
|
||
end
|
||
if ~isempty(id); break; end
|
||
pause(0.01) % wait for the data to be prepared
|
||
end
|
||
elseif iscell(output_package)
|
||
thread_id = 1;
|
||
id = output_package{thread_id}.id;
|
||
volData_small = output_package{thread_id}.volData_small;
|
||
timing = output_package{thread_id}.timing;
|
||
else
|
||
disp('FAILED ?? ')
|
||
keyboard
|
||
end
|
||
|
||
if isempty(volData_small)
|
||
warning('ASTRA projection probably failed')
|
||
keyboard
|
||
end
|
||
if isa(volData_small, 'shm')
|
||
% load data from shared memory
|
||
[s,volData_small] = volData_small.attach;
|
||
s.protected = false; % release shared memory
|
||
elseif ~isnumeric(volData_small)
|
||
keyboard
|
||
end
|
||
|
||
|
||
% write back to the full array stored in RAM
|
||
positions = zeros(size(volData_small,3),2)+cfg_all{id}.volume_shift(1:2);
|
||
indices = (1:size(volData_small,3)) + cfg_all{id}.volume_shift(3);
|
||
% use a MEX code to speed it up
|
||
utils.add_to_3D_projection(volData_small, volData,positions,indices, true, false);
|
||
|
||
|
||
end
|
||
|
||
|
||
|
||
function prepared_block = prepare_block(id, projData, cfg, cfg_small, vectors,CoR_offset, split, Nproj_sub, Nvol_sub ,r, varargin)
|
||
|
||
[x,y,z,angle_block_id] = ind2sub(split,id);
|
||
block_pos = [x,y,z];
|
||
|
||
shift = zeros(3,1);
|
||
for n = 1:3
|
||
%% find optimal shift of the subvolume
|
||
if mod(split(n),2)==1 %% odd
|
||
shift(n) = (block_pos(n) - ceil(split(n)/2))*Nvol_sub(n);
|
||
else
|
||
shift(n) = (block_pos(n) - split(n)/2-1/2)*Nvol_sub(n);
|
||
end
|
||
end
|
||
Nangle_per_blocks = ceil(cfg.iProjAngles / split(4));
|
||
|
||
%% for splitting to angular blocks using close angles, it makes ASTRA faster
|
||
angle_ids = (1+(angle_block_id-1)*Nangle_per_blocks:min(cfg.iProjAngles, angle_block_id*Nangle_per_blocks));
|
||
cfg_small.iProjAngles = length(angle_ids);
|
||
Nproj_sub(3) = cfg_small.iProjAngles;
|
||
|
||
|
||
%% find optimal shift of the subvolume
|
||
deform_fields_small = get_subdeform_fields(block_pos, r.deformation_fields, split);
|
||
|
||
CoR_offset = [CoR_offset, zeros(cfg.iProjAngles,1)] ;
|
||
|
||
% shift the sub-projections off center to create single
|
||
% large projection after assembling, if shift == 0 =>
|
||
% projections will be rotationally centered
|
||
vec = vectors;
|
||
|
||
% apply optimal shift
|
||
projection_shift = cfg.pixel_scale(1)^2.*[vec(:,10:12)*shift, vec(:,7:9)*shift, zeros(cfg.iProjAngles,1)];
|
||
projection_shift = bsxfun(@plus, projection_shift, [cfg.iProjV/2-cfg_small.iProjV/2,cfg.iProjU/2-cfg_small.iProjU/2,0] + CoR_offset);
|
||
|
||
% calculate subpixel shifts
|
||
projection_shift_subpix = projection_shift - round(projection_shift) - CoR_offset;
|
||
projection_shift = round(projection_shift);
|
||
|
||
% apply subpixel shifts
|
||
vec(:,4:6) = vec(:,4:6) - ...
|
||
( bsxfun(@times,vectors(:,7:9),projection_shift_subpix(:,2))+ ...
|
||
bsxfun(@times,vectors(:,10:12),projection_shift_subpix(:,1)));
|
||
% just store shifts for later
|
||
cfg_small.volume_shift = (block_pos-1).*Nvol_sub;
|
||
cfg_small.projection_shift = projection_shift;
|
||
|
||
% move data using custom made MEX routine
|
||
% allocate a small sub array
|
||
projData_small = zeros(Nproj_sub, 'single');
|
||
|
||
if r.use_shared_memory
|
||
s = shm();
|
||
try
|
||
s.allocate(projData_small) % attach the shared memory
|
||
catch
|
||
keyboard
|
||
end
|
||
[s, projsmall_shm] = s.attach();
|
||
% === write data =====
|
||
% use self-made MEX OMP function to move the data s
|
||
utils.get_from_3D_projection(projsmall_shm, projData, projection_shift(angle_ids,1:2), angle_ids);
|
||
% detach the shared memory
|
||
projData_small = s;
|
||
s.detach;
|
||
else
|
||
% use custom made MEX OMP function to move the data
|
||
utils.get_from_3D_projection(projData_small, projData, projection_shift(angle_ids,1:2), angle_ids);
|
||
end
|
||
|
||
% return only vector for the use angles
|
||
vec = vec(angle_ids,:);
|
||
|
||
prepared_block = {projData_small, cfg_small, vec, r.split_sub, varargin{:}, ...
|
||
'deformation_fields', deform_fields_small, 'verbose', 0, 'GPU', []};
|
||
|
||
end
|
||
|
||
function [vol, timing,block_id] = run_partial_projector(prepared_block, block_id, thread_id, GPU_list, verbose)
|
||
|
||
t0 = tic;
|
||
gpu_id = GPU_list(thread_id);
|
||
|
||
% let matlab to choose which GPU use
|
||
gpu = gpuDevice();
|
||
if ~isempty(GPU_list) && gpu.Index ~= gpu_id
|
||
gpuDevice(gpu_id); % avoid unneeded initalization
|
||
end
|
||
t_init = toc(t0);
|
||
|
||
timing = [t_init, 0,0,0,0,0];
|
||
t = tic;
|
||
if isa(prepared_block{1}, 'shm')
|
||
% data are attached to shared memory
|
||
[s,projData_small] = prepared_block{1}.attach();
|
||
else
|
||
% data are given directly to the worker
|
||
projData_small = prepared_block{1};
|
||
end
|
||
timing(2) = toc(t);
|
||
t = tic;
|
||
timing(3) = toc(t);
|
||
% call the next level abstraction
|
||
t = tic;
|
||
is_remote = ~isempty(getCurrentTask());
|
||
|
||
|
||
% call the next level abstraction around ASTRA wrapper
|
||
try
|
||
for ii = 1:2
|
||
try
|
||
vol = astra.Atx_partial(projData_small, prepared_block{2:end},'verbose', max(is_remote, verbose));
|
||
break
|
||
catch err
|
||
if ii == 2
|
||
rethrow(err)
|
||
else
|
||
warning('Atx_partial failed, trying again ... ')
|
||
end
|
||
end
|
||
end
|
||
|
||
catch err
|
||
qkeyboard
|
||
gpu = gpuDevice
|
||
reset(gpu)
|
||
fprintf('Error on GPU %i / %i', gpu.Index, gpuDeviceCount)
|
||
disp( getReport(err, 'extended', 'hyperlinks', 'on'))
|
||
% projData_small = single([]);
|
||
rethrow(err)
|
||
end
|
||
%%%%%%%%%%%
|
||
timing(4) = toc(t);
|
||
t = tic;
|
||
vol = gather(vol); % move to RAM
|
||
timing(5) = toc(t);
|
||
|
||
t = tic;
|
||
if isscalar(prepared_block{1})
|
||
% data are distributed to shared memory
|
||
s.detach();
|
||
s = shm(true);
|
||
s.upload(vol);
|
||
vol = s;
|
||
end
|
||
timing(6) = toc(t);
|
||
end
|
||
|
||
|
||
function deformation_fields_sub = get_subdeform_fields(block_pos, deform_fields, split)
|
||
% crop the deformation field only int the region of interest
|
||
if isempty(deform_fields)
|
||
deformation_fields_sub = {};
|
||
return
|
||
end
|
||
if prod(split)==1
|
||
deformation_fields_sub = deform_fields;
|
||
return
|
||
end
|
||
|
||
|
||
for ii = 1:3
|
||
N_full = size(deform_fields{1,ii});
|
||
N_small = ceil( N_full ./ reshape(split(1:3),[],1)');
|
||
for kk = 1:3
|
||
ind_def{kk} = (1+(block_pos(kk)-1)*N_small(kk)) : min(N_full(kk), (split(kk))*N_small(kk));
|
||
end
|
||
for jj = 1:2
|
||
deformation_fields_sub{jj,ii} = deform_fields{jj,ii}(ind_def{:});
|
||
end
|
||
end
|
||
|
||
warning('Splitting of deformation field is not supported / recommended')
|
||
|
||
|
||
end
|
||
|
||
|
||
function volData = Atx_angle_split(projData, cfg, vectors, split, GPU, verbose)
|
||
%% simplified version for multiGPU reconstruction of volumes that are small enough
|
||
% simply split the projections by angles and apply each projection
|
||
% block on a separated GPU
|
||
Nangles = cfg.iProjAngles;
|
||
N_GPU = length(GPU);
|
||
Nblocks = ceil(Nangles / N_GPU);
|
||
if length(split) == 4
|
||
split(4) = max(1, split(4) / N_GPU);
|
||
end
|
||
if verbose; utils.progressbar(1,4); end
|
||
for id = 1:N_GPU
|
||
ind{id} = 1+(id-1)*Nblocks:min(id*Nblocks,Nangles);
|
||
projData_shm{id} = shm();
|
||
projData_shm{id} = tomo.get_from_array(projData, projData_shm{id}, ind{id});
|
||
volData_shm{id} = [];
|
||
end
|
||
if verbose; utils.progressbar(2,4); end
|
||
parfor(id = 1:N_GPU, N_GPU)
|
||
t = getCurrentTask();
|
||
gpu_id = GPU(1+mod(t.ID-1, N_GPU));
|
||
% let parfor to choose which GPU use
|
||
gpu = gpuDevice();
|
||
if ~isempty(GPU) && gpu.Index ~= gpu_id
|
||
gpuDevice(gpu_id); % avoid unneeded initalization
|
||
end
|
||
[s,projData_block] = projData_shm{id}.attach();
|
||
volData_block = astra.Atx_partial(projData_block, cfg, vectors(ind{id},:),split, 'keep_on_GPU', true, 'verbose',0);
|
||
volData_block = gather(volData_block);
|
||
s.upload(volData_block);
|
||
s.protected = true;
|
||
volData_shm{id} = s;
|
||
end
|
||
if verbose; utils.progressbar(3,4); end
|
||
volData = single(0);
|
||
for id = 1:N_GPU
|
||
[s,volData_block] = volData_shm{id}.attach();
|
||
volData = volData + volData_block;
|
||
s.free;
|
||
end
|
||
if verbose; utils.progressbar(4,4); end
|
||
end
|