Skip to content
Snippets Groups Projects
Commit cfd6062d authored by Stefan Schmiederer's avatar Stefan Schmiederer
Browse files

Fixed mean_vals bug; trivial formatting. 


Signed-off-by: default avatarStefan Schmiederer <sschmied@esrf.fr>
parent 758bccd3
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
1_preprocessing/gtCreateAbsLive.m
+ 44
− 66
View file @ cfd6062d
......@@ -33,7 +33,7 @@ disp('gtCreateAbsLive.m')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% special case for running interactively in current directory
if ~exist('workingdirectory','var')
workingdirectory=pwd;
workingdirectory = pwd;
end
cd(workingdirectory);
......@@ -58,19 +58,6 @@ waitfortranslated = false;
% waitfortranslated=parameters.prep.correct_drift_iteration;
% end
% how many images in this scan (for renumbering of interlaced scans)
if strcmp(acq.type, '360degree')
totproj = 2*acq.nproj;
elseif strcmp(acq.type, '180degree')
totproj = acq.nproj;
else
Mexc = MException('ACQ:invalid_parameter', ...
'Unknown type of scan. Check parameters.acq.type!');
throw(Mexc);
end
interlacedscan = parameters.acq.interlaced_turns + 1;
imageperturn = totproj / (interlacedscan);
%check if flatfielding is required
if (isfield(prep,'prior_flatfield') && prep.prior_flatfield == true)
do_flatfield = false;
......@@ -90,42 +77,26 @@ if do_flatfield
darkend_name = fullfile(datadir, 'darkend0000.edf');
if ~exist(darkend_name, 'file')
dark_name = fullfile(datadir, 'dark.edf');
d = edf_wait_read(dark_name, prep.bbox, waitfortranslated);
dark = edf_wait_read(dark_name, prep.bbox, waitfortranslated);
else
d = edf_wait_read(darkend_name, prep.bbox, waitfortranslated)/acq.ndark;
dark = edf_wait_read(darkend_name, prep.bbox, waitfortranslated) / acq.ndark;
end
end
%%%%%%%%%%%%%%%
% 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
%%%%%%%%%%%%%%%%
refname2old = [];
refname1old = [];
int0 = [];
% required for fullbeam normalisation with parallel processing
if strcmp(prep.normalisation, 'fullbeam') && first_img~=0
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)-d;
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)-d;
im = im./ref;
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);
......@@ -147,60 +118,53 @@ for ii=first_img:last_img
name = fullfile(datadir, sprintf('%s%04d.edf', acq.name, ii));
[im, info] = edf_wait_read(name, prep.bbox, waitfortranslated);
im = im -d;
% %ak - this is repeated later....
% % keep a record of the mean intensity
% mean_vals(ii+1)=gtImgMeanValue(im(:,[1:prep.margin,end-prep.margin+1:end]));
im = im - dark;
% check if we need ot update the reference images
index1 = floor(ii/acq.refon)*acq.refon;
index1 = floor(ii / acq.refon) * acq.refon;
refname1 = fullfile(datadir, sprintf('refHST%04d.edf', index1) );
% take mono_tune into account wl 01/2008
if ( isfield(acq,'mono_tune') && acq.mono_tune ~= 0 ... % sabine if mono_tune
&& mod(ceil(ii/acq.refon)*acq.refon,acq.mono_tune) == 0 )
index2 = ceil(ii/acq.refon)*acq.refon-1;
if ( isfield(acq, 'mono_tune') && acq.mono_tune ~= 0 ... % sabine if mono_tune
&& mod(ceil(ii / acq.refon) * acq.refon, acq.mono_tune) == 0 )
index2 = ceil(ii / acq.refon) * acq.refon - 1;
index2 = max(0, index2);
refname2 = fullfile(datadir, sprintf('refHST%04d.edf', index2));
else
index2 = ceil(ii/acq.refon)*acq.refon;
index2 = ceil(ii / acq.refon) * acq.refon;
refname2 = fullfile(datadir, sprintf('refHST%04d.edf', index2));
end
if (~strcmp(refname1, refname1old) || ~strcmp(refname2, refname2old))
ref1 = edf_wait_read(refname1, prep.bbox, waitfortranslated)-d;
ref2 = edf_wait_read(refname2, prep.bbox, waitfortranslated)-d;
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;
w = mod(ii, acq.refon) / acq.refon;
ref = (1 - w) .* ref1 + (w .* ref2);
% flatfield correction
im = im./ref;
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);
[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'
% should be the standard way to do it - uses prep.margin pixels left and right from the sample BoundingBox
% n1=acq.bb(1)-prep.bbox(1)+1;
% n2=prep.bbox(1)+prep.bbox(3)-acq.bb(1)-acq.bb(3);
% mean_vals(ii+1)=gtImgMeanValue(im(:,[n1:n1+prep.margin-1,end-n2-prep.margin+1:end-n2]));
% 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) );
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;
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
% 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);
......@@ -208,14 +172,14 @@ for ii=first_img:last_img
inti = gtImgMeanValue(im);
im = (int0/inti)*im;
end
mean_vals(ii+1) = gtImgMeanValue(im);
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;
mean_vals(ii + 1) = 1;
end
% apply a 2D median filter to reduce noise
im = gtImgMedianFilter2(im,prep.filtsize);
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);
......@@ -228,10 +192,24 @@ for ii=first_img:last_img
absname = fullfile(acq.dir, '1_preprocessing', 'abs', sprintf('abs%04d.edf', ii) );
info.datatype = 'float32';
edf_write(im, absname, info);
end % for
end % end of for-loop
fname = fullfile(acq.dir, '1_preprocessing', 'mean.mat');
save(fname,'mean_vals');
end % end of function
%%%%%%%%%%%%%%%
% 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
%%%%%%%%%%%%%%%%
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment