Skip to content
Snippets Groups Projects
Select Git revision
  • charles_rmn
  • new_diffractometer
  • Some-additional-functions-from-Zheheng
  • 11-add-calculation-of-2DIGM-and-IGM-components
  • 2-add-TRD-fwd-proj
  • master default protected
  • 10-create-twin-cluster-fwd-proj
  • 1-generic-data-loading-layer
  • 7-notebook-preprocessing-and-segmentation
  • 6-py-bkg-rm
  • 8-virtualenv
  • 7-pytest
  • 5-git-commit-hooks
  • 5-git-commit-hooks-2
  • patch-1
  • dev/python-bkg-rm
16 results
Blame Permalink
gtPreprocessing.m 19.65 KiB 
function gtPreprocessing()
% GTPREPROCESSING  Handles the preprocessing of images
%     gtPreprocessing()
%     -----------------
%
%     Creates median refs (refHST)
%     Creates abs/ext/full images
%     Corrects drifts and rotation axis to centre of full/ext/abs images
%     Runs from the top level of the analysis directory (where the parameters.mat file is)
%     gtSequenceMedianRefs uses parameters.mat rather than .xml for info
%
%     Version 003 30-09-2013 by LNervo
%       Created gtPrepDefaultParameters to clean a bit the function
%
%     Version 002 08-11-2011 by LNervo
%       *** Cleanest version (first attempt) ***
%       Why double calculation or parameters settings??
%
%     Version 001 29-04-2011 by AKing
%       creation of dark.edf image
%       calculation of sample bounding box and centre of rotation
%       calculation of sample drifts via a Quali scan


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% recompile binaries - not often needed.  But more often than it should be...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

check = inputwdefault('Recompile functions for OAR? Not normally needed! [y/n]', 'n');
if strcmpi(check,'y')
    disp('recompiling all functions')
    funcsToCompile = { ...
        'gtSequenceMedianRefs', ...
        'gtCreateAbsLive', ...
        'gtAbsMedianLive', ...
        'gtMovingMedianLive', ...
        'gtCreateFullLive' };
    gtExternalCompileFunctions(funcsToCompile{:});
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loading parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
parameters = [];
list       = [];
if exist('parameters.mat', 'file')
    load('parameters.mat');
else
    disp('Parameters file not found - cd to the analysis directory')
    disp('If you have not already run gtSetup, you should do this first')
    gtError('parameters.mat not found!')
end
if exist('list.mat','file')
    load('list.mat');
else
    [~,list]=make_parameters(2);
end

parameters_name = fullfile(parameters.acq.dir,'parameters.mat');

% now prep contains totproj
prep = gtPrepDefaultParameters(parameters.acq);
check = inputwdefault('Do you want to reset all the parameters for preprocessing? [y/n]', 'n');
if strcmpi(check,'y')
    parameters.prep = gtAddMatFile(parameters.prep, prep, true, false, false);
    save(parameters_name,'parameters');
    disp('Preprocessing parameters have been reset to default values')
end

% total number of projections
totproj = prep.totproj;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% check preprocessing parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
parameters = gtCheckParameters(parameters, 'prep', 'verbose', true);
save(parameters_name,'parameters');

check = inputwdefault('Do you want to verify/change the (default) preprocessing parameters? [y/n]', 'n');
if strcmpi(check,'y')
    
    bad = true;
    while (bad)
        parameters.prep = gtModifyStructure(parameters.prep, list.prep, [1 2], 'Preprocessing parameters:');
        bad = false;
        % test values
        if mod(totproj, parameters.prep.absint) ~= 0
            bad = true;
            disp('absint should be a factor of the number of images per turn!')
            list.prep{findValueIntoCell(list.prep(:,1), 'absint'),2} = '!!! Moving median interval, direct beam (images, mod(total projection, absint)=0) !!!';
            parameters.prep.absint   = prep.absint;
        end
        if mod(totproj, parameters.prep.fullint) ~= 0
            bad = true;
            disp('fullint should be a factor of the number of images per turn!')
            list.prep{findValueIntoCell(list.prep(:,1), 'fullint'),2} = '!!! Moving median interval, diffracted image (images, mod(total projection, fullint)=0) !!!';
            parameters.prep.fullint  = prep.fullint;
        end
        if mod(parameters.prep.absrange/parameters.prep.absint, 2) ~= 0
            bad = true;
            disp('absrange should be an even multiple of absint!')
            list.prep{findValueIntoCell(list.prep(:,1), 'absrange'),2} = '!!! Moving median range, direct beam (images, n*2*absint) !!!';
            parameters.prep.absrange = 10*parameters.prep.absint;
        end
        if mod(parameters.prep.fullrange/parameters.prep.fullint, 2) ~= 0
            bad = true;
            disp('fullrange should be an even multiple of fullint!')
            list.prep{findValueIntoCell(list.prep(:,1), 'fullrange'),2} = '!!! Moving median range, diffracted image (images, n*2*fullint) !!!';
            parameters.prep.fullrange = 10*parameters.prep.fullint;
        end
    end
    save(parameters_name,'parameters');
    disp('Parameters have been updated and saved on disk')
end % end verify parameters

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set parameters file correct for start of preprocessing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% after the copying has been done, the collection_dir in parameters
% points to an empty directory. Collection_dir_old becomes a list of the collection_dir from
% which the gtSetup has been launched. Directories are separated by ':'
% Then collection_dir points to the /0_rawdata/Orig folder
parameters.acq.collection_dir_old = [parameters.acq.collection_dir ':' parameters.acq.collection_dir_old];
parameters.acq.collection_dir     = fullfile(parameters.acq.dir, '0_rawdata', 'Orig');

% reset renumbering to false
parameters.prep.renumbered = prep.renumbered;

save(parameters_name,'parameters');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create dark.edf file of scan
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
darkname    = fullfile(parameters.acq.dir,'0_rawdata',parameters.acq.name,'dark.edf');
darkendname = fullfile(parameters.acq.dir,'0_rawdata',parameters.acq.name,'darkend0000.edf');
if ~exist(darkname,'file')
    disp('Creating dark.edf image of scan.')
    [dark,info] = edf_read(darkendname);
    dark = dark/parameters.acq.ndark;
    info.datatype = 'float32';
    edf_write(dark,darkname,info);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% acq/prep parameters determined from images
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if processing "live", may not have all images required.  Use sensible
% defaults or compromises, that can be refined later, and give warnings!
% dangerous to allow this step to be skipped, as it could lead to applying
% a single rotation axis shift more than once

% special case processing for scan where there is no direct beam image
if (parameters.acq.no_direct_beam)
    disp('%%% no direct beam - manual input of bounding boxes %%')
    %similarly for bounding box - manual to allow zeros to be input for offset
    %scan, but can put a bounding box around the glow at the centre of a
    %taper image, for example
    check = inputwdefault('Do you need to define a direct beam BoundingBox graphically? [y/n]', 'n');
    if strcmpi(check,'y')
        parameters.acq.bb    = gtFindDirectBeamBB(parameters.acq);
        parameters.acq.bbdir = parameters.acq.bb;
        parameters.prep.bbox = parameters.acq.bb;
    else
        parameters.acq.bb    = inputwdefault('Insert the "sample bounding box" (without squared brackets)',num2str(parameters.acq.bb));
        parameters.acq.bb    = str2num(parameters.acq.bb);
        parameters.prep.bbox = parameters.acq.bb;
        parameters.acq.bbdir = parameters.acq.bb;
    end
    check = inputwdefault('Is it like a taper scan with the sample still in the images? [y/n]', 'y');
    if strcmpi(check,'y')
        % define the sample envelope - assume centred on direct beam
        labgeo_tmp = gtGeoSamEnvFromAcq(parameters.labgeo, parameters.acq);
        parameters.labgeo.samenvtop = labgeo_tmp.samenvtop;
        parameters.labgeo.samenvbot = labgeo_tmp.samenvbot;
        parameters.labgeo.samenvrad = labgeo_tmp.samenvrad;
    end
else
    % find direct beam bounding box
    if isempty(parameters.acq.bbdir)
        check = 'y';
    else
        check = inputwdefault(['Do you want to redefine the direct beam BoundingBox (currently ' sprintf('[%d %d %d %d]', parameters.acq.bbdir) ')? [y/n]'], 'n');
    end
    if strcmpi(check,'y')
        parameters.acq.bbdir = gtFindDirectBeamBB(parameters.acq);
    end
    save(parameters_name,'parameters');

    bbdir = parameters.acq.bbdir;

    % determine rotation axis
    % "u" is the horizontal pixel coordinate
    if isempty(parameters.acq.rotu)
        check='y';
    else
        rot = parameters.acq.rotu;
        check = inputwdefault(['Do you want to (re)-calculate the Rotation Axis position (currently ' sprintf('%0.1f',rot) ')? [y/n]'], 'n');
    end
    if strcmpi(check,'y')
        parameters.acq.rotu = gtFindRotationAxis(parameters.acq,bbdir);
        parameters.acq.rotx = parameters.acq.rotu;
        parameters.labgeo.detrefpos(2) = -(parameters.acq.rotu - parameters.acq.xdet / 2) * parameters.acq.pixelsize;
    end
    save(parameters_name,'parameters');

    %determine the limits of sample horizontally for acq.bb
    %this should be centred on the rotation axis (parameters.acq.rotu)
    if isempty(parameters.acq.bb)
        check = 'y';
    else
        check = inputwdefault(['Do you want to (re)-define the Sample BoundingBox (currently ' sprintf('[%d %d %d %d]',parameters.acq.bb) ')? [y/n]'], 'n');
    end
    if strcmpi(check,'y')
        parameters.acq.bb  = gtFindSampleEdges(parameters.acq,parameters.prep,bbdir);
        parameters.prep.bbox = parameters.acq.bb+[-parameters.prep.margin 0 2*parameters.prep.margin 0];
        
        % define the sample envelope - assume centred on direct beam
        labgeo_tmp = gtGeoSamEnvFromAcq(parameters.labgeo, parameters.acq);
        parameters.labgeo.samenvtop = labgeo_tmp.samenvtop;
        parameters.labgeo.samenvbot = labgeo_tmp.samenvbot;
        parameters.labgeo.samenvrad = labgeo_tmp.samenvrad;        
    end
    save(parameters_name,'parameters');
end % end no_direct_beam case

if isempty(parameters.acq.maxradius)
    check='y';
else
    check = inputwdefault('Do you want to (re)-define the Active Area of the detector? [y/n]', 'n');
end
if strcmpi(check,'y')
    parameters.acq.maxradius = gtFindActiveArea(parameters);
end

% Check labgeo parameters
header = 'Check setup geometry parameters';
parameters.labgeo = gtModifyStructure(parameters.labgeo, list.labgeo, 1, header);

save(parameters_name,'parameters');
disp(' ')
disp('Parameters have been saved in file:')
disp(parameters_name)
disp(' ')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate sample drifts
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% this is the time to analyse a calibration or quali scan to determine
% sample drifts
ydrift=zeros(1, totproj);  % default values for drift
zdrift=zeros(1, totproj);

check = inputwdefault('Do you want to perform sample drift analysis for this measurement? [y/n]', 'n');
if strcmpi(check,'y')
    disp('Quali scan analysis')
    showcorr = inputwdefault('Show image correlations? [y/n]', 'y');
    if strcmpi(showcorr, 'y')
        showcorr = true;
    else
        showcorr = false;
    end
    if isfield(parameters, 'calib') && isfield(parameters.calib, 'zdrift') && isfield(parameters.calib, 'ydrift')
        disp('parameters file already has calibration scan data')
        check = inputwdefault('Repeat calibration scan analysis? [y/n]', 'y');
        if strcmpi(check, 'y')
            disp('calling gtSetupCalibrationQuali.m')
            warning('gtSetupCalibrationQuali:matlab:v1','* caution * not sure if coordinate system changes or non-centred direct beam bb could affect this')
            parameters = gtSetupCalibrationQuali(parameters.acq.bb,showcorr);

            % output is automatically saved into parameters file
            % as parameters.calib.ydrift and parameters.calib.zdrift
            % modified function to also return updated parameters
        end
    else
        disp('calling gtSetupCalibrationQuali.m')
        warning('gtSetupCalibrationQuali:matlab:v2','* caution * not sure if coordinate system changes or non-centred direct beam bb could affect this')
        parameters = gtSetupCalibrationQuali(parameters.acq.bb,showcorr);
    end
    ydrift=parameters.calib.ydrift;
    zdrift=parameters.calib.zdrift;
end % sample drift analysis

% shifts will only be applied after flatfield correction of full, abs and
% ext images !
if (~parameters.acq.no_direct_beam)
    % We will shift images at the end of preprocessing, no shifting before flatfielding...
    % with arbitary geometry, no need to force rotation axis to centre
    parameters.prep.udrift           = ydrift;
    parameters.prep.udriftabs        = ydrift;
    parameters.prep.vdrift           = zdrift;
    % to correct final reference group, add one more value
    parameters.prep.udrift(end+1)    = parameters.prep.udrift(end);
    parameters.prep.udriftabs(end+1) = parameters.prep.udriftabs(end);
    parameters.prep.vdrift(end+1)    = parameters.prep.vdrift(end);
end

save(parameters_name,'parameters');
save('backup.mat','parameters');
disp(' ')
disp('Parameters have been saved in file:')
disp(parameters_name)
disp(' ')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% launch preprocessing on OAR
% shifts will be applied in gtCreateFullLive
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('Ready to submit preprocessing jobs to OAR')
disp('can submit all, or select invidual jobs')

% Determine OAR parameters.
% Running "live" - wall time should be ~ acquisition time.  Running on a
% dataset that has already been collected - wall time should be estimate
% data treatment time.
% take into account mono tuning
if parameters.acq.mono_tune > 0
    extra_refs = ceil(totproj / (parameters.acq.refon * parameters.acq.mono_tune)) * parameters.acq.nref;
else
    extra_refs = 0;
end
if parameters.acq.interlaced_turns == 0
    filename = sprintf('%s%d0*.edf',parameters.acq.name,totproj/100);
else
    filename = sprintf('%s*_%d0*.edf',parameters.acq.name,totproj/100);
end
if ~isempty( dir( fullfile(parameters.acq.dir,'0_rawdata','Orig', filename ) ) )
    extra_proj = length( dir( fullfile(parameters.acq.dir,'0_rawdata','Orig', filename) ) );
else
    extra_proj = 0;
end
totdarks  = 1;
if parameters.acq.mono_tune > 0
    totqualis = (totproj/parameters.acq.refon)+1 + totproj/(parameters.acq.refon * parameters.acq.mono_tune);
else
    totqualis = (totproj/parameters.acq.refon)+1;
end
totrefs   = (((totproj/parameters.acq.refon)+1)*parameters.acq.nref);
totfiles  = totproj + totrefs + extra_refs + totqualis + totdarks + extra_proj; % extra images at the end if some option is activated in fasttomosetup

nfiles     = length(dir( fullfile(parameters.acq.dir,'0_rawdata','Orig',[parameters.acq.name '*.edf']) ));
nextrarefs = length(dir( fullfile(parameters.acq.dir,'0_rawdata','Orig','ref*_*99*.edf') ));
nrefs      = length(dir( fullfile(parameters.acq.dir,'0_rawdata','Orig','ref*_*00*.edf') ));
nqualis    = length(dir( fullfile(parameters.acq.dir,'0_rawdata','Orig','quali_*.edf') ));
ndarks     = length(dir( fullfile(parameters.acq.dir,'0_rawdata','Orig','dark*.edf') ));

if (nextrarefs + nrefs + nqualis + ndarks + nfiles) ~= totfiles
    disp('Something is wrong with the number of images in 0_rawdata/Orig... Please check')
end

tottime=totfiles*(parameters.acq.count_time+0.4)*1.3; %1.3 safety factor
% worst case E2V readout time

if nfiles>(totfiles*0.8)
    % finished - or nearly finished
    OARtime=(totproj/30)*15; % guess 15 seconds per job to be conservative
else
    % running - need to run until end of scan
    OARtime=tottime;
end

OARtmp=ceil(OARtime/60); % minutes
OARm=mod(OARtmp, 60); % minutes
OARh=(OARtmp-OARm)/60; % hours
walltime=sprintf('%02d:%02d:00', OARh, OARm);

OAR_parameters.walltime='6:00:00';
disp(['hard coded walltime to 6 hours instead of the estimated ' walltime ' hours...'])

disp('Numbers of OAR jobs launched are hardcoded')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
checkall=inputwdefault('Launch all jobs? [y/n]', 'y');
if strcmpi(checkall,'y')
    delimg = inputwdefault('Do you want to delete all the existing preprocessed images? [y/n]', 'y');
    if strcmpi(delimg,'y')
        delete(fullfile('0_rawdata',parameters.acq.name,'refHST*.edf'))
        delete(fullfile('1_preprocessing','abs','*.edf'))
        delete(fullfile('1_preprocessing','full','*.edf'))
        delete(fullfile('1_preprocessing','ext','*.edf'))
    end
    % gtSequenceMedianRefs : make refHST in raw data
    OAR_make('gtSequenceMedianRefs', 0, 0, 1, [parameters.acq.dir ' overwrite'], true,...
        'walltime',OAR_parameters.walltime);
    disp('gtSequenceMedianRefs launched')

    if (~parameters.acq.no_direct_beam)
        % gtCreateAbsLive
        OAR_make('gtCreateAbsLive', 0, totproj-1, 10, parameters.acq.dir, true,...
            'walltime',OAR_parameters.walltime);
        disp('gtCreateAbsLive launched')
        % gtAbsMedianLive
        OAR_make('gtAbsMedianLive', 0, totproj, 1, parameters.acq.dir, true,...
            'walltime',OAR_parameters.walltime);
        disp('gtAbsMedianLive launched')
    end

    % gtMovingMedianLive
    OAR_make('gtMovingMedianLive', 0, totproj, 1, parameters.acq.dir, true,...
        'walltime',OAR_parameters.walltime);
    disp('gtMovingMedianLive launched')
    % gtCreateFullLive
    OAR_make('gtCreateFullLive', 0, totproj-1, 30, parameters.acq.dir, true,...
        'walltime',OAR_parameters.walltime);
    disp('gtCreateFullLive launched')

else % or check before launching each job
    delimg = inputwdefault('Do you want to delete the existing corresponding images to the relative job? [y/n]', 'y');
    if strcmpi(delimg,'y')
        do_delete = true;
    else
        do_delete = false;
    end
    check = inputwdefault('Launch gtSequenceMedianRefs? [y/n]', 'y');
    if strcmpi(check,'y')
        if (do_delete), delete(fullfile('0_rawdata',parameters.acq.name,'refHST*.edf')); end
        OAR_make('gtSequenceMedianRefs', 0, 0, 1, [parameters.acq.dir ' overwrite'], true,...
            'walltime',OAR_parameters.walltime)
        disp('gtSequenceMedianRefs launched')
    end
    if (~parameters.acq.no_direct_beam)
        check = inputwdefault('Launch gtCreateAbsLive? [y/n]', 'y');
        if strcmpi(check,'y')
            if (do_delete), delete(fullfile('1_preprocessing','abs','abs*.edf')); end
            OAR_make('gtCreateAbsLive', 0, totproj-1, 15, parameters.acq.dir, true,...
                'walltime',OAR_parameters.walltime);
            disp('gtCreateAbsLive launched')
        end
        check = inputwdefault('Launch gtAbsMedianLive? [y/n]', 'y');
        if strcmpi(check,'y')
            if (do_delete), delete(fullfile('1_preprocessing','abs','med*.edf')); end
            OAR_make('gtAbsMedianLive', 0, totproj, 1, parameters.acq.dir, true,...
                'walltime',OAR_parameters.walltime);
            disp('gtAbsMedianLive launched')
        end
    end
    check = inputwdefault('Launch gtMovingMedianLive? [y/n]', 'y');
    if strcmpi(check,'y')
        if (do_delete), delete(fullfile('1_preprocessing','full','med*.edf')); end
        OAR_make('gtMovingMedianLive', 0, totproj, 1, parameters.acq.dir, true,...
            'walltime',OAR_parameters.walltime);
        disp('gtMovingMedianLive launched')
    end
    check = inputwdefault('Launch gtCreateFullLive? [y/n]', 'y');
    if strcmpi(check,'y')
        if (do_delete), delete(fullfile('1_preprocessing','full','full*.edf'));
            delete(fullfile('1_preprocessing','ext','*.edf')); end
        OAR_make('gtCreateFullLive', 0, totproj-1, 30, parameters.acq.dir, true,...
            'walltime',OAR_parameters.walltime);
        disp('gtCreateFullLive launched')
    end

end % end launch jobs with OAR

end % end of function