Commit fbf3d249 authored by j.schaeper's avatar j.schaeper
Browse files

Merge branch 'master' of gwdg-gitlab-stud:irp/holotomotoolbox

parents e355c173 0ee2bd84
Pipeline #103308 passed with stage
in 3 minutes and 17 seconds
function varargout = fftfreq(N, dx, computeMeshgrid)
function varargout = fftfreq(N, dx)
% FFTFREQ creates a grid of Fourier-frequencies corresponding to the sampling points
% of an n-dimensional FFT of given size.
%
% ``varargout = fftfreq(N, dx, computeMeshgrid)``
% ``varargout = fftfreq(N, dx)``
%
%
% Parameters
......@@ -12,16 +12,12 @@ function varargout = fftfreq(N, dx, computeMeshgrid)
% dx : float or tuple of floats, optional
% Underyling spacing in real-space, possibly different along each grid-dimension.
% 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
% -------
% varargout : tuple of numerical arrays
% 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
% --------
......@@ -64,19 +60,11 @@ dx = dx(:).' .* ones([1,numel(N)]);
ndim = numel(N);
% Grid in Fourier space
xi = cell(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)) ) );
end
% Optionally assemble ndgrid (occupies more memory!)
if nargin == 3 && computeMeshgrid
[xi{:}] = ndgrid(xi{:});
elseif ndim > 1
for jj = 1:ndim
xi{jj} = reshape( xi{jj}, [ones([1,jj-1]), N(jj), ones([1,ndim-jj])] );
end
end
% Unpack cell-array
varargout = xi;
......
......@@ -138,7 +138,7 @@ function kernel = fresnelPropagationKernel_fourier(N, fresnelNumbers, gpuFlag)
ndim = numel(N);
xi = cell([ndim,1]);
[xi{:}] = fftfreq(N, 1, false);
[xi{:}] = fftfreq(N, 1);
kernel = 1;
for dim = 1:ndim
xi{dim} = gpuArrayIf(xi{dim}, gpuFlag);
......@@ -154,7 +154,7 @@ function kernel = fresnelPropagationKernel_chirp(N, fresnelNumbers, gpuFlag)
ndim = numel(N);
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)));
for dim = 1:ndim
x{dim} = gpuArrayIf(x{dim}, gpuFlag);
......@@ -170,7 +170,7 @@ function kernel = fresnelPropagationKernel_chirpLimited(N, fresnelNumbers, gpuFl
ndim = numel(N);
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)));
for dim = 1:ndim
x{dim} = gpuArrayIf(x{dim}, gpuFlag);
......@@ -196,7 +196,7 @@ function kernel = fresnelPropagationKernel_pixel2pixel(N, fresnelNumbers, gpuFla
ndim = numel(N);
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;
for dim = 1:ndim
x{dim} = gpuArrayIf(x{dim}, gpuFlag);
......
......@@ -28,8 +28,8 @@ function [averages, radii] = angularAverage(im)
%
% N = [1024,1024];
% fresnel = 1e-2;
% [xi_1,xi_2] = fftfreq(N, 1, true);
% ctf = fftshift(sin((xi_1.^2 + xi_2.^2)/(4*pi*fresnel)));
% xi2 = fftfreqNormSq(N, 1);
% ctf = fftshift(sin(xi2/(4*pi*fresnel)));
% [ctfAveraged, radii] = angularAverage(ctf);
%
% figure('name', 'image'); imagesc(ctf); colorbar;
......
......@@ -29,7 +29,7 @@ end
ndim = numel(sizeHolo);
xi = cell([ndim,1]);
[xi{:}] = fftfreq(sizeHolo, 1, false);
[xi{:}] = fftfreq(sizeHolo, 1);
XiSqDivFresnel = 0;
for jj = 1:ndim
XiSqDivFresnel = XiSqDivFresnel + gpuArrayIf(xi{jj}, gpuFlag).^2 ./ (2*pi^2*fresnelNumbers(jj));
......
......@@ -6,8 +6,7 @@ function sinoCorrected = ringremove(sino,settings)
% 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
% 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
% wavelet-based approach are taken from the original paper
% by Muench et al. :cite:`Muench_OE_2009`, denoted as 'wavelet'.
%
% Parameters
% ----------
......@@ -38,7 +37,7 @@ function sinoCorrected = ringremove(sino,settings)
% Highest decomposition level in the wavelet decomposition when using method
% 'wavelet'.
% 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'.
%
% Returns
......@@ -123,7 +122,6 @@ settings = completeStruct(settings, defaults);
switch settings.method
case 'additive'
numAngles = size(sino,2);
% averaging of sinogram along the angle
sinoMean = mean(sino,2);
......@@ -140,50 +138,43 @@ switch settings.method
end
% identification of horizontal stripes
correctionProfile = sinoMean - sinoMeanSmoothed;
% creation of 2d correction mask
correction2D = repmat(correctionProfile,[1 numAngles]);
% sinogram correction
sinoCorrected = sino - settings.additive_Strength*correction2D;
sinoCorrected = sino - settings.additive_Strength * correctionProfile;
case 'wavelet'
numAngles = size(sino,2);
sizeProjs = size(sino,1);
% wavelet decomposition
for decompIdx = 1:settings.wavelet_DecNum
[sino,horzPart{decompIdx},vertPart{decompIdx},diagPart{decompIdx}] = ...
dwt2(sino,settings.wavelet_Wname);
% perform multilevel wavelet decomposition
H = cell(settings.wavelet_DecNum, 1);
V = cell(settings.wavelet_DecNum, 1);
D = cell(settings.wavelet_DecNum, 1);
S = cell(settings.wavelet_DecNum, 1);
for l = 1:settings.wavelet_DecNum
S{l} = size(sino);
[sino, H{l}, V{l}, D{l}] = dwt2(sino,settings.wavelet_Wname);
end
% FFT transform of horizontal frequency bands
for decompIdx = 1:settings.wavelet_DecNum
% filter horizontal frequency bands
for l = 1:settings.wavelet_DecNum
% FFT
horzFourier = fftshift(fft(horzPart{decompIdx}'));
[my,mx] = size(horzFourier);
horzFourier = fft(H{l}');
% damping of horizontal stripe information
dampGauss = 1 - exp(-(-floor(my/2):-floor(my/2) + my - 1).^2/(2*settings.wavelet_Sigma^2));
horzFourier = horzFourier.*repmat(dampGauss',1,mx);
m = size(horzFourier,1);
k = fftfreq(m, 2*pi/m);
dampGauss = 1 - exp(-1/2 * k.^2 / settings.wavelet_Sigma^2 );
horzFourier = horzFourier .* dampGauss;
% inverse FFT
horzPart{decompIdx} = ifft(ifftshift(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);
H{l} = ifft(horzFourier)';
end
% make sure that corrected sinogram has the same size as the original sinogram
if size(sinoCorrected,2) == numAngles + 1
sinoCorrected = sinoCorrected(:,1:end - 1);
end
if size(sinoCorrected,1) == sizeProjs + 1
sinoCorrected = sinoCorrected(1:end - 1,:);
end
% wavelet reconstruction
sinoCorrected = sino;
for l = settings.wavelet_DecNum:-1:1
sinoCorrected = sinoCorrected(1:size(V{l},1),1:size(V{l},2));
sinoCorrected = idwt2(sinoCorrected,H{l}, V{l}, D{l},settings.wavelet_Wname, S{l});
end
otherwise
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