gtReconstructGrains.m

function varargout = gtReconstructGrains(first, last, workingdirectory, phaseID, parameters, varargin)
% GTRECONSTRUCTGRAINS  Launch a series of 3D ART (ASTRA) reconstruction on a GPU machine
%     gtReconstructGrains(first, last, workingdirectory, phaseID, [parameters], [varargin])
%     --------------------------------------------------------------------------------------------------
%
%     If submitted on a GPU machine it will work locally - otherwise it will
%     submit the job via OAR to a machine with GPU (needed for 3D ASTRA code)
%
%     INPUT:
%       first            = <int>     first grainid
%       last             = <int>     last grainid
%       workingdirectory = <string>  working directory
%       phaseID          = <int>     phase ID
%       parameters       = <struct>  {parameters.mat}
%       varargin         =           parameters.rec
%           example: 'list', [1, 3, 4:10]
%           'display', true
%     Version 001 

    if (~gtCheckGpu())
        error([mfilename 'no_GPU'], ...
            'You should submit this job on a GPU machine ! Exiting...');
    end

    currentDir = pwd;
    cd(workingdirectory);

    if (isdeployed)
        global GT_DB %#ok<NUSED,TLEV>
        global GT_MATLAB_HOME %#ok<NUSED,TLEV>
        load('workspaceGlobal.mat');
        first   = str2double(first);
        last    = str2double(last);
        phaseID = str2double(phaseID);
    end

    if (~exist('parameters','var') || isempty(parameters))
        parameters = gtLoadParameters();
    end

    if (isfield(parameters.rec, 'grains'))
        rec_gr = parameters.rec.grains;
    else
        % Old style parameters
        warning('Old "rec" parameters')
        rec_gr = parameters.rec;
        rec_gr.algorithm = 'SIRT';
    end
    rec_gr.display = false;
    if (~isempty(varargin))
        rec_gr = parse_pv_pairs(rec_gr, varargin);
        parameters.rec.grains = rec_gr;
    end

    samplefile = fullfile(parameters.acq.dir, '4_grains', 'sample.mat');
    sample = GtSample.loadFromFile(samplefile);

    if (phaseID < 1 || phaseID > numel(sample.phases))
        error([mfilename ':wrong_phase_number'], ...
            'Number of phases should be between [1, %d], but got %d', ...
            numel(sample.phases), phaseID)
    end
    num_grains = sample.phases{phaseID}.getNumberOfGrains();

    if isempty(rec_gr.list)
        grains_list = first : last;
    else
        grains_list = rec_gr.list;
        % in case the code is deployed, we have to convert to strings
        if ischar(grains_list)
            grains_list = sscanf(grains_list, '%d');
        end
    end

    % Filtering out deselected/deactivated/bad grains
    active_grains = sample.phases{phaseID}.getSelected(grains_list);

    grains_list = grains_list(active_grains);
    t = GtThreshold();
    
    fprintf('Reconstructing grains: ');
    c = tic();
    for ii = 1 : numel(grains_list)
        gr_id = grains_list(ii);
        num_chars = fprintf('%d/%d (gr: %d)', ii, numel(grains_list), gr_id);

        if (gr_id < 1 || gr_id > num_grains)
            error([mfilename ':wrong_grain_number'], ...
                'Number of grains for phase %d should be between [1, %d], but got %d', ...
                phaseID, num_grains, gr_id)
        end

        switch(upper(rec_gr.algorithm))
            case {'SIRT', '3DTV'}
                rec_vol = gtAstraReconstructGrain(gr_id, phaseID, parameters);
            case {'6DL1', '6DLS', '6DTV', '6DTVL1'}
                rec_vol = gtReconstructGrainOrientation(gr_id, phaseID, parameters);
        end
        if (~isdeployed)
            try
                t.singleGrainAutoThreshold(phaseID, gr_id);
                fprintf(repmat('\b', [1 num_chars]))
            catch
                sprintf('Segmentation of grain %d failed - check reconstruction!', gr_id)
            end
            if rec_gr.display
                GtVolView(rec_vol.intensity)
            end
        end
    end
    fprintf('Done (in %f seconds).\n', toc(c))
    if nargout > 0
        varargout{1} = rec_vol;
    end
    cd(currentDir);
end