function mean_vals = gtCreateAbsLive(first_img, last_img, workingdirectory)
% GTCREATEABSLIVE Flatfield correction, normalisation and 2D medianfiltering to the absorption part
% mean_vals = gtCreateAbsLive(first_img, last_img, workingdirectory)
% ------------------------------------------------------------------
%
% When running it manually for a scan with, for example, 3600 images,
% choose first_img=0 and last_img=3604 to process the extra few images
% at the end of the scan.
%
%
% Version 003 10-04-2014 by PReischig
% Added capacity to process extra images at end of scan (for checks).
%
% Version 002 10-11-2011 by LNervo
% Add comments and clean syntax and warnings
%
% Sub-functions:
%[sub]- sfPutbb
%[sub]- sfGetbb
%
% Version 001
% if the scan has been done as many, broken acquistions, it is easier to the
% flat-fielding for each partial sequence using sequence_flatfield. In this
% case, add a field to parameters: prep.prior_flatfield.
% add the possibility to apply a horizontal shift to the images:
% prep.ydrift
% live! will wait for files to appear during the acquistion as a single
% condor job
%
% add test to wait for images to be translated by gtApplyDrifts before
% continuing - in edf_wait_read
% from prep.correct_drift
%
%
disp('gtCreateAbsLive.m')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% special case for running interactively in current directory
if ~exist('workingdirectory','var')
workingdirectory = pwd;
end
cd(workingdirectory);
parameters = gtLoadParameters();
acq = parameters.acq;
prep = parameters.prep;
datadir = fullfile(acq.dir, '0_rawdata', acq.name);
if (isdeployed)
global GT_DB %#ok<TLEV,NUSED>
global GT_MATLAB_HOME %#ok<NUSED,TLEV>
load('workspaceGlobal.mat');
first_img = str2double(first_img);
last_img = str2double(last_img);
end
% do not wait for images to be shifted by gtApplyDrifts - shift them after
% flatfield correction has been performed
waitfortranslated = false;
% if strcmp(prep.correct_drift, 'required')
% waitfortranslated=parameters.prep.correct_drift_iteration;
% end
%check if flatfielding is required
if (isfield(prep,'prior_flatfield') && prep.prior_flatfield == true)
do_flatfield = false;
else
do_flatfield = true;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do create abs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read the summed dark image and divide it by the number of frames
% we read in the boundingbox prep.bbox, which is slightly bigger than acq.bb
% the extra margin will be used for intensity normalisation of the images
% the saved absorption images will be of size acq.bb
if (do_flatfield)
darkend_name = fullfile(datadir, 'darkend0000.edf');
if (~exist(darkend_name, 'file'))
dark_name = fullfile(datadir, 'dark.edf');
dark = edf_wait_read(dark_name, prep.bbox, waitfortranslated);
else
dark = edf_wait_read(darkend_name, prep.bbox, waitfortranslated) / acq.ndark;
end
end
refname2old = [];
refname1old = [];
int0 = [];
% required for fullbeam normalisation with parallel processing
if (strcmp(prep.normalisation, 'fullbeam') && (first_img ~= 0))
% get the intensity of the first_img abs image
if (do_flatfield)
refname1 = fullfile(datadir, sprintf('refHST%04d.edf', 0) );
ref = edf_wait_read(refname1, prep.bbox, waitfortranslated) - dark;
name = fullfile(datadir, sprintf('%s%04d.edf', acq.name, 0) );
im = edf_wait_read(name, prep.bbox, waitfortranslated) - dark;
im = im ./ ref;
else
name = fullfile(datadir, sprintf('%sflat_%04d.edf', acq.name,0) );
im = edf_wait_read(name, prep.bbox, waitfortranslated);
end
int0 = gtImgMeanValue(im);
end
% no. of images in dataset
if strcmpi(acq.type, '360degree')
totproj = 2 * acq.nproj;
elseif strcmpi(acq.type, '180degree')
totproj = acq.nproj;
else
gtError('ACQ:invalid_parameter', ...
'Unknown type of scan. Check parameters.acq.type!');
end
% now loop through all the images: we do the following processing steps:
% - flatfield correction
% - normalization of the beam intensity outside the sample to be equal to
% one: allows to remove horizontal intensity variations in the images
% - we keep a record of the mean intensity in the image in the vector mean_vals(ii+1)
gtDBConnect();
for ii = first_img:last_img
if (do_flatfield) % if doing the flatfield in this step
% read image
% !!! is this not a problem with interlaced scans??
name = fullfile(datadir, sprintf('%s%04d.edf', acq.name, ii));
[im, info] = edf_wait_read(name, prep.bbox, waitfortranslated);
im = im - dark;
% check if we need to update the reference images
index1 = floor(ii / acq.refon) * acq.refon;
refname1 = fullfile(datadir, sprintf('refHST%04d.edf', index1) );
if (ii > (totproj-1))
if ~strcmp(refname1, refname1old)
ref = edf_wait_read(refname1, prep.bbox, waitfortranslated) - dark;
refname1old = refname1;
end
else
% if mono_tune = 3 and refon = 300, we have
% 0000, 0300, 0600, 0899, 0900, 1200, 1500, 1799, 1800, .......
if acq.mono_tune ~= 0
check = mod(ceil(ii / acq.refon), acq.mono_tune) == 0;
if check
index2 = ceil(ii / acq.refon) * acq.refon - 1;
index2 = max(0, index2);
refname2 = fullfile(datadir, sprintf('refHST%04d.edf', index2));
elseif mod(ceil(ii / acq.refon) * acq.refon, acq.mono_tune) == 0
index2 = ceil(ii / acq.refon) * acq.refon;
refname2 = fullfile(datadir, sprintf('refHST%04d.edf', index2));
end
else
index2 = ceil(ii / acq.refon) * acq.refon;
refname2 = fullfile(datadir, sprintf('refHST%04d.edf', index2));
end
fprintf('%d %d %d\n',index1, index2, ii);
if (~strcmp(refname1, refname1old) || ~strcmp(refname2, refname2old))
ref1 = edf_wait_read(refname1, prep.bbox, waitfortranslated) - dark;
ref2 = edf_wait_read(refname2, prep.bbox, waitfortranslated) - dark;
refname1old = refname1;
refname2old = refname2;
end
w = mod(ii, acq.refon) / acq.refon;
ref = (1 - w) .* ref1 + (w .* ref2);
end
% flatfield correction
im = im ./ ref;
else
%read flat fielded image
name = fullfile(datadir, sprintf('%sflat_%04d.edf', acq.name, ii));
[im, info] = edf_wait_read(name, prep.bbox, waitfortranslated);
end %end of if do flat field
%Now do further, special normalisation
% keep a record of the mean intensity, nomralise image intensity
switch (prep.normalisation)
case 'margin'
% standard way to do the normalisation
% calculate the vertical mean intensity profile outside the sample
% and renormalize this area to one
margins = [1:prep.margin, (size(im, 2) - prep.margin + 1):size(im, 2)];
mean_vals(ii + 1) = gtImgMeanValue( im(:, margins) );
mim = mean(im(:, margins), 2);
mim = repmat(mim, 1, size(im, 2));
im = im ./ mim;
case 'fullbeam'
% ak - mod to use central part of the image (within the sample) to normalise
if ii == 0
% Take the first image as the reference
int0 = gtImgMeanValue(im);
else
inti = gtImgMeanValue(im);
im = (int0/inti)*im;
end
mean_vals(ii + 1) = gtImgMeanValue(im);
case 'none'
% do nothing - mean_vals(ii) is obsolete and will not be used any more
mean_vals(ii + 1) = 1;
end
% apply a 2D median filter to reduce noise
im = gtImgMedianFilter2(im, prep.filtsize);
% we will adjust to the final size (acq.bb) after having shifted the images...
dummy = zeros(acq.ydet, acq.xdet);
im = sfPutbb(dummy, prep.bbox, im);
im = sfGetbb(im, acq.bb);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write out the abs images
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
absname = fullfile(acq.dir, '1_preprocessing', 'abs', sprintf('abs%04d.edf', ii) );
info.datatype = 'float32';
edf_write(im, absname, info);
end
fname = fullfile(acq.dir, '1_preprocessing', 'mean.mat');
save(fname,'mean_vals');
end
%%%%%%%%%%%%%%%%%%%%%%%
% SUB-FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%
% two little helper functions to put and get subimages defined by a
% gt_boundingbox from a bigger image
function a = sfGetbb(img, gtbb)
a = img(gtbb(2):(gtbb(2) + gtbb(4) - 1), gtbb(1):(gtbb(1) + gtbb(3) - 1));
end
function b = sfPutbb(img, gtbb, roi)
b = img;
b(gtbb(2):(gtbb(2) + gtbb(4) - 1), gtbb(1):(gtbb(1) + gtbb(3) - 1)) = roi;
end