gtSequenceMedianRefs.m

function gtSequenceMedianRefs(~, ~, workingdirectory, overwrite_flag)
% GTSEQUENCEMEDIANREFS  Generate median image for each set of reference images in a scan
%     gtSequenceMedianRefs(first, last, workingdirectory, overwrite_flag)
%     -------------------------------------------------------------------
%       Works on reference files found in current directory, and only
%       overwrites if explicitly told to.
%
%     Version 003 22-02-2007 by GJ
%       takes arguments 'offset',-1 or 1 : will look for offsets below (or above)
%       those in the xml file
%       'overwrite' : will overwrite pre-existing files
%
%     Version 002 15-06-2007 by AKing
%       final modifications - offests only -1, but not at every reference group.
%       Always wait, but skip if the next reference group appears
%       overwrite_flag should be the string 'overwrite' to overwrite references
%
%       THIS SHOULD BE COMPILED WITH ID19_MCC TO PUT IN THE RIGHT PLACE
%       app.overwrite=false;  %single arg...
%       app.wait=true;  %always wait (see above)
%       app.offset=0;   deal with -1 and 0  (so 0, 50, 99, 100, 150 etc)
%       app=parse_pv_pairs(app,varargin);
%       offsets=sort(unique([app.offset 0 app.offset]));
%
%       gt version - if parameters.prep.correct_drift=='required', function must
%       wait for 'istranslated'==true in the edf header
%
%       modify to work from the parameters file (which is created from the .xml)
%
%       propagate edf header info to the refHST created
%
%

disp('gtSequenceMedianRefs.m  -  graintracking specific version')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load('parameters.mat');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~exist('workingdirectory', 'var') || isempty(workingdirectory)
    workingdirectory=pwd;
end
cd(workingdirectory)

if exist('parameters.mat', 'file')
    load('parameters.mat');
    acq=parameters.acq;
else
    disp('gtSequenceMedianRefs should be launched in the analysis directory (where the parameters.mat file is)')
    disp('quitting...')
    return
end

% overwrite preexisting refs
if ~exist('overwrite_flag','var')
  overwrite_flag='';
end

% should we wait for images to be shifted by gtApplyDrifts?
waitfortranslated=false;
% if strcmp(parameters.prep.correct_drift, 'required')
%     waitfortranslated=parameters.prep.correct_drift_iteration;
% end

%images per turn
if strcmp(acq.type, '360degree')
    nimages=2*acq.nproj;
elseif strcmp(acq.type, '180degree')
    nimages=acq.nproj;
else
    Mexc = MException('ACQ:invalid_parameter', ...
                      'Unknown type of scan. Check parameters.acq.type!');
    throw(Mexc);
end

%offsets - for situations where the mono is tuned at certain reference
%groups
offsets=[-1 0];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% move into the directory with the images
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd([acq.dir '/0_rawdata/' acq.name]);


for n=0:acq.refon:nimages
  for offset=offsets

    if (n+offset<0)
      continue
    end

    fname_refC=sprintf('refHST%04d.edf',n+offset);
    if exist(fname_refC,'file') && ~strcmpi(overwrite_flag, 'overwrite')
      fprintf('%s already exists - skipping\n',fname_refC);
      continue
    end

    fprintf('Looking for references at %04d\n',n+offset);
    ref=zeros(acq.nref,acq.ydet,acq.xdet);

    skip=false;

    %determine next reference group after this one - if it appears, move on
    if offset==-1
      fname_next=sprintf('ref0000_%04d.edf', n);
    else
      fname_next=[];
    end

    for m=0:acq.nref-1
      fname_ref=sprintf('ref%04d_%04d.edf',m,n+offset);

      if ~exist(fname_ref,'file')
        if isempty(fname_next)
          fprintf('Waiting for %s...',fname_ref)
        else
          fprintf('Waiting for %s or %s...',fname_ref, fname_next)
        end
        %modif andy and sabine to avoid unfinished jobs in condor queue
        i=0;
        while ~exist(fname_ref,'file') && ~exist(fname_next,'file')
            i=i+1;
          pause(2)
          if i>2000 % timeout
              fprintf('\nTimeout. Exiting...\n')
			  cd(parameters.acq.dir);
              return
          end
        end
        fprintf('Found file.\n')
      end

      if exist(fname_ref,'file')
        fprintf('\tReading image %d of %d\r',m+1,acq.nref)
        ref(m+1,:,:) = edf_wait_read(fname_ref, [], waitfortranslated);
        info=edf_info(fname_ref);
      elseif exist(fname_next,'file')
        skip=true;
        break
      end

    end

    if skip==false
      fprintf('Now performing median...\n')
      if (parameters.acq.nref == 1)
          refC = squeeze(ref(1,:,:));
      else
          refC=squeeze(median(ref));
      end
      info.datatype='uint16';
	  fprintf('Saving file: %s\n', fname_refC)
      edf_write(refC,fname_refC, info);
    end

  end
end

cd(parameters.acq.dir);

end % end function