Skip to content
Snippets Groups Projects
Commit da0ac768 authored by Nicola Vigano's avatar Nicola Vigano
Browse files

Small formatting changes to pre-processing functions 


Signed-off-by: default avatarNicola Vigano <nicola.vigano@esrf.fr>
parent 295a1d1f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
1_preprocessing/gtCopy.m
+ 129
− 150
View file @ da0ac768
......@@ -26,174 +26,153 @@ function gtCopy(fnamein, fnameout, fnameout_orig, parameters, move_data_flag)
%
%
fprintf('gtCopy of "%s"\n', fnamein);
fprintf('gtCopy of "%s"\n', fnamein);
[~, ~, extension] = fileparts(fnamein);
switch (extension)
case '.edf'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% image corrections for edf files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
copyEdf(fnamein, fnameout, parameters);
otherwise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% non-images files: just copy
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
try
[~, msg] = copyfile(fnamein, fnameout); disp(msg);
catch mexc
gtPrintException(mexc);
mexc1 = GtException('SETUP:copy_error', ...
'Not able to copy "%s" to: "%s".', fnamein, fnameout);
mexc1 = mexc1.addCause(mexc);
throw(mexc1);
end
end
[~, ~, extension] = fileparts(fnamein);
switch (extension)
case '.edf'
% image corrections for edf files
copyEdf(fnamein, fnameout, parameters);
otherwise
% non-images files: just copy
try
[~, msg] = copyfile(fnamein, fnameout); disp(msg);
catch mexc
gtPrintException(mexc);
mexc1 = GtException('SETUP:copy_error', ...
'Not able to copy "%s" to: "%s".', fnamein, fnameout);
mexc1 = mexc1.addCause(mexc);
throw(mexc1);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% move or copy the original file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (move_data_flag)
try
[~, msg] = movefile(fnamein, fnameout_orig); disp(msg);
catch mexc
% if no permission to mv, copy instead
errorMsg = ['Not able to move "' fnamein '" to: "' fnameout_orig '".'];
gtPrintException(mexc, errorMsg);
disp('Trying to just copy...');
% move or copy the original file
if (move_data_flag)
try
[~, msg] = copyfile(fnamein, fnameout_orig); disp(msg);
catch mexc2
gtPrintException(mexc2);
mexc1 = GtException('SETUP:copy_error', ...
'Not able to move or copy "%s" to: "%s".', fnamein, fnameout);
mexc1 = mexc1.addCause(mexc);
mexc1 = mexc1.addCause(mexc2);
throw(mexc1);
[~, msg] = movefile(fnamein, fnameout_orig); disp(msg);
catch mexc
% if no permission to mv, copy instead
errorMsg = ['Not able to move "' fnamein '" to: "' fnameout_orig '".'];
gtPrintException(mexc, errorMsg);
disp('Trying to just copy...');
try
[~, msg] = copyfile(fnamein, fnameout_orig); disp(msg);
catch mexc2
gtPrintException(mexc2);
mexc1 = GtException('SETUP:copy_error', ...
'Not able to move or copy "%s" to: "%s".', fnamein, fnameout);
mexc1 = mexc1.addCause(mexc);
mexc1 = mexc1.addCause(mexc2);
throw(mexc1);
end
end
end
end
end %end of function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sub-functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function copyEdf(fnamein, fnameout, parameters)
[img, info] = edf_wait_read(fnamein);
[img, info] = edf_wait_read(fnamein);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ccd specific corrections
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (ismember(lower(parameters.acq.sensortype), {'kodak4mv1', 'frelon 4m', 'frelon hd4m'}))
%add fields to info
if ( (isfield(info, 'isccdcorrected') && strcmpi(info.isccdcorrected, 'true')) ...
|| (isfield(info, 'israw') && strcmpi(info.israw, 'false')) )
disp('Kodak corrections already applied');
else
disp('applying Kodak chip corrections');
info.ccdmodel = 'Kodak4Mv1';
info.isccdcorrected = true;
% do the kodak chip corrections
if parameters.acq.xdet==2048 && parameters.acq.ydet==2048
midline = round(parameters.acq.ydet/2);
img = [img(1:midline,:); ...
mean([img(midline,:); img(midline+1,:)]); ...
img(midline+1:end-1,:) ];
% im=im.^0.96;
elseif parameters.acq.xdet==1024 && parameters.acq.ydet==1024
disp('assuming 2x2 binning')
% scale up, add the line, scale down
img = imresize(img, 2);
midline = 1024;
img = [img(1:midline,:); ...
mean([img(midline,:); img(midline+1,:)]); ...
img(midline+1:end-1,:) ];
img = imresize(img, 0.5);
% ccd specific corrections
if (ismember(lower(parameters.acq.sensortype), {'kodak4mv1', 'frelon 4m', 'frelon hd4m'}))
%add fields to info
if ( (isfield(info, 'isccdcorrected') && strcmpi(info.isccdcorrected, 'true')) ...
|| (isfield(info, 'israw') && strcmpi(info.israw, 'false')) )
disp('Kodak corrections already applied');
else
disp('applying Kodak chip corrections');
info.ccdmodel = 'Kodak4Mv1';
info.isccdcorrected = true;
% do the kodak chip corrections
if (parameters.acq.xdet == 2048 && parameters.acq.ydet == 2048)
midline = round(parameters.acq.ydet/2);
img = [img(1:midline,:); ...
mean([img(midline,:); img(midline+1,:)]); ...
img(midline+1:end-1,:) ];
% im=im.^0.96;
elseif (parameters.acq.xdet == 1024 && parameters.acq.ydet == 1024)
disp('assuming 2x2 binning')
% scale up, add the line, scale down
img = imresize(img, 2);
midline = 1024;
img = [img(1:midline,:); ...
mean([img(midline,:); img(midline+1,:)]); ...
img(midline+1:end-1,:) ];
img = imresize(img, 0.5);
end
end
else
%add fields to info
info.ccdmodel = parameters.acq.sensortype;
end
else
%add fields to info
info.ccdmodel = parameters.acq.sensortype;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% distortion correction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (strcmp(parameters.acq.distortion, 'none'))
disp('No correction applied.');
elseif (isfield(info, 'isundistorted') && strcmpi(info.isundistorted, true))
disp('Image has already been corrected.');
else
try
% get distortion map - passed in or read
fprintf('Reading distortion map from file: %s\n', ...
parameters.acq.distortion);
detector_size = [];
detector_size.x = parameters.acq.true_detsizeu;
detector_size.y = parameters.acq.true_detsizev;
if ( detector_size.x == parameters.acq.xdet ...
&& detector_size.y == parameters.acq.ydet )
fprintf('- Using full detector..');
distortion = distortion_read(parameters.acq.distortion, detector_size);
else
fprintf('- Using a ROI of the detector (%d, %d)..', ...
parameters.acq.xdet, parameters.acq.ydet)
% ROI is of the form [xmin ymin x_size y_size]
roi = [ parameters.acq.detroi_u_off parameters.acq.detroi_v_off ...
parameters.acq.xdet parameters.acq.ydet ];
distortion = distortion_read(parameters.acq.distortion, detector_size, roi);
% distortion correction
if (strcmp(parameters.acq.distortion, 'none'))
disp('No correction applied.');
elseif (isfield(info, 'isundistorted') && strcmpi(info.isundistorted, true))
disp('Image has already been corrected.');
else
try
% get distortion map - passed in or read
fprintf('Reading distortion map from file: %s\n', ...
parameters.acq.distortion);
detector_size = [];
detector_size.x = parameters.acq.true_detsizeu;
detector_size.y = parameters.acq.true_detsizev;
if ( detector_size.x == parameters.acq.xdet ...
&& detector_size.y == parameters.acq.ydet )
fprintf('- Using full detector..');
distortion = distortion_read(parameters.acq.distortion, detector_size);
else
fprintf('- Using a ROI of the detector (%d, %d)..', ...
parameters.acq.xdet, parameters.acq.ydet)
% ROI is of the form [xmin ymin x_size y_size]
roi = [ parameters.acq.detroi_u_off parameters.acq.detroi_v_off ...
parameters.acq.xdet parameters.acq.ydet ];
distortion = distortion_read(parameters.acq.distortion, detector_size, roi);
end
fprintf(' Done.\nApplying distortion correction..');
img = distortion_correct(img, distortion);
info.isundistorted = true;
fprintf(' Done.\n');
catch mexc
errorMsg = [ 'Not able to apply correction from file: "' ...
parameters.acq.distortion '".' ];
gtPrintException(mexc, errorMsg);
end
fprintf(' Done.\nApplying distortion correction..');
img = distortion_correct(img, distortion);
info.isundistorted = true;
fprintf(' Done.\n');
catch mexc
errorMsg = [ 'Not able to apply correction from file: "' ...
parameters.acq.distortion '".' ];
gtPrintException(mexc, errorMsg);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% turn images for gt coordinate system
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (isfield(parameters.acq, 'rotation_axis') ...
&& strcmp(parameters.acq.rotation_axis, 'horizontal') )
if (~isfield(info, 'gtrotated') || strcmp(info.gtrotated, 'false'))
if strcmpi(parameters.acq.rotation_direction, 'clockwise')
disp('Rotating 90 degrees counterclockwise (horizontal rotation axis).');
info.gtrotated = true;
img = rot90(img, -1); % 90 degrees
elseif strcmpi(parameters.acq.rotation_direction, 'counterclockwise')
disp('Rotating 90 degrees clockwise (horizontal rotation axis).');
info.gtrotated = true;
img = rot90(img); % -90 degrees
else
disp('Options are: "counterclockwise" and "clockwise".');
% turn images for gt coordinate system
if (isfield(parameters.acq, 'rotation_axis') ...
&& strcmp(parameters.acq.rotation_axis, 'horizontal') )
if (~isfield(info, 'gtrotated') || strcmp(info.gtrotated, 'false'))
if strcmpi(parameters.acq.rotation_direction, 'clockwise')
disp('Rotating 90 degrees counterclockwise (horizontal rotation axis).');
info.gtrotated = true;
img = rot90(img, -1); % 90 degrees
elseif strcmpi(parameters.acq.rotation_direction, 'counterclockwise')
disp('Rotating 90 degrees clockwise (horizontal rotation axis).');
info.gtrotated = true;
img = rot90(img); % -90 degrees
else
disp('Options are: "counterclockwise" and "clockwise".');
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% flip images for gt coordinate system
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if parameters.acq.flip_images==true
if (~isfield(info, 'gtlrflipped') || strcmp(info.gtlrflipped, 'false') || info.gtlrflipped==false)
disp('Flipping LR for gt coordinate system.');
info.gtlrflipped = true;
img = fliplr(img);
% flip images for gt coordinate system
if (parameters.acq.flip_images)
if (~isfield(info, 'gtlrflipped') || strcmpi(info.gtlrflipped, 'false') || ~info.gtlrflipped)
disp('Flipping LR for gt coordinate system.');
info.gtlrflipped = true;
img = fliplr(img);
end
else
disp('Any flip has been applied to the images.')
end
else
disp('Any flip has been applied to the images.')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write the corrected image and info to the new location
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
edf_write(img, fnameout, info);
end % end copyEdf
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write the corrected image and info to the new location
edf_write(img, fnameout, info);
end
1_preprocessing/gtCopyCorrectUndistortCondor.m
+ 3
− 4
View file @ da0ac768
......@@ -19,8 +19,8 @@ function gtCopyCorrectUndistortCondor(~, ~, workingdirectory)
%
if isdeployed
global GT_MATLAB_HOME
global GT_DB
global GT_MATLAB_HOME %#ok<TLEV,NUSED>
global GT_DB %#ok<TLEV,NUSED>
load('workspaceGlobal.mat')
end
......@@ -28,8 +28,7 @@ end
% load and read data from parameters.mat
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cd(workingdirectory)
parameters = [];
load('parameters.mat');
parameters = gtLoadParameters();
source_dir = parameters.acq.collection_dir;
dest_dir = fullfile(parameters.acq.dir, '0_rawdata', parameters.acq.name);
......
1_preprocessing/gtCreateAbsLive.m
+ 47
− 51
View file @ da0ac768
......@@ -4,7 +4,7 @@ function 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
% choose first_img=0 and last_img=3604 to process the extra few images
% at the end of the scan.
%
%
......@@ -33,8 +33,6 @@ function mean_vals = gtCreateAbsLive(first_img, last_img, workingdirectory)
%
%
parameters=[];
disp('gtCreateAbsLive.m')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% get parameters
......@@ -45,14 +43,15 @@ if ~exist('workingdirectory','var')
end
cd(workingdirectory);
load('parameters.mat');
parameters = gtLoadParameters();
acq = parameters.acq;
prep = parameters.prep;
datadir = fullfile(acq.dir, '0_rawdata', acq.name);
if (isdeployed)
global GT_DB
global GT_MATLAB_HOME
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);
......@@ -83,7 +82,7 @@ end
% the saved absorption images will be of size acq.bb
if (do_flatfield)
darkend_name = fullfile(datadir, 'darkend0000.edf');
if ~exist(darkend_name, 'file')
if (~exist(darkend_name, 'file'))
dark_name = fullfile(datadir, 'dark.edf');
dark = edf_wait_read(dark_name, prep.bbox, waitfortranslated);
else
......@@ -97,9 +96,9 @@ 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
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) );
......@@ -112,7 +111,6 @@ if strcmp(prep.normalisation, 'fullbeam') && first_img ~= 0
int0 = gtImgMeanValue(im);
end
% no. of images in dataset
if strcmpi(acq.type, '360degree')
totproj = 2 * acq.nproj;
......@@ -120,10 +118,9 @@ elseif strcmpi(acq.type, '180degree')
totproj = acq.nproj;
else
gtError('ACQ:invalid_parameter', ...
'Unknown type of scan. Check parameters.acq.type!');
'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
......@@ -135,7 +132,7 @@ 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??
% !!! 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);
......@@ -145,41 +142,41 @@ for ii = first_img:last_img
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
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
......@@ -229,14 +226,13 @@ 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 % end of for-loop
end
fname = fullfile(acq.dir, '1_preprocessing', 'mean.mat');
save(fname,'mean_vals');
end % end of function
end
%%
%%%%%%%%%%%%%%%%%%%%%%%
% SUB-FUNCTIONS
%%%%%%%%%%%%%%%%%%%%%%%
......
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