Skip to content
GitLab
Commit e091450c authored by smaretz's avatar smaretz
Browse files

minor stuff

parent eff03c5d
Hide whitespace changes
Inline Side-by-side
examples/template_tomo_analysis_Julia.m
View file @ e091450c
......@@ -2,6 +2,7 @@
clear all; close all; clc
set(0,'DefaultFigureColor','w');
%% parameters of the scan
spec = 'julia';
year = '2019';
......@@ -22,8 +23,8 @@ numAngles = numProjs+1;
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)
%% Paths for tomographic data and toolboxes to be loaded (modify if loading of toolboxes or data fails)
if ispc
toolboxPath = 'S:/Projects_X-ray_Imaging/holotomotoolbox/functions';
astraPath = 'S:/Projects_X-ray_Imaging/ASTRAToolbox/windows';
......@@ -33,7 +34,7 @@ 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/'];
saveDir = fullfile('/tmp/', filePrefix);
end
% change if not using julia
......@@ -90,6 +91,7 @@ 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(getFileName(number));
......@@ -97,6 +99,7 @@ raw0 = imageReader(getFileName(number));
showImage(raw0);
colormap gray
%% image numbers
clear emptynumbers;
......@@ -116,6 +119,7 @@ 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)
......@@ -171,13 +175,17 @@ parfor indAngle = 1:numAngles
thetas(indAngle) = wwas(getDumpName(number+(indImage-1)),rotMotor);
end
%% remove hotpixel errors
parfor indProj = 1:numAngles
projs(:,:,indProj) = removeOutliers(projs(:,:,indProj));
disp(indProj);
end
showImage(projs(:,:,1))
%% show random projection to check results
showImage(projs(:,:,randi(size(projs,3))));
%% align rotation axis automatically
proj0Degrees = projs(:,:,1);
......@@ -200,7 +208,7 @@ detectorTilt = tiltAngle/2;
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.')
title('Check that images cancel each other out in the center! If not, automatic rotation-axis alignment has possibly failed.')
%% if automatic alignment fails: pre-align rotation axis manually (for fine alignment use tomographic reconstruction)
%
......@@ -223,6 +231,7 @@ parfor indProj = 1:size(projs,3)
projs(:,:,indProj) = shiftRotateImage(projs(:,:,indProj), [0,shiftAxis], detectorTilt);
end
%% determine geometric parameters
% source-sample distance
z01 = wwas(getDumpName(numFirstangle),'sx') + z01Corr;
......@@ -251,12 +260,13 @@ 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);colormap gray
figure('name', 'Original image'); showImage(projExample); colormap gray;
settingsPhaseRecon = struct;
settingsPhaseRecon.padx = 500;
......@@ -265,12 +275,12 @@ settingsPhaseRecon.reg_alpha = 0.005;
% Reconstruction with MBA-algorithm
projRecMBA = phaserec_mba(projExample, F, settingsPhaseRecon);
figure('name', 'MBA-reconstruction'); showImage(projRecMBA);colormap gray
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);colormap gray
figure('name', 'BAC-reconstruction'); showImage(projRecBAC); colormap gray;
%% Reconstruct all projections angles with BAC using the determined settings
parfor indProj = 1:numAngles
......@@ -281,11 +291,12 @@ end
%% show stack of reconstructed projections
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
figure('name', 'Original slice'); showImage(sliceOriginal); colormap gray;
% additive approach for ring removal
settingsRingremove = ringremove;
......@@ -296,7 +307,7 @@ settingsRingremove.additive_GaussSigma = 1;
sinoAdditive = ringremove(sino,settingsRingremove);
sliceAdditive = astraFBP(sinoAdditive,thetas);
figure('name', 'Additive ring removal'); showImage(sliceAdditive);colormap gray
figure('name', 'Additive ring removal'); showImage(sliceAdditive); colormap gray;
% wavelet-based approach for ring removal
settingsRingremove.method = 'wavelet';
......@@ -306,15 +317,23 @@ settingsRingremove.wavelet_Sigma = 1;
sinoWavelet = ringremove(sino,settingsRingremove);
sliceWavelet = astraFBP(sinoWavelet,thetas);
figure('name', 'Wavelet-based ring removal'); showImage(sliceWavelet);colormap gray
figure('name', 'Wavelet-based ring removal'); showImage(sliceWavelet); colormap gray;
%% OPTIONAL: Perform ring-removal with determined method and settings (only reasonable if sliceAdditive or sliceWavelet looks cleaner than sliceOriginal)
settingsRingremove.method = 'additive';
parfor sliceIdx = 1:size(projs,1)
disp(sliceIdx);
projs(sliceIdx,:,:) = ringremove(squeeze(projs(sliceIdx,:,:)), settingsRingremove);
end
%% TOMOGRAPHIC RECONSTRUCTION (WITH OPTIONAL RING-REMOVAL)
%% First determine suitable reconstruction parameters by trial-reconstruction of a central tomographic slice
%% TOMOGRAPHIC RECONSTRUCTION
%% First determine suitable reconstruction parameters by trial-reconstruction of the central tomographic slice
slicenumber = ceil(size(projs,1)/2);
sinoCrop = 0; % Number of pixels cropped from the edges of the sinogram before reconstruction
sinoPad = ceil(0.25*size(projs,2)); % Number of pixels padded to the edges of the sinogram before reconstruction
sliceSize = size(projs,2); % Size of the reconstructed slices
sinoCrop = 0; % Number of pixels cropped from the edges of the sinogram before reconstruction
sinoPad = 0; % Number of pixels padded to the edges of the sinogram before reconstruction
sliceSize = size(projs,2); % Size of the reconstructed slices
settingsTomoRec = astraFDK;
settingsTomoRec.outputSize = sliceSize;
......@@ -323,53 +342,40 @@ settingsTomoRec.numSlices = 1; % reconstruct only one slice
settingsTomoRec.offset = 0; % position of the slice relative to the equator of the circular cone beam
% geometry. The value 0 corresponds to reconstructing the central slice
% Ring-removal settings: see documentation of ringremove for details
doRingRemoval = true;
settingsRingremove.method = 'additive';
% Perform reconstruction with optional ring-removal
sino = projs(slicenumber,:,:);
if doRingRemoval
sino(:) = ringremove(squeeze(sino), settingsRingremove);
end
sliceOriginal = astraFDK(padarray(sino(:,sinoCrop+1:end-sinoCrop,:), [0,sinoPad,0], 'replicate'), thetas, z01, z02, dx*1000 , settingsTomoRec);
sliceTest = astraFDK(padarray(projs(slicenumber,sinoCrop+1:end-sinoCrop,:), [0,sinoPad,0], 'replicate'), thetas, z01, z02, dx*1000, settingsTomoRec);
% View results to see whether parameters have to be adjusted
showImage(-sliceOriginal); colormap gray;
%% Optional ring removal for all slices with determined parameters
parfor sliceIdx = 1:size(projs,1)
disp(sliceIdx);
projs(sliceIdx,:,:) = ringremove(squeeze(projs(sliceIdx,:,:)), settingsRingremove);
end
showImage(sliceTest); colormap gray;
%% Tomographic reconstruction with the determined parameters
%% Tomographic reconstruction of all slices using the determined parameters
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);
saveraw(slices, fullfile(saveDir, 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);
saveraw(slices, fullfile(saveDir, savePrefix));
......
Supports Markdown
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