Skip to content

GitLab

Commit 80035f27 authored by mtoeppe's avatar mtoeppe
Browse files

update templates

parent 08b463f9
Hide whitespace changes
Inline Side-by-side
examples/template_tomo_analysis_GINIX.m
View file @ 80035f27
......@@ -13,48 +13,37 @@ filePrefix = '12816-6_Kleinhirn3_tomo7';
%% Add holotomotoolbox to path (modify if loading fails)
%% (modify if loading fails)
if ispc
toolboxPath = 'S:/Projects_X-ray_Imaging/holotomotoolbox/functions';
elseif isunix
toolboxPath = '/IRP/AG_Salditt/Projects_X-ray_Imaging/holotomotoolbox/functions';
end
if ~exist(toolboxPath,'dir')
error(['Assigned path to the holotomotoolbox ''', toolboxPath, ''' does not exist.']);
end
addpath(genpath(toolboxPath));
%% Add ASTRA-toolbox to path
if ispc
astraPath = 'S:/Projects_X-ray_Imaging/ASTRAToolbox/windows';
elseif isunix
astraPath = '/IRP/AG_Salditt/Projects_X-ray_Imaging/ASTRAToolbox/linux';
end
if ~exist(astraPath,'dir')
error(['Assigned path to the ASTRA-toolbox ''', astraPath, ''' does not exist.']);
end
addpath(genpath(astraPath));
%% Set up data paths
if ispc
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);
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);
......
examples/template_tomo_analysis_Julia.m 100644 → 100755
View file @ 80035f27
......@@ -2,16 +2,15 @@
clear all; close all; clc
set(0,'DefaultFigureColor','w');
%% parameters of the scan
spec = 'julia';
year = '2018';
runname = '2018_02_21';
year = '2019';
runname = '2019_05_02';
detector = 'argos';
filePrefix = 'Hirn7_rechts_tomo1';
filePrefix = 'HH_TAB1_C3-3_tomo1_unten';
rotMotor = 'szrot';
angRange = 180;
angRange = 181;
numProjs = 1000;
numAverage = 1;
......@@ -20,72 +19,50 @@ numDarks = 10;
numAngles = numProjs+1;
%% source-detector distance and correction factor for source-sample distance
z02 = 183.017;
z02 = 185.74;
z01Corr = 0; % in mm; correction factor for source-sample distance wrt to sx motor position
%% Add holotomotoolbox to path (modify if loading fails)
if ispc
toolboxPath = 'S:/Projects_X-ray_Imaging/holotomotoolbox/functions';
astraPath = 'S:/Projects_X-ray_Imaging/ASTRAToolbox/windows';
prePath = 'J:/';
saveDir = fullfile('T:/Julialabor Tomo', runname, filePrefix, 'Slices');
elseif isunix
toolboxPath = '/IRP/AG_Salditt/Projects_X-ray_Imaging/holotomotoolbox/functions';
astraPath = '/IRP/AG_Salditt/Projects_X-ray_Imaging/ASTRAToolbox/linux';
prePath = '/IRP/Labdata/AG_Salditt/julialab/';
saveDir = ['/tmp/'];
end
% change if not using julia
dataPath = fullfile(prePath, year, 'julia' , runname, 'detectors', detector, filePrefix);
if ~exist(toolboxPath,'dir')
error(['Assigned path to the holotomotoolbox ''', toolboxPath, ''' does not exist.']);
end
addpath(genpath(toolboxPath));
%% Add ASTRA-toolbox to path
if ispc
astraPath = 'S:/Projects_X-ray_Imaging/holotomotoolbox/ASTRAToolbox/windows';
elseif isunix
computerName = char(java.net.InetAddress.getLocalHost.getHostName);
% If name is of the form 'santa0.roentgen.physik.uni-goettingen.de',
% only take prefix until first point:
pointIdx = strfind(computerName, '.');
if ~isempty(pointIdx)
computerName = computerName(1:pointIdx(1)-1);
end
astraPath = ['/IRP/AG_Salditt/Projects_X-ray_Imaging/ASTRAToolbox/', computerName];
end
if ~exist(astraPath,'dir')
error(['Assigned path to the ASTRA-toolbox ''', astraPath, ''' does not exist.']);
end
addpath(genpath(astraPath));
%% Set up data paths
if ispc
prePath = 'J:/';
elseif isunix
prePath = '/IRP/Labdata/AG_Salditt/julialab/';
end
switch lower(spec)
case 'mona'
prePath = fullfile(prePath, year, '/mona/');
case 'julia'
prePath = fullfile(prePath, year, '/julia/');
end
dataPath = fullfile(prePath, runname, '/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
%% 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.edf.txt';
fileEnding = '_%04i.tif'; dumpEnding = '_%04i.tif.txt';
dx = 0.54e-6;
case 'talos'
fileEnding = '_%04i.tif'; dumpEnding = '_%04i.tif.txt';
......@@ -107,18 +84,19 @@ switch(lower(detector))
dx = 6.5e-6;
end
fileName = [dataPath,'/',filePrefix,fileEnding];
dumpName = [dataPath,'/dump/',filePrefix,dumpEnding];
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 = numFlats+1;
raw0 = imageReader(sprintf(fileName,number));
raw0 = imageReader(getFileName(number));
showImage(raw0);
colormap gray
%% image numbers
clear emptynumbers;
......@@ -138,13 +116,12 @@ emptynumbers{2} = numLastangle+1:numLastangle+numFlats;
% darkfield images after tomographic scan
darknumbers = 2*numFlats+numAngles*numAverage+1:2*numFlats+numAngles*numAverage+numDarks;
%% read in empty images
emptys = cell(numel(emptynumbers),1);
for indEmpty = 1:numel(emptynumbers)
emptyTmp = zeros(size(raw0,1),size(raw0,2),numFlats);
for indImg = 1:numFlats
emptyTmp(:,:,indImg) = imageReader(sprintf(fileName,emptynumbers{indEmpty}(indImg)));
emptyTmp(:,:,indImg) = imageReader(getFileName(emptynumbers{indEmpty}(indImg)));
fprintf('.');
end
% take median of all emptys
......@@ -152,31 +129,29 @@ for indEmpty = 1:numel(emptynumbers)
fprintf('\n');
end
%% read in darkfields
darkTmp = zeros(size(raw0,1),size(raw0,2),numDarks);
for indImg = 1:numDarks
darkTmp(:,:,indImg) = imageReader(sprintf(fileName,darknumbers(indImg)));
darkTmp(:,:,indImg) = imageReader(getFileName(darknumbers(indImg)));
fprintf('.');
end
% take median of all darkfield images
dark = median(darkTmp,3);
fprintf('\n');
%% read in tomographic scan data
binningFactor = 2; % factor for binning of detector pixels, factor 1 does nothing
projs = zeros(size(raw0,1)/binningFactor,size(raw0,2)/binningFactor,numAngles,'single');
thetas = zeros(numAngles,1);
for indAngle = 1:numAngles
parfor indAngle = 1:numAngles
fprintf('Reading angle %4i of %4i\n',indAngle,numAngles);
% read images and tka emdian over all images
% read images and tka median over all images
number = rawnumbers(indAngle);
imgTmp = zeros(size(raw0,1),size(raw0,2),numAverage);
for indImage = 1:numAverage
imgTmp(:,:,indImage) = imageReader(sprintf(fileName,number+(indImage-1)));
imgTmp(:,:,indImage) = imageReader(getFileName(number+(indImage-1)));
end
% average over all images
raw = median(imgTmp,3);
......@@ -193,10 +168,9 @@ for indAngle = 1:numAngles
end
% read tomographic angle from motor position
thetas(indAngle) = wwas(sprintf(dumpName,number+(indImage-1)),rotMotor);
thetas(indAngle) = wwas(getDumpName(number+(indImage-1)),rotMotor);
end
%% remove hotpixel errors
parfor indProj = 1:numAngles
projs(:,:,indProj) = removeOutliers(projs(:,:,indProj));
......@@ -205,7 +179,6 @@ end
showImage(projs(:,:,1))
%% align rotation axis automatically
proj0Degrees = projs(:,:,1);
[~, idx180Degrees] = min(abs((thetas-thetas(1))-180));
......@@ -224,10 +197,11 @@ shiftAxis = shift(2)/2;
detectorTilt = tiltAngle/2;
% check result
figure; showImage(shiftRotateImage(proj0Degrees, [0,shiftAxis], detectorTilt) - fliplr(shiftRotateImage(proj180Degrees, [0,shiftAxis], detectorTilt)))
showImage(shiftRotateImage(proj0Degrees, [0,shiftAxis], detectorTilt) - fliplr(shiftRotateImage(proj180Degrees, [0,shiftAxis], detectorTilt)))
colormap gray
caxis([-0.4 0.4])
title('Check that images cancel each other out! If not, automatic rotation-axis alignment has possibly failed.')
%% if automatic alignment fails: pre-align rotation axis manually (for fine alignment use tomographic reconstruction)
%
% % manually determined parameters
......@@ -243,18 +217,15 @@ title('Check that images cancel each other out! If not, automatic rotation-axis
% showImage(checkRotaxis(proj0Degrees,proj180Degrees,settings));
% caxis([-0.4 0.4])
%% correct for position of rotation axis and detector tilt (only if known!! otherwise find out manually, see bottom of this script)
parfor indProj = 1:size(projs,3)
disp(indProj);
projs(:,:,indProj) = shiftRotateImage(projs(:,:,indProj), [0,shiftAxis], detectorTilt);
end
%% determine geometric parameters
% source-sample distance
z01 = wwas(sprintf(dumpName,numFirstangle),'sx') + z01Corr;
z01 = wwas(getDumpName(numFirstangle),'sx') + z01Corr;
% effective pixel size
M = z02/z01;
......@@ -280,29 +251,26 @@ lambda = 12.398/E*1e-10;
F = dx_eff^2/(lambda*z_eff);
fprintf('Fresnel number: %4.5f\n',F);
%% PHASE RECONSTRUCTION
%% First determine suitable method and parameters for an exemplary projection
projIdx = 1; %randi(size(projs,3));
projExample = projs(:,:,projIdx);
figure('name', 'Original image'); showImage(projExample);
figure('name', 'Original image'); showImage(projExample);colormap gray
settingsPhaseRecon = struct;
settingsPhaseRecon.padx = 500;
settingsPhaseRecon.pady = 500;
settingsPhaseRecon.reg_alpha = 0.01;
settingsPhaseRecon.reg_alpha = 0.005;
% Reconstruction with MBA-algorithm
projRecMBA = phaserec_mba(projExample, F, settingsPhaseRecon);
figure('name', 'MBA-reconstruction'); showImage(projRecMBA);
figure('name', 'MBA-reconstruction'); showImage(projRecMBA);colormap gray
% Reconstruction with BAC-algorithm
settingsPhaseRecon.reg_gamma = 1;
projRecBAC = phaserec_bac(projExample, settingsPhaseRecon);
figure('name', 'BAC-reconstruction'); showImage(projRecBAC);
figure('name', 'BAC-reconstruction'); showImage(projRecBAC);colormap gray
%% Reconstruct all projections angles with BAC using the determined settings
parfor indProj = 1:numAngles
......@@ -310,16 +278,35 @@ parfor indProj = 1:numAngles
projs(:,:,indProj) = phaserec_bac(projs(:,:,indProj), settingsPhaseRecon);
end
%% show stack of reconstructed projections
showStack(projs,settings);
showStack(projs(11:end-10,11:end-10,:));
%% Determine settings for ringremoval
sino = squeeze(projs(ceil(size(projs,1)/2),:,:));
sliceOriginal = astraFBP(sino,thetas);
figure('name', 'Original slice'); showImage(sliceOriginal);colormap gray
% additive approach for ring removal
settingsRingremove = ringremove;
settingsRingremove.method = 'additive';
settingsRingremove.additive_SmoothMethod = 'mean';
settingsRingremove.additive_WindowSize = 7;
settingsRingremove.additive_GaussSigma = 1;
%% Take logarithm to convert intensities to absorption values
projs = -log(projs);
sinoAdditive = ringremove(sino,settingsRingremove);
sliceAdditive = astraFBP(sinoAdditive,thetas);
figure('name', 'Additive ring removal'); showImage(sliceAdditive);colormap gray
% wavelet-based approach for ring removal
settingsRingremove.method = 'wavelet';
settingsRingremove.wavelet_DecNum = 3;
settingsRingremove.wavelet_Wname = 'sym8';
settingsRingremove.wavelet_Sigma = 1;
sinoWavelet = ringremove(sino,settingsRingremove);
sliceWavelet = astraFBP(sinoWavelet,thetas);
figure('name', 'Wavelet-based ring removal'); showImage(sliceWavelet);colormap gray
%% TOMOGRAPHIC RECONSTRUCTION (WITH OPTIONAL RING-REMOVAL)
%% First determine suitable reconstruction parameters by trial-reconstruction of a central tomographic slice
......@@ -338,18 +325,17 @@ settingsTomoRec.offset = 0; % position of the slice relative to the equator
% Ring-removal settings: see documentation of ringremove for details
doRingRemoval = true;
settingsRingremove = ringremove;
settingsRingremove.method = 'additive';
% Perform reconstruction with optional ring-removal
sino = projs(slicenumber,:,:);
if doRingRemoval
sino(:) = ringremove(squeeze(sino), settingsRingremove);
end
slice = astraFDK(padarray(sino(:,sinoCrop+1:end-sinoCrop,:), [0,sinoPad,0], 'replicate'), thetas, z01, z02, dx*1000 , settingsTomoRec);
sliceOriginal = astraFDK(padarray(sino(:,sinoCrop+1:end-sinoCrop,:), [0,sinoPad,0], 'replicate'), thetas, z01, z02, dx*1000 , settingsTomoRec);
% View results to see whether parameters have to be adjusted
figure; showImage(slice); colormap bone;
showImage(-sliceOriginal); colormap gray;
%% Reconstruct all slices with determined parameters
......@@ -365,34 +351,30 @@ settingsTomoRec.offset = 0;
settingsTomoRec.numSlices = size(projs,1); % all slices
slices = astraFDK(padarray(projs(:,sinoCrop+1:end-sinoCrop,:), [0,sinoPad,0], 'replicate'), thetas, z01, z02, dx*1000 , settingsTomoRec);
%% show the reconstructed slices
fprintf('... show reconstructed slices... \n' ) ;
showStack(slices);
%% save results as raw
if exist(saveDir,'dir') ~= 7
mkdir(saveDir)
end
cd(saveDir)
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);
%% %%%%%%%%%%%%%%%%%%%%% OPTIONAL: manually determine shift and tilt of rotation axis %%%%%%%%%%%%%%%%%%%%%%%
%% align rotation axis manually (for fine alignment use tomographic reconstruction)
proj0Degrees = projs(:,:,1);
......
Markdown is supported
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