Commit 8e3c2a53 authored by p.jhagema's avatar p.jhagema
Browse files

Merge branch 'master' of gitlab.gwdg.de:irp/holotomotoolbox

parents bef31492 1fd83a4b
Pipeline #108649 passed with stage
in 1 minute and 54 seconds
...@@ -4,7 +4,6 @@ set(0,'DefaultFigureColor','w'); ...@@ -4,7 +4,6 @@ set(0,'DefaultFigureColor','w');
%% parameters of the scan %% parameters of the scan
spec = 'julia';
year = '2019'; year = '2019';
runname = '2019_05_02'; runname = '2019_05_02';
detector = 'argos'; detector = 'argos';
...@@ -29,16 +28,15 @@ if ispc ...@@ -29,16 +28,15 @@ if ispc
toolboxPath = 'S:/Projects_X-ray_Imaging/holotomotoolbox/functions'; toolboxPath = 'S:/Projects_X-ray_Imaging/holotomotoolbox/functions';
astraPath = 'S:/Projects_X-ray_Imaging/ASTRAToolbox/windows'; astraPath = 'S:/Projects_X-ray_Imaging/ASTRAToolbox/windows';
prePath = 'J:/'; prePath = 'J:/';
saveDir = fullfile('T:/Julialabor Tomo', runname, filePrefix, 'Slices');
elseif isunix elseif isunix
toolboxPath = '/IRP/AG_Salditt/Projects_X-ray_Imaging/holotomotoolbox/functions'; toolboxPath = '/home/AG_Salditt/Projects_X-ray_Imaging/holotomotoolbox/functions';
astraPath = '/IRP/AG_Salditt/Projects_X-ray_Imaging/ASTRAToolbox/linux'; astraPath = '/home/AG_Salditt/Projects_X-ray_Imaging/ASTRAToolbox/linux';
prePath = '/IRP/Labdata/AG_Salditt/julialab/'; prePath = '/home/Labdata/AG_Salditt/julialab/';
saveDir = fullfile('/tmp/', filePrefix);
end end
% change if not using julia % change if not using julia
dataPath = fullfile(prePath, year, 'julia' , runname, 'detectors', detector, filePrefix); dataPath = fullfile(prePath, year, 'julia' , runname, 'detectors', detector, filePrefix);
saveDir = fullfile(prePath, year, 'julia', runname, 'analysis', filePrefix, 'Slices');
if ~exist(toolboxPath,'dir') if ~exist(toolboxPath,'dir')
error(['Assigned path to the holotomotoolbox ''', toolboxPath, ''' does not exist.']); error(['Assigned path to the holotomotoolbox ''', toolboxPath, ''' does not exist.']);
...@@ -155,7 +153,7 @@ parfor indAngle = 1:numAngles ...@@ -155,7 +153,7 @@ parfor indAngle = 1:numAngles
number = rawnumbers(indAngle); number = rawnumbers(indAngle);
imgTmp = zeros(size(raw0,1),size(raw0,2),numAverage); imgTmp = zeros(size(raw0,1),size(raw0,2),numAverage);
for indImage = 1:numAverage for indImage = 1:numAverage
imgTmp(:,:,indImage) = imageReader(getFileName(number+(indImage-1))); imgTmp(:,:,indImage) = imageReader(getFileName(number+(indImage-1))); %#ok<PFBNS>
end end
% average over all images % average over all images
raw = median(imgTmp,3); raw = median(imgTmp,3);
...@@ -246,12 +244,7 @@ z12 = (z02-z01)/1000; ...@@ -246,12 +244,7 @@ z12 = (z02-z01)/1000;
z_eff = z12/M; z_eff = z12/M;
% main energy peak of the x-ray spectrum % main energy peak of the x-ray spectrum
switch lower(spec) E = 9.25; % gallium
case 'mona'
E = 17.5; % molybdenum
case 'julia'
E = 9.25; % gallium
end
% wavelength of the main energy peak % wavelength of the main energy peak
lambda = 12.398/E*1e-10; lambda = 12.398/E*1e-10;
......
function varargout = fftfreq(N, dx, computeMeshgrid) function varargout = fftfreq(N, dx)
% FFTFREQ creates a grid of Fourier-frequencies corresponding to the sampling points % FFTFREQ creates a grid of Fourier-frequencies corresponding to the sampling points
% of an n-dimensional FFT of given size. % of an n-dimensional FFT of given size.
% %
% ``varargout = fftfreq(N, dx, computeMeshgrid)`` % ``varargout = fftfreq(N, dx)``
% %
% %
% Parameters % Parameters
...@@ -12,16 +12,12 @@ function varargout = fftfreq(N, dx, computeMeshgrid) ...@@ -12,16 +12,12 @@ function varargout = fftfreq(N, dx, computeMeshgrid)
% dx : float or tuple of floats, optional % dx : float or tuple of floats, optional
% Underyling spacing in real-space, possibly different along each grid-dimension. % Underyling spacing in real-space, possibly different along each grid-dimension.
% Default = ones(size(N)). % Default = ones(size(N)).
% computeMeshgrid : tuple of non-negative integers, optional
% Return Fourier-frequency grid as meshgrid (case: true) or merely as coordinate vectors
% for the different dimensions? Default = false.
% %
% %
% Returns % Returns
% ------- % -------
% varargout : tuple of numerical arrays % varargout : tuple of numerical arrays
% numel(N) vectors of length N(1), N(2),... that span the Fourier-grid % numel(N) vectors of length N(1), N(2),... that span the Fourier-grid
% OR, if computeMeshgrid == true, the meshgrid corresponding to these vectors
% %
% See also % See also
% -------- % --------
...@@ -64,19 +60,18 @@ dx = dx(:).' .* ones([1,numel(N)]); ...@@ -64,19 +60,18 @@ dx = dx(:).' .* ones([1,numel(N)]);
ndim = numel(N); ndim = numel(N);
% Grid in Fourier space % Grid in Fourier space
xi = cell(1, ndim);
for jj = 1:ndim for jj = 1:ndim
xi{jj} = ifftshift( ( -floor(0.5*N(jj)) : floor(0.5*(N(jj)-1)) ).' * ( 2*pi / (N(jj)*dx(jj)) ) ); xi{jj} = ifftshift( ( -floor(0.5*N(jj)) : floor(0.5*(N(jj)-1)) ).' * ( 2*pi / (N(jj)*dx(jj)) ) );
end end
% Optionally assemble ndgrid (occupies more memory!) if ndim > 1
if nargin == 3 && computeMeshgrid
[xi{:}] = ndgrid(xi{:});
elseif ndim > 1
for jj = 1:ndim for jj = 1:ndim
xi{jj} = reshape( xi{jj}, [ones([1,jj-1]), N(jj), ones([1,ndim-jj])] ); xi{jj} = reshape( xi{jj}, [ones([1,jj-1]), N(jj), ones([1,ndim-jj])] );
end end
end end
% Unpack cell-array % Unpack cell-array
varargout = xi; varargout = xi;
......
...@@ -47,8 +47,8 @@ if nargin < 2 ...@@ -47,8 +47,8 @@ if nargin < 2
dx = ones(size(N)); dx = ones(size(N));
end end
xi = cell([numel(N),1]); xi = cell(numel(N),1);
[xi{:}] = fftfreq(N, dx, false); [xi{:}] = fftfreq(N, dx);
XiSq = xi{1}.^2; XiSq = xi{1}.^2;
for dim = 2:numel(xi) for dim = 2:numel(xi)
XiSq = XiSq + xi{dim}.^2; XiSq = XiSq + xi{dim}.^2;
......
...@@ -138,7 +138,7 @@ function kernel = fresnelPropagationKernel_fourier(N, fresnelNumbers, gpuFlag) ...@@ -138,7 +138,7 @@ function kernel = fresnelPropagationKernel_fourier(N, fresnelNumbers, gpuFlag)
ndim = numel(N); ndim = numel(N);
xi = cell([ndim,1]); xi = cell([ndim,1]);
[xi{:}] = fftfreq(N, 1, false); [xi{:}] = fftfreq(N, 1);
kernel = 1; kernel = 1;
for dim = 1:ndim for dim = 1:ndim
xi{dim} = gpuArrayIf(xi{dim}, gpuFlag); xi{dim} = gpuArrayIf(xi{dim}, gpuFlag);
...@@ -154,7 +154,7 @@ function kernel = fresnelPropagationKernel_chirp(N, fresnelNumbers, gpuFlag) ...@@ -154,7 +154,7 @@ function kernel = fresnelPropagationKernel_chirp(N, fresnelNumbers, gpuFlag)
ndim = numel(N); ndim = numel(N);
x = cell([ndim,1]); x = cell([ndim,1]);
[x{:}] = fftfreq(N, 2*pi./N, false); % fft-shifted grid with spacing 1 [x{:}] = fftfreq(N, 2*pi./N); % fft-shifted grid with spacing 1
kernel = prod(sqrt(abs(fresnelNumbers)) .* exp((-1i*pi/4) .* sign(fresnelNumbers))); kernel = prod(sqrt(abs(fresnelNumbers)) .* exp((-1i*pi/4) .* sign(fresnelNumbers)));
for dim = 1:ndim for dim = 1:ndim
x{dim} = gpuArrayIf(x{dim}, gpuFlag); x{dim} = gpuArrayIf(x{dim}, gpuFlag);
...@@ -170,7 +170,7 @@ function kernel = fresnelPropagationKernel_chirpLimited(N, fresnelNumbers, gpuFl ...@@ -170,7 +170,7 @@ function kernel = fresnelPropagationKernel_chirpLimited(N, fresnelNumbers, gpuFl
ndim = numel(N); ndim = numel(N);
x = cell([ndim,1]); x = cell([ndim,1]);
[x{:}] = fftfreq(N, 2*pi./N, false); % fft-shifted grid with spacing 1 [x{:}] = fftfreq(N, 2*pi./N); % fft-shifted grid with spacing 1
kernel = prod(sqrt(abs(fresnelNumbers)) .* exp((-1i*pi/4) .* sign(fresnelNumbers))); kernel = prod(sqrt(abs(fresnelNumbers)) .* exp((-1i*pi/4) .* sign(fresnelNumbers)));
for dim = 1:ndim for dim = 1:ndim
x{dim} = gpuArrayIf(x{dim}, gpuFlag); x{dim} = gpuArrayIf(x{dim}, gpuFlag);
...@@ -196,7 +196,7 @@ function kernel = fresnelPropagationKernel_pixel2pixel(N, fresnelNumbers, gpuFla ...@@ -196,7 +196,7 @@ function kernel = fresnelPropagationKernel_pixel2pixel(N, fresnelNumbers, gpuFla
ndim = numel(N); ndim = numel(N);
x = cell([ndim,1]); x = cell([ndim,1]);
[x{:}] = fftfreq(N, 2*pi./N, false); % fft-shifted grid with spacing 1 [x{:}] = fftfreq(N, 2*pi./N); % fft-shifted grid with spacing 1
kernel = 1; kernel = 1;
for dim = 1:ndim for dim = 1:ndim
x{dim} = gpuArrayIf(x{dim}, gpuFlag); x{dim} = gpuArrayIf(x{dim}, gpuFlag);
......
...@@ -28,8 +28,8 @@ function [averages, radii] = angularAverage(im) ...@@ -28,8 +28,8 @@ function [averages, radii] = angularAverage(im)
% %
% N = [1024,1024]; % N = [1024,1024];
% fresnel = 1e-2; % fresnel = 1e-2;
% [xi_1,xi_2] = fftfreq(N, 1, true); % xi2 = fftfreqNormSq(N, 1);
% ctf = fftshift(sin((xi_1.^2 + xi_2.^2)/(4*pi*fresnel))); % ctf = fftshift(sin(xi2/(4*pi*fresnel)));
% [ctfAveraged, radii] = angularAverage(ctf); % [ctfAveraged, radii] = angularAverage(ctf);
% %
% figure('name', 'image'); imagesc(ctf); colorbar; % figure('name', 'image'); imagesc(ctf); colorbar;
......
...@@ -76,9 +76,9 @@ end ...@@ -76,9 +76,9 @@ end
padsize = floor((outputSize - imageSize)/2); padsize = floor((outputSize - imageSize)/2);
padsize2 = outputSize - (imageSize + 2 * padsize); padsize2 = outputSize - (imageSize + 2 * padsize);
result = padarray(im, padsize, 'both', padval); result = padarray(im, padsize, padval, 'both');
if any(padsize2 > 0) if any(padsize2 > 0)
result = padarray(result, padsize2, 'pre', padval); result = padarray(result, padsize2, padval, 'pre');
end end
end end
...@@ -154,10 +154,10 @@ laplaceKernel = -fftfreqNormSq([ny, nx]); ...@@ -154,10 +154,10 @@ laplaceKernel = -fftfreqNormSq([ny, nx]);
switch (settings.reg_type) switch (settings.reg_type)
case 'none' case 'none'
% avoid division by 0 % avoid division by 0
laplaceKernel(0,0) = 1; laplaceKernel(1,1) = 1;
invLaplaceFilt = 1 ./ laplaceKernel; invLaplaceFilt = 1 ./ laplaceKernel;
% set zero frequency (singularity) to 0 % set zero frequency (singularity) to 0
invLaplaceFilt(0,0) = 0; invLaplaceFilt(1,1) = 0;
case 'mba' case 'mba'
invLaplaceFilt = -1 ./ (-laplaceKernel + settings.reg_alpha); invLaplaceFilt = -1 ./ (-laplaceKernel + settings.reg_alpha);
case 'tikhonov' case 'tikhonov'
......
...@@ -29,7 +29,7 @@ end ...@@ -29,7 +29,7 @@ end
ndim = numel(sizeHolo); ndim = numel(sizeHolo);
xi = cell([ndim,1]); xi = cell([ndim,1]);
[xi{:}] = fftfreq(sizeHolo, 1, false); [xi{:}] = fftfreq(sizeHolo, 1);
XiSqDivFresnel = 0; XiSqDivFresnel = 0;
for jj = 1:ndim for jj = 1:ndim
XiSqDivFresnel = XiSqDivFresnel + gpuArrayIf(xi{jj}, gpuFlag).^2 ./ (2*pi^2*fresnelNumbers(jj)); XiSqDivFresnel = XiSqDivFresnel + gpuArrayIf(xi{jj}, gpuFlag).^2 ./ (2*pi^2*fresnelNumbers(jj));
......
...@@ -6,8 +6,7 @@ function sinoCorrected = ringremove(sino,settings) ...@@ -6,8 +6,7 @@ function sinoCorrected = ringremove(sino,settings)
% Removes ring artifacts in tomographic reconstructions by reducing the amount of % Removes ring artifacts in tomographic reconstructions by reducing the amount of
% horizontal lines in the given sinogram. Ring removal is either based on the method % horizontal lines in the given sinogram. Ring removal is either based on the method
% proposed by Ketcham et al. :cite:`Ketcham_PS_2006`, denoted as 'additive', or on the method proposed % proposed by Ketcham et al. :cite:`Ketcham_PS_2006`, denoted as 'additive', or on the method proposed
% by Muench et al. :cite:`Muench_OE_2009`, denoted as 'wavelet'. Large parts of the code for the % by Muench et al. :cite:`Muench_OE_2009`, denoted as 'wavelet'.
% wavelet-based approach are taken from the original paper
% %
% Parameters % Parameters
% ---------- % ----------
...@@ -38,7 +37,7 @@ function sinoCorrected = ringremove(sino,settings) ...@@ -38,7 +37,7 @@ function sinoCorrected = ringremove(sino,settings)
% Highest decomposition level in the wavelet decomposition when using method % Highest decomposition level in the wavelet decomposition when using method
% 'wavelet'. % 'wavelet'.
% wavelet_Sigma : Default = 1 % wavelet_Sigma : Default = 1
% Standard deviation for the Gaussian filtering of the horizontal part of the % Width of the Gaussian filter kernel for the horizontal part of the
% sinogram in Fourier space when using method 'wavelet'. % sinogram in Fourier space when using method 'wavelet'.
% %
% Returns % Returns
...@@ -123,7 +122,6 @@ settings = completeStruct(settings, defaults); ...@@ -123,7 +122,6 @@ settings = completeStruct(settings, defaults);
switch settings.method switch settings.method
case 'additive' case 'additive'
numAngles = size(sino,2);
% averaging of sinogram along the angle % averaging of sinogram along the angle
sinoMean = mean(sino,2); sinoMean = mean(sino,2);
...@@ -140,50 +138,43 @@ switch settings.method ...@@ -140,50 +138,43 @@ switch settings.method
end end
% identification of horizontal stripes % identification of horizontal stripes
correctionProfile = sinoMean - sinoMeanSmoothed; correctionProfile = sinoMean - sinoMeanSmoothed;
% creation of 2d correction mask
correction2D = repmat(correctionProfile,[1 numAngles]);
% sinogram correction % sinogram correction
sinoCorrected = sino - settings.additive_Strength*correction2D; sinoCorrected = sino - settings.additive_Strength * correctionProfile;
case 'wavelet' case 'wavelet'
numAngles = size(sino,2);
sizeProjs = size(sino,1); % perform multilevel wavelet decomposition
H = cell(settings.wavelet_DecNum, 1);
% wavelet decomposition V = cell(settings.wavelet_DecNum, 1);
for decompIdx = 1:settings.wavelet_DecNum D = cell(settings.wavelet_DecNum, 1);
[sino,horzPart{decompIdx},vertPart{decompIdx},diagPart{decompIdx}] = ... S = cell(settings.wavelet_DecNum, 1);
dwt2(sino,settings.wavelet_Wname);
for l = 1:settings.wavelet_DecNum
S{l} = size(sino);
[sino, H{l}, V{l}, D{l}] = dwt2(sino,settings.wavelet_Wname);
end end
% FFT transform of horizontal frequency bands % filter horizontal frequency bands
for decompIdx = 1:settings.wavelet_DecNum for l = 1:settings.wavelet_DecNum
% FFT % FFT
horzFourier = fftshift(fft(horzPart{decompIdx}')); horzFourier = fft(H{l}');
[my,mx] = size(horzFourier);
% damping of horizontal stripe information % damping of horizontal stripe information
dampGauss = 1 - exp(-(-floor(my/2):-floor(my/2) + my - 1).^2/(2*settings.wavelet_Sigma^2)); m = size(horzFourier,1);
horzFourier = horzFourier.*repmat(dampGauss',1,mx); k = fftfreq(m, 2*pi/m);
dampGauss = 1 - exp(-1/2 * k.^2 / settings.wavelet_Sigma^2 );
horzFourier = horzFourier .* dampGauss;
% inverse FFT % inverse FFT
horzPart{decompIdx} = ifft(ifftshift(horzFourier))'; H{l} = ifft(horzFourier)';
end
% wavelet reconstruction
sinoCorrected = sino;
for decompIdx = settings.wavelet_DecNum:-1:1
sinoCorrected = sinoCorrected(1:size(vertPart{decompIdx},1),1:size(vertPart{decompIdx},2));
sinoCorrected = idwt2(sinoCorrected,horzPart{decompIdx},vertPart{decompIdx},...
diagPart{decompIdx},settings.wavelet_Wname);
end end
% make sure that corrected sinogram has the same size as the original sinogram % wavelet reconstruction
if size(sinoCorrected,2) == numAngles + 1 sinoCorrected = sino;
sinoCorrected = sinoCorrected(:,1:end - 1); for l = settings.wavelet_DecNum:-1:1
end sinoCorrected = sinoCorrected(1:size(V{l},1),1:size(V{l},2));
if size(sinoCorrected,1) == sizeProjs + 1 sinoCorrected = idwt2(sinoCorrected,H{l}, V{l}, D{l},settings.wavelet_Wname, S{l});
sinoCorrected = sinoCorrected(1:end - 1,:); end
end
otherwise otherwise
error('Choose between ringremove based on "additive" or "wavelet" correction for setting "method"'); error('Choose between ringremove based on "additive" or "wavelet" correction for setting "method"');
......
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