Commit 9c5700a4 authored by Simon Maretzke's avatar Simon Maretzke
Browse files

bisschen entschlackt... nur ein Vorschlag

parent 6a3ce55e
Pipeline #100446 passed with stage
in 1 minute and 50 seconds
......@@ -482,7 +482,7 @@ end
%% TOMOGRAPHIC RECONSTRUCTION (WITH OPTIONAL RING-REMOVAL)
%% First determine suitable reconstruction parameters by trial-reconstructions of a single tomographic slices
%% 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
......@@ -494,7 +494,6 @@ settingsTomoRec.outputSize = sliceSize;
doRingRemoval = false;
settingsRingremove = ringremove;
% Perform reconstruction with optional ring-removal
sino = squeeze(projs(slicenumber,:,:));
if doRingRemoval
......@@ -502,7 +501,6 @@ if doRingRemoval
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;
......@@ -520,7 +518,7 @@ slices = astraFBP(padarray(projs(:,sinoCrop+1:end-sinoCrop,:), [0, sinoPad], 're
%% show slices
%% show the reconstructed slices
fprintf('... show reconstructed slices... \n' ) ;
showStack(slices);
......
......@@ -2,7 +2,8 @@
clear all; close all; clc
set(0,'DefaultFigureColor','w');
%% details of the scan
%% parameters of the scan
spec = 'julia';
year = '2018';
runname = '2018_02_21';
......@@ -22,36 +23,65 @@ numAngles = numProjs+1;
z02 = 183.017;
z01Corr = 0; % in mm; correction factor for source-sample distance wrt to sx motor position
%% setup paths & toolbox
%% Add holotomotoolbox to path (modify if loading fails)
if ispc
prefix = 'J:/';
toolboxPath = 'S:/Projects_X-ray_Imaging/holotomotoolbox/functions';
elseif isunix
prefix = '/homegroups/Labdata/AG_Salditt/julialab/';
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/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 = [prefix,year,'/mona/'];
prePath = fullfile(prePath, year, '/mona/');
case 'julia'
prePath = [prefix,year ,'/julia/'];
otherwise
prePath = prefix;
prePath = fullfile(prePath, year, '/julia/');
end
% path to the detectors folder
prePath = [prePath,runname,'/detectors/'];
dataPath = fullfile(prePath, runname, '/detectors/', detector, filePrefix);
% path to the toolbox folder
tbPath = 'H:/Toolboxes/holotomotoolbox';
addpath(genpath(tbPath));
if ~exist(dataPath,'dir')
error(['Assumed data path ''', dataPath, ''' does not exist. Scan parameters might be wrong. Please fix the path manually.']);
end
% path to the ASTRA toolbox
astraPath = 'H:/Toolboxes/toolbox/astra-toolbox-FDK-windows';
addpath(genpath(astraPath))
%% set up detector
%% Set up image-reader for the specified detector
imageReader = getImageReaderIRP(detector);
dataPath = [prePath,detector,'/',filePrefix];
switch(lower(detector))
case 'argos'
......@@ -80,13 +110,15 @@ end
fileName = [dataPath,'/',filePrefix,fileEnding];
dumpName = [dataPath,'/dump/',filePrefix,dumpEnding];
%% view one projection to check that everything is correct
%% Read one image to check if path is correct
number = numFlats+1;
raw0 = imageReader(sprintf(fileName,number));
showImage(raw0);
colormap gray
%% image numbers
clear emptynumbers;
......@@ -119,6 +151,7 @@ for indEmpty = 1:numel(emptynumbers)
fprintf('\n');
end
%% read in darkfields
darkTmp = zeros(size(raw0,1),size(raw0,2),numDarks);
for indImg = 1:numDarks
......@@ -129,6 +162,7 @@ end
dark = median(darkTmp,3);
fprintf('\n');
%% read in tomographic scan data
binningFactor = 1; % factor for binning of detector pixels, factor 1 does nothing
projs = zeros(size(raw0,1)/binningFactor,size(raw0,2)/binningFactor,numAngles,'single');
......@@ -161,42 +195,39 @@ for indAngle = 1:numAngles
thetas(indAngle) = wwas(sprintf(dumpName,number+(indImage-1)),rotMotor);
end
%% remove hotpixel errors
projsFilt = projs;
parfor indProj=1:numAngles
projsFilt(:,:,indProj) = removeOutliers(projs(:,:,indProj));
parfor indProj = 1:numAngles
projs(:,:,indProj) = removeOutliers(projs(:,:,indProj));
disp(indProj);
end
showImage(projsFilt(:,:,1))
showImage(projs(:,:,1))
%% align rotation axis automatically
firstProj = projsFilt(:,:,1);
lastProj = projsFilt(:,:,round(180/(angRange/numProjs)+1));
% cut away edges that might prevent proper image registration
cutTopBottom = 0;
cutLeftRight = 0;
%% 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 = registerImages;
settings = alignImages;
settings.registrationMode = 'shiftRotation';
settings.registrationAccuracy = 100;
settings.maxIterations = 10;
[shift,tiltAngle] = registerImages(firstProj(1+cutTopBottom:end-cutTopBottom,1+cutLeftRight:end-cutLeftRight),...
fliplr(lastProj(1+cutTopBottom:end-cutTopBottom,1+cutLeftRight:end-cutLeftRight)),settings);
shiftAxis = -shift(2)/2;
detectorTilt = -tiltAngle/2;
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
showImage(shiftRotateImage(firstProj,-shift/2,-tiltAngle/2) - fliplr(shiftRotateImage(lastProj,-shift/2,-tiltAngle/2)))
title('Check that images cancel each other out! If not: repeat procedure with more iterations or different registrationAccuracy!')
caxis([-0.4 0.4])
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.')
%% if automatic alignment fails: pre-align rotation axis manually (for fine alignment use tomographic reconstruction)
% firstProj = projsFilt(:,:,1);
% lastProj = projsFilt(:,:,round(180/(angRange/numProjs)+1));
%% if automatic alignment fails: pre-align rotation axis manually (for fine alignment use tomographic reconstruction)
%
% % manually determined parameters
% shiftAxis = 0;
......@@ -208,9 +239,18 @@ caxis([-0.4 0.4])
% settings.tiltAngle = detectorTilt;
%
% % check result
% showImage(checkRotaxis(firstProj,lastProj,settings));
% 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;
......@@ -237,224 +277,125 @@ lambda = 12.398/E*1e-10;
% Fresnel number
F = dx_eff^2/(lambda*z_eff);
fprintf('Fresnel number: %4.5f\n',F)
fprintf('Fresnel number: %4.5f\n',F);
%% find parameters for phase retrieval - starting with MBA step
Itmp = projsFilt(:,:,1);
% settings for phase retrieval
padx = 500;
pady = 500;
regAlpha = 0.01;
mbasettings = phaserec_mba;
mbasettings.padx = padx;
mbasettings.pady = pady;
mbasettings.reg_alpha = regAlpha;
%% PHASE RECONSTRUCTION
%% First determine suitable method and parameters for an exemplary projection
projIdx = 1; %randi(size(projs,3));
projExample = projs(:,:,:,projIdx);
% reconstruction via MBA to remove edge enhancement
phi = phaserec_mba(Itmp,F,mbasettings);
figure('name', 'Original image'); showImage(projExample);
% show result
showImage(phi)
settingsRecon = struct;
settingsRecon.padx = 500;
settingsRecon.pady = 500;
settingsRecon.reg_alpha = 0.01;
%% sharpen the image via the BAC algorithm
regGamma = 1; % relative to 1/(2*pi)
% Reconstruction with MBA-algorithm
projRecMBA = phaserec_mba(projExample, F, settingsRecon);
figure('name', 'MBA-reconstruction'); showImage(projRecMBA);
% Parameter
bacsettings = phaserec_bac;
bacsettings.padx = padx;
bacsettings.pady = pady;
bacsettings.reg_alpha = regAlpha;
bacsettings.reg_gamma = regGamma;
% Reconstruction with BAC-algorithm
settingsRecon.reg_gamma = 1;
projRecBAC = phaserec_bac(projExample, settingsRecon);
figure('name', 'BAC-reconstruction'); showImage(projRecBAC);
% reconstruction of intensity in the object plane
ampl = phaserec_bac(Itmp,bacsettings);
% show result
showImage(ampl)
%% compare to original intensity
showImage(Itmp)
%% reconstruct all projections
projsRec = zeros(size(projs),'single');
%% Reconstruct all projections angles with BAC using the determined settings
parfor indProj = 1:numAngles
Itmp = projsFilt(:,:,indProj);
% Parameter
bacsettings = struct;
bacsettings.padx = padx;
bacsettings.pady = pady;
bacsettings.reg_alpha = regAlpha;
bacsettings.reg_gamma = regGamma;
ampl = phaserec_bac(Itmp,bacsettings );
projsRec(:,:,indProj) = ampl;
fprintf('.');
if(rand < 0.02)
fprintf('\n');
end
disp(indProj);
projs(:,:,indProj) = phaserec_bac(projs(:,:,indProj), bacsettings);
end
%% show stack of reconstructed projections
settings = showStack;
settings.step = 10;
settings.caxis = -1;
showStack(projsRec,settings);
%% reconstruct single slice to inspect severeness of ring artifacts
sino = squeeze(projsRec(ceil(size(projsRec,1)/2),:,:));
sino = circshiftSubPixel(sino,[shiftAxis 0]);
slice = iradon(sino,thetas);
showImage(slice);
caxis(ans)
% zoom(2);
%% additive approach for ring removal
settingsAdditive = ringremove;
settingsAdditive.method = 'additive';
settingsAdditive.additive_SmoothMethod = 'mean';
settingsAdditive.additive_WindowSize = 20;
settingsAdditive.additive_GaussSigma = 5;
sinoAdditive = ringremove(sino,settingsAdditive);
sliceAdditive = iradon(sinoAdditive,thetas);
showImage(sliceAdditive);
caxis(ans)
% zoom(2);
%% wavelet-based approach for ring removal
settingsWavelet = ringremove;
settingsWavelet.method = 'wavelet';
settingsWavelet.wavelet_DecNum = 3;
settingsWavelet.wavelet_Wname = 'sym8';
settingsWavelet.wavelet_Sigma = 1;
sinoWavelet = ringremove(sino,settingsWavelet);
sliceWavelet = iradon(sinoWavelet,thetas);
showImage(sliceWavelet);
caxis(ans)
% zoom(2)
%% remove rings in all projections
% choose your favourite approach
algorithmToUse = 1; % 1: additive approach, 2: wavelet-based approach
projsClean = projsRec;
if algorithmToUse == 1
parfor indProj = 1:size(projsClean,1)
sino = squeeze(projsRec(indProj,:,:));
projsClean(indProj,:,:) = ringremove(sino,settingsAdditive);
end
showStack(projs,settings);
%% Take logarithm to convert intensities to absorption values
projsPad = -log(projs);
%% TOMOGRAPHIC RECONSTRUCTION (WITH OPTIONAL RING-REMOVAL)
%% First determine suitable reconstruction parameters by trial-reconstruction of a central 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 = astraFDK;
settingsTomoRec.outputSize = sliceSize;
settingsTomoRec.shortScan = 1; % short scan should be enabled when the scan range is smaller than 360 degree
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 = false;
settingsRingremove = ringremove;
% Perform reconstruction with optional ring-removal
sino = squeeze(projs(slicenumber,:,:));
if doRingRemoval
sino = ringremove(sino, settingsRingremove);
end
if algorithmToUse == 2
parfor indProj = 1:size(projsClean,1)
sino = squeeze(projsRec(indProj,:,:));
projsClean(indProj,:,:) = ringremove(sino,settingsWavelet);
slice = astraFDK(padarray(sino(sinoCrop+1:end-sinoCrop,:), [sinoPad,0], 'replicate'), thetas, z01, z02, dx*1000 , 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
%% pad the projections to reduce region-of-interest tomography artifacts
padval = 500;
projsPad = padarray(-log(projsClean),[0 padval 0],'replicate');
settingsTomoRec.offset = 0;
settingsTomoRec.numSlices = size(projs,1); % all slices
slices = astraFBP(padarray(projs(:,sinoCrop+1:end-sinoCrop,:), [0, sinoPad], 'replicate'), thetas, settingsTomoRec);
%% correct for position of rotation axis and detector tilt (only if known!! otherwise go to end of script and find out manually)
projsAstra = projsPad;
for indProj = 1:size(projsPad,3)
projsAstra(:,:,indProj) = shiftRotateImage(projsPad(:,:,indProj),[0 shiftAxis],detectorTilt);
end
%% reconstruction with the ASTRA toolbox CUDA implementation of the FDK algorithm
settingsAstra = astraFDK;
% size of the reconstructed slices
settingsAstra.outputSize = 1500;
%% show the reconstructed slices
fprintf('... show reconstructed slices... \n' ) ;
showStack(slices);
% short scan should be enabled when the scan range is smaller than 360 degree
settingsAstra.shortScan = 1;
% number of reconstructed slices
settingsAstra.numSlices = 1;
% settingsAstra.numSlices=size(projsAstra,1); % all slices
%% save results as raw
savePrefix = [filePrefix, '_reconstructedSlices'];
saveraw(slices, savePrefix);
% if you want to reconstruct the whole volume set offset to 0, otherwise it
% indicates the position relative to the central line
settingsAstra.offset = -500;
if(settingsAstra.numSlices>1 && settingsAstra.offset ~=0)
reply = input(sprintf('You are about to reconstruct %i slices. Do you want to keep the slice offset of %i? Y/N [Y]:',...
settingsAstra.numSlices,settingsAstra.offset),'s');
if isempty(reply)
reply = 'Y';
end
if strcmpi(reply,'N')
settingsAstra.offset = input('New offset: ');
end
end
% do the reconstruction
tic
slicesAstra = astraFDK(projsAstra,thetas,z01,z02,dx*1000,settingsAstra);
toc
% show central slice
showImage(-slicesAstra(:,:,ceil(size(slicesAstra,3)/2)))
% zoom(2)
%% show slices
stepSize = 20; % step size between displayed slices
sliceStart = 1; % first slice to be shown
sliceEnd = 1200; % last slice to be shown (maximum: min(size(slicesAstra))
for indSlice = sliceStart:stepSize:sliceEnd
settings = showOrthoslices;
settings.cmap = flipud(gray);
settings.caxis = [-0.4 0.7]*1e-3;
settings.figIdx = 1;
settings.sliceIndices = [indSlice indSlice indSlice];
showOrthoslices(slicesFiltered,settings);
drawnow
end
%% 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);
%% lowpass filtering to increase SNR (CAUTION: increases PSF of the imaging system)
gaussSigma = 1;
slicesFiltered = imgaussfilt3(slicesAstra,gaussSigma);
%% save the unfiltered slices
saveDir = ['T:\Julialabor Tomo\',runname,'\',filePrefix,'\Slices\'];
if exist(saveDir,'dir') ~= 7
mkdir(saveDir)
end
%% save filtered results as raw
savePrefix = [filePrefix, '_reconstructedSlices_filtered_sigma=',num2str(filterSigma)];
saveraw(slices, savePrefix);
cd(saveDir)
tic
saveraw(slicesAstra,filePrefix);
toc
%% save the filtered slices
saveDir = ['T:\Julialabor Tomo\',runname,'\',filePrefix,'\Slices\'];
if exist(saveDir,'dir') ~= 7
mkdir(saveDir)
end
cd(saveDir)
tic
saveraw(slicesFiltered,[filePrefix,'_filtered_sigma=',num2str(gaussSigma)]);
toc
%% %%%%%%%%%%%%%%%%%%%%% optional: manually determine shift and tilt of rotation axis %%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%% OPTIONAL: manually determine shift and tilt of rotation axis %%%%%%%%%%%%%%%%%%%%%%%
%% align rotation axis manually (for fine alignment use tomographic reconstruction)
firstProj = projsFilt(:,:,1);
lastProj = projsFilt(:,:,round(180/(angRange/numProjs)+1));
proj0Degrees = projs(:,:,1);
[~, idx180Degrees] = min(abs((thetas-thetas(1))-180));
proj180Degrees = projs(:,:,idx180Degrees);
% manually determined parameters
shiftx = 0;
......@@ -464,7 +405,7 @@ settings.shiftx = shiftx;
settings.shifty = shifty;
% check result
showImage(checkRotaxis(firstProj,lastProj,settings));
showImage(checkRotaxis(proj0Degrees,proj180Degrees,settings));
caxis([-0.4 0.4])
%% align rotation axis via tomographic reconstruction
......
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