Skip to content

GitLab

Commit 1d52afdc authored by Leon Merten Lohse's avatar Leon Merten Lohse
Browse files

Update template_tomo_analysis_GINIX.m

parent 8e3c2a53
Pipeline #109194 passed with stage
in 1 minute and 39 seconds
Hide whitespace changes
Inline Side-by-side
examples/template_tomo_analysis_GINIX.m
View file @ 1d52afdc
%% Start from scratch
clear 'variables'; clc;
set(0,'DefaultFigureColor','w');
%% Parameters of the scan
year = '2016';
runname = 'run56_LTP';
date = '20161116';
detector = 'pirra';
filePrefix = '12816-6_Kleinhirn3_tomo7';
%% (modify if loading fails)
if ispc
toolboxPath = 'S:/Projects_X-ray_Imaging/holotomotoolbox/functions';
astraPath = 'S:/Projects_X-ray_Imaging/ASTRAToolbox/windows';
prePath = 'S:/Messzeiten/';
elseif isunix
toolboxPath = '/IRP/AG_Salditt/Projects_X-ray_Imaging/holotomotoolbox/functions';
astraPath = '/IRP/AG_Salditt/Projects_X-ray_Imaging/ASTRAToolbox/linux';
prePath = '/IRP/AG_Salditt/Messzeiten/';
end
dataPath = fullfile(prePath, year, 'GINIX', runname, '/data/current/raw', date, 'detectors', detector, filePrefix);
if ~exist(dataPath,'dir')
error(['Assumed data path ''', dataPath, ''' does not exist. Scan parameters might be wrong. Please fix the path manually.']);
end
if ~exist(toolboxPath,'dir')
error(['Assigned path to the holotomotoolbox ''', toolboxPath, ''' does not exist.']);
end
if ~exist(astraPath,'dir')
error(['Assigned path to the ASTRA-toolbox ''', astraPath, ''' does not exist.']);
end
%% Add holotomotoolbox to path
addpath(genpath(toolboxPath));
%% Add ASTRA-toolbox to path
addpath(genpath(astraPath));
%% Set up image-reader for the specified detector
imageReader = getImageReaderIRP(detector);
switch(lower(detector))
case 'argos'
fileEnding = '_%04i.tif'; dumpEnding = '_%04i.tif.txt';
dx = 0.54e-6;
case 'talos'
fileEnding = '_%04i.tif'; dumpEnding = '_%04i.tif.txt';
dx = 75e-6;
case 'pirra'
fileEnding = '_%i.tif'; dumpEnding = '_%04i.tif.txt';
dx = 6.5e-6;
case 'thyia'
fileEnding = '_%i.tif'; dumpEnding = '_%04i.tif.txt';
dx = 4.5e-6;
case 'iris'
fileEnding = '_%04i.img'; dumpEnding = '_%04i.img.txt';
dx = 6.5e-6;
case 'timepix'
fileEnding = '_%04i.raw'; dumpEnding = '_%04i.raw.txt';
dx = 55e-6;
case 'zyla'
fileEnding = '_%04i.tif'; dumpEnding = '_%04i.tif.txt';
dx = 6.5e-6;
end
fileName = [filePrefix,fileEnding];
dumpName = [filePrefix,dumpEnding];
getFileName = @(number) fullfile(dataPath, sprintf(fileName,number));
getDumpName = @(number) fullfile(dataPath,'dump',sprintf(dumpName,number));
%% Read one image to check if path is correct
number = 100;
raw0 = imageReader(getFileName(number));
showImage(raw0);
colormap gray
%% setup the current run - edit here
fprintf('... add specifications of setup... \n' ) ;
angularRange = 180;
E = 8000; % energy in eV
lambda = 12.398/E*1e-10; % wavelength in m
% detector distance
z02 = 5047; % in mm
% calibration motor vs real distance to wg
xmotor = 'stx';
delta_x = 0; % width of the waveguide (silicon channels)
corr_x = 0; % width of the vertical 1dWG (crossed multilayer)
corr_y = 0; % width of the horizontal 1dWG (crossed multilayer)
%% set up file numbers - edit here
fprintf('... set up filenumbers... \n' ) ;
numDistances = 4; % number of defocus distances
numAngles = 1501; % number of projections
numFlats = 50; % number of flat fields
numDarks = 50; % number of dark fields
offset=84; % number of images prior to flatfields (wg_align)
% Compute file numbers of projections, flatfields and darkfields - do not edit
numImages=offset+2*numFlats+numAngles; %number of images per distance scan
flatFilenum = zeros([2,numDistances,numFlats]);
projFilenum = zeros([numDistances,numAngles]);
for i=1:numDistances
flatFilenum(1,i,:) = offset+1+(i-1)*numImages:offset+numFlats+(i-1)*numImages;
projFilenum(i,:) = flatFilenum(1,i,end)+1:flatFilenum(1,i,end)+numAngles;
flatFilenum(2,i,:) = projFilenum(i,end)+1:projFilenum(i,end)+numFlats;
end
darknum = flatFilenum(2,end,end)+1:flatFilenum(2,end,end)+numDarks;
%% Check ranges of projection-, flat-field- and dark-field images
fprintf('... please check flat, projection and dark... \n' ) ;
for iii = 1 : numDistances
tmp = imageReader(getFileName(flatFilenum(1,iii,1)));
showImage(tmp);
title('Is this a flatfield? (Y- enter) (N - Ctr+c)')
colormap gray
pause()
tmp = imageReader(getFileName(projFilenum(iii,1)));
showImage(tmp);
title('Is this a projection? (Y- enter) (N - Ctr+c)')
colormap gray
pause()
end
tmp = imageReader(getFileName(darknum(1)));
showImage(tmp);
title('Is this a darkfield? (Y- enter) (N - Ctr+c)')
colormap gray
pause()
%% read flat fields - do not edit
fprintf('... read flatfields... \n' ) ;
empty = zeros([size(raw0,1),size(raw0,2),numDistances,2], 'single');
for indDistance = 1:numDistances
for indEmpty = 1:size(flatFilenum,1)
emptyTmp = zeros([size(raw0,1),size(raw0,2),numFlats], 'single');
for indImg = 1:numFlats
emptyTmp(:,:,indImg) = imageReader(getFileName(flatFilenum(indEmpty,indDistance,indImg)));
fprintf('.');
end
% take median of all emptys
empty(:,:,indDistance,indEmpty) = median(emptyTmp,3);
fprintf('\n');
end
end
%% read dark - do not edit
fprintf('... read darkfields... \n' ) ;
darkTmp = zeros(size(raw0,1),size(raw0,2),numDarks);
for indImg = 1:numDarks
darkTmp(:,:,indImg) = imageReader(getFileName(darknum(indImg)));
fprintf('.');
end
% take median of all darkfield images
dark = median(darkTmp,3);
fprintf('\n');
%% calculate geometry parameters - do not edit
fprintf('... calculate geometry parameters... \n' ) ;
% Read defocus distances from dump files
stx = zeros([1,numDistances]);
for distIdx = 1:numDistances
stx(distIdx) = wwas(getDumpName(projFilenum(distIdx,2)),xmotor);
end
% calculate derived parameters: (suffixes 'x' and 'y' account for
% possible astigmatism)
%
% true defocus distances including corrections
z01x = stx+corr_x+delta_x;
z01y = stx+corr_y+delta_x;
%
% distances between object and detetor
z12x = z02-z01x;
z12y = z02-z01y;
%
% magnifications
Mx = z02./z01x;
My = z02./z01y;
%
% effective distance between object and detetor
z_effx = z12x./Mx;
z_effy = z12y./My;
%
% effective pixel size
dxeffx = dx./Mx;
dxeffy = dx./My;
%
% Fresnel numbers in the effective geometry
Fx = dxeffx.^2./(z_effx*lambda);
Fy = dxeffy.^2./(z_effy*lambda);
% Magnifications and Fresnel numbers for all defocus distances (columns),
% possibly different along the x- and y-dimensions of the images (rows)
M = [My(:).'; Mx(:).'];
F = [Fy(:).'; Fx(:).'];
%% Read, rescale and align holograms
% Which distances should be read?
distToUse = 1:numDistances;
% Settings for alignment of the different defocused images. See
% documentation of rescaleDefocusSeries for details and additional options
alignSettings = rescaleDefocusSeries;
alignSettings.alignMethod = 'dftregistration';
alignSettings.medfilt = true; % median filter
alignSettings.sigmaHighpass = 50; % high pass filter
alignSettings.alignReconstructedImages = true; % should the holograms or single distance CTF reconstructed images be sued for registration?
% settings for the preliminary single-distance CTF-reconstruction used for aligning the images
ctfsettingsAlign = phaserec_ctf;
ctfsettingsAlign.betaDeltaRatio = 1/10;
ctfsettingsAlign.lim2 = 0.1;
alignSettings.ctfsettings = ctfsettingsAlign;
% no changes below this point
fprintf('... reading holograms... \n' ) ;
holograms = zeros([size(raw0,1), size(raw0,2), numel(distToUse), numAngles], 'single');
thetas = zeros([numAngles,1]);
parfor indProj = 1:numAngles
fprintf('projection: %i \n',indProj)
% read in all distances
I = zeros([size(raw0,1), size(raw0,2), length(distToUse)]);
idx = 0;
for indDistance = distToUse
idx = idx+1;
% interpolate linearly between the two empty images
a = (numAngles-indProj)/(numAngles-1);
emptyInterpolated = a*empty(:,:,indDistance,1)+(1-a)*empty(:,:,indDistance,2);
% read in data and perform empty beam correction
tmp = imageReader(getFileName(projFilenum(indDistance,indProj)));
tmpCorr=(tmp-dark)./(emptyInterpolated-dark);
% remove hot pixels
I(:,:,idx) = removeOutliers(tmpCorr);
end
% read rotation angle
if indProj~=1
thetas(indProj) = wwas(getDumpName(projFilenum(distToUse(1),indProj)),'stzrot');
else
thetas(indProj)=0;
end
% rescale and align the images recorded at different defocus distances
IAligned = rescaleDefocusSeries(I, F(:,distToUse), M(:,distToUse), alignSettings);
% Store aligned hologram
holograms(:,:,:,indProj) = IAligned;
end
% Compute rescaled Fresnel-numbers (cannot be done in parfor-loop)
[~, fresnelNumbers] = rescaleDefocusSeries(ones([2,2,numel(distToUse)]), F(:,distToUse), M(:,distToUse), alignSettings);
%% Show all holograms for a random tomographic angle to check results
projIdx = randi(numAngles);
for jj = 1:size(holograms,3)
showImage(holograms(:,:,jj,projIdx));
colormap gray
pause(1)
end
%% Show full tomographic data set for first distance to check results
figure; showStack(squeeze(holograms(:,:,1,:)));
%% OPTIONAL: delete corrupted holograms due to missing wedge from data
missingWedge = []; % Enter range of corrupted holograms here
holograms(:,:,:,missingWedge) = [];
thetas(missingWedge) = [];
%% OPTIONAL: high-pass filtering of holograms to reduce possible artifacts from empty-beam correction
sigmaHighpass = 100;
parfor holoIdx = 1:size(holograms,3)*size(holograms,4)
holo = holograms(:,:,holoIdx);
holo = mean(holo(:)) + (holo - imgaussfilt(holo, sigmaHighpass));
holograms(:,:,holoIdx) = holo;
end
%% OPTIONAL: removal of vertical and/or horizontal stripes in the holograms
% See documentation of removeStripes for details
stripeSettings = removeStripes;
stripeSettings.direction = 'vertical'; % direction of profile removal: 'vertical', 'horizontal', 'both'
stripeSettings.method = 'additive'; % method of stripe removal: 'additive' or 'multiplicative'
stripeSettings.rangeTop = 200; % range to average over at top of the image
stripeSettings.rangeBottom = 200; % range to average over at bottom of the image
stripeSettings.rangeLeft = 200; % range to average over at left edge of the image
stripeSettings.rangeRight = 200; % range to average over at right edge of the image
parfor holoIdx = 1:size(holograms,3)*size(holograms,4)
holo = holograms(:,:,holoIdx);
holo = removeStripes(holo, stripeSettings);
holograms(:,:,holoIdx) = holo;
end
%% OPTIONAL: save preprocessed holograms as raw (+ save setup parameters as mat)
save('-v7', [filePrefix, '_paramaters.mat'], 'thetas', 'fresnelNumbers');
savePrefix = [filePrefix, '_preprocessedHolograms'];
saveraw(holograms, savePrefix);
%% OPTIONAL: Load holograms - insert correct filenames, data type and size
parameterFile = 'parameterFile.mat';
hologramsFile = 'hologramsFile.raw';
hologramsDataType = 'single';
hologramsSize = [2048,2048,4,1501];
load(parameterFile);
clear 'holograms';
holograms = readraw(hologramsFile, hologramsDataType, hologramsSize);
%% OPTIONAL: find correct Fresnelnumber by comparing radially averaged PSF of first hologram to expected CTF: change 'scale' until minima coincide
fprintf('... check Fresnel number... \n' ) ;
projIdx = 1;
distIdx = 1;
scale = 1;
fresnelNumberExpected = scale * fresnelNumbers(1,distIdx);
settings.betaDeltaRatio = 0;
checkFresnelNumber(holograms(:,:,distIdx,projIdx), fresnelNumberExpected, settings);
%% PHASE RECONSTRUCTION
%% First determine suitable method and parameters for a random projection
projIdx = randi(size(holograms,4));
holos = double(holograms(:,:,:,projIdx));
N = [size(holos,1),size(holos,2)];
settingsRecon = struct;
settingsRecon.lim1 = 1e-3; % Regularization parameter for low frequencies
settingsRecon.lim2 = 3e-2; % Regularization parameter for high frequencies
settingsRecon.betaDeltaRatio = 1/100; % Estimated ratio between absorption and phase shifts
settingsRecon.padx = 128; % Horizontal and...
settingsRecon.pady = 128; % vertical padding: use value > 0 if strong artifacts occur at the image boundaries
%settingsRecon.maxPhase = 0; % Uncomment to impose negative phases
% Uncomment to impose support constraint via a binary mask (choice below =
% square support of half the aspect length of the full field-of-view)
%settingsRecon.support = padarray(true(N/2), N/4);
% CTF-reconstruction (linear)
phaseRecCTF = phaserec_ctf(holos, fresnelNumbers, settingsRecon);
figure('name', 'CTF-reconstruction'); showImage(phaseRecCTF);
% Nonlinear reconstruction: if the result differs significantly from
% the CTF-recon (if so, it looks hopefully better...), then it should be
% considered to reconstruct all holograms with the nonlinear approach
phaseRecNL = phaserec_nonlinTikhonov(holos, fresnelNumbers, settingsRecon);
figure('name', 'Nonlinear reconstruction'); showImage(phaseRecNL);
%% Reconstruct holograms for all angles with the chosen method and settings
method = 'phaserec_ctf'; % Replace by 'phaserec_nonlinTikhonov' to reconstruct nonlinearly
projs = zeros([size(holograms,1), size(holograms,2), size(holograms,4)], 'single');
for projIdx = 1:size(projs,3)
disp(projIdx);
projs(:,:,projIdx) = feval(method, double(holograms(:,:,:,projIdx)), fresnelNumbers, settingsRecon);
end
%% OPTIONAL: save projections as raw
savePrefix = [filePrefix, '_reconstructedProjections'];
saveraw(projs, savePrefix);
%% OPTIONAL: Load projections - insert correct filenames, data type and size
parameterFile = 'parameterFile.mat';
projectionsFile = 'projectionsFile.raw';
projectionsDataType = 'single';
hologramsSize = [2048,2048,1501];
load(parameterFile);
clear 'projs';
projs = readraw(projectionsFile, projectionsDataType, projectionsSize);
%% align rotation axis automatically
proj0Degrees = projs(:,:,1);
[~, idx180Degrees] = min(abs((thetas-thetas(1))-180));
proj180Degrees = projs(:,:,idx180Degrees);
% perform automatic alignment of shift and rotation (might fail for too large cone angles)
settings = alignImages;
settings.registrationMode = 'shiftRotation';
settings.registrationAccuracy = 20;
settings.sigmaLowpass = 0;
settings.medfilt = true;
settings.sigmaHighpass = 50; % High-pass filtering before alignment to avoid
% errors due to low-frequency artifacts
[~,shift,tiltAngle,~] = alignImages(fadeoutImage(proj0Degrees), fadeoutImage(fliplr(proj180Degrees)), settings);
shiftAxis = shift(2)/2;
detectorTilt = tiltAngle/2;
% check result
figure; showImage(shiftRotateImage(proj0Degrees, [0,shiftAxis], detectorTilt) - fliplr(shiftRotateImage(proj180Degrees, [0,shiftAxis], detectorTilt)))
title('Check that images cancel each other out! If not, automatic rotation-axis alignment has possibly failed.')
%% correct for position of rotation axis and detector tilt (only if known!! otherwise find out manually)
parfor indProj = 1:size(projs,3)
disp(indProj);
projs(:,:,indProj) = shiftRotateImage(projs(:,:,indProj), [0,shiftAxis], detectorTilt);
end
%% TOMOGRAPHIC RECONSTRUCTION (WITH OPTIONAL RING-REMOVAL)
%% First determine suitable reconstruction parameters by trial-reconstructions of a single tomographic slice
slicenumber = ceil(size(projs,1)/2);
sinoCrop = 50; % Number of pixels cropped from the edges of the sinogram before reconstruction
sinoPad = 500; % Number of pixels padded to the edges of the sinogram before reconstruction
sliceSize = 1800; % Size of the reconstructed slices
settingsTomoRec.outputSize = sliceSize;
% Ring-removal settings: see documentation of ringremove for details
doRingRemoval = false;
settingsRingremove = ringremove;
% Perform reconstruction with optional ring-removal
sino = squeeze(projs(slicenumber,:,:));
if doRingRemoval
sino = ringremove(sino, settingsRingremove);
end
slice = astraFBP(padarray(sino(sinoCrop+1:end-sinoCrop,:), [sinoPad,0], 'replicate'), thetas, settingsTomoRec);
% View results to see whether parameters have to be adjusted
figure; showImage(slice); colormap bone;
%% Reconstruct all slices with determined parameters
% optional ring removal for all slices
if doRingRemoval
parfor sliceIdx = 1:size(projs,1)
projs(sliceIdx,:,:) = ringremove(squeeze(projs(sliceIdx,:,:)), settingsRingremove);
end
end
slices = astraFBP(padarray(projs(:,sinoCrop+1:end-sinoCrop,:), [0, sinoPad], 'replicate'), thetas, settingsTomoRec);
%% show the reconstructed slices
fprintf('... show reconstructed slices... \n' ) ;
showStack(slices);
%% save results as raw
savePrefix = [filePrefix, '_reconstructedSlices'];
saveraw(slices, savePrefix);
%% OPTIONAL: filtering with Gaussian for better SNR
fprintf('... filter slices with a Gaussian function... \n' ) ;
filterSigma = 1; % standard deviation in pixels
slices = imgaussfilt3(slices,filterSigma);
%% save filtered results as raw
savePrefix = [filePrefix, '_reconstructedSlices_filtered_sigma=',num2str(filterSigma)];
saveraw(slices, savePrefix);
%% Start from scratch
clear 'variables'; clc;
set(0,'DefaultFigureColor','w');
%% Parameters of the scan
year = '2016';
runname = 'run56_LTP';
date = '20161116';
detector = 'pirra';
filePrefix = '12816-6_Kleinhirn3_tomo7';
%% (modify if loading fails)
if ispc
toolboxPath = 'S:/Projects_X-ray_Imaging/holotomotoolbox/functions';
astraPath = 'S:/Projects_X-ray_Imaging/ASTRAToolbox/windows';
prePath = 'S:/Messzeiten/';
elseif isunix
toolboxPath = '/IRP/AG_Salditt/Projects_X-ray_Imaging/holotomotoolbox/functions';
astraPath = '/IRP/AG_Salditt/Projects_X-ray_Imaging/ASTRAToolbox/linux';
prePath = '/IRP/AG_Salditt/Messzeiten/';
end
dataPath = fullfile(prePath, year, 'GINIX', runname, '/data/current/raw', date, 'detectors', detector, filePrefix);
if ~exist(dataPath,'dir')
error(['Assumed data path ''', dataPath, ''' does not exist. Scan parameters might be wrong. Please fix the path manually.']);
end
if ~exist(toolboxPath,'dir')
error(['Assigned path to the holotomotoolbox ''', toolboxPath, ''' does not exist.']);
end
if ~exist(astraPath,'dir')
error(['Assigned path to the ASTRA-toolbox ''', astraPath, ''' does not exist.']);
end
%% Add holotomotoolbox to path
addpath(genpath(toolboxPath));
%% Add ASTRA-toolbox to path
addpath(genpath(astraPath));
%% Set up image-reader for the specified detector
imageReader = getImageReaderIRP(detector);
switch(lower(detector))
case 'argos'
fileEnding = '_%04i.tif'; dumpEnding = '_%04i.tif.txt';
dx = 0.54e-6;
case 'talos'
fileEnding = '_%04i.tif'; dumpEnding = '_%04i.tif.txt';
dx = 75e-6;
case 'pirra'
fileEnding = '_%i.tif'; dumpEnding = '_%04i.tif.txt';
dx = 6.5e-6;
case 'thyia'
fileEnding = '_%i.tif'; dumpEnding = '_%04i.tif.txt';
dx = 4.5e-6;
case 'iris'
fileEnding = '_%04i.img'; dumpEnding = '_%04i.img.txt';
dx = 6.5e-6;
case 'timepix'
fileEnding = '_%04i.raw'; dumpEnding = '_%04i.raw.txt';
dx = 55e-6;
case 'zyla'
fileEnding = '_%04i.tif'; dumpEnding = '_%04i.tif.txt';
dx = 6.5e-6;
end
fileName = [filePrefix,fileEnding];
dumpName = [filePrefix,dumpEnding];
getFileName = @(number) fullfile(dataPath, sprintf(fileName,number));
getDumpName = @(number) fullfile(dataPath,'dump',sprintf(dumpName,number));