Commit b3975a7e authored by Simon Maretzke's avatar Simon Maretzke
Browse files

Impose stronger regularization for spatial frequencies beyond the numerical...

Impose stronger regularization for spatial frequencies beyond the numerical aperture to prevent high-frequency artifacts in the case of deeply holographic data
parent 177797b5
Pipeline #133078 passed with stage
in 1 minute and 46 seconds
...@@ -236,7 +236,8 @@ F = NonlinearPhaseContrastOperator(N, fresnelNumbers, parOp); ...@@ -236,7 +236,8 @@ F = NonlinearPhaseContrastOperator(N, fresnelNumbers, parOp);
% Tikhonov functional(f) = || F(f) - holograms ||_2^2 + 2*|| regWeights * fft2(f) ||_2^2, % Tikhonov functional(f) = || F(f) - holograms ||_2^2 + 2*|| regWeights * fft2(f) ||_2^2,
% where regWeights contains the same Fourier-space regularization weights as used in phaserec_ctf % where regWeights contains the same Fourier-space regularization weights as used in phaserec_ctf
% (note that "*" denotes composition of Operators/functionals in the following statements) % (note that "*" denotes composition of Operators/functionals in the following statements)
regWeightsTimes2 = 2*ctfRegWeights(N, mean(fresnelNumbers,2), settings.lim1, settings.lim2, useGPU); regWeightsTimes2 = 2*ctfRegWeights(N, mean(fresnelNumbers,2), settings.lim1, settings.lim2, ...
useGPU, 2*numHolos);
TikhonovFunctional = NormSqFunctional * (F - holograms) ... TikhonovFunctional = NormSqFunctional * (F - holograms) ...
+ NormSqFunctional(@(f) real(ifft2(regWeightsTimes2 .* fft2(f)))); + NormSqFunctional(@(f) real(ifft2(regWeightsTimes2 .* fft2(f))));
......
function regWeights = ctfRegWeights(sizeHolo, fresnelNumbers, lim1, lim2, gpuFlag) function regWeights = ctfRegWeights(sizeHolo, fresnelNumbers, lim1, lim2, gpuFlag, lim3)
% Computes a regularization weights in Fourier space, suitable for regularized inversion % Computes a regularization weights in Fourier space, suitable for regularized inversion
% of the CTF and similar reconstruction methods. The main idea, originally proposed by % of the CTF and similar reconstruction methods. The main idea, originally proposed by
% Peter Cloetens, is to create a mask that smoothly transitions from one (typically low) % Peter Cloetens, is to create a mask that smoothly transitions from one (typically low)
...@@ -32,11 +32,27 @@ xi = cell([ndim,1]); ...@@ -32,11 +32,27 @@ xi = cell([ndim,1]);
[xi{:}] = fftfreq(sizeHolo, 1); [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)); xi{jj} = gpuArrayIf(xi{jj}, gpuFlag);
XiSqDivFresnel = XiSqDivFresnel + xi{jj}.^2 ./ (2*pi^2*fresnelNumbers(jj));
end end
sigma = (sqrt(2) - 1) / 2; sigma = (sqrt(2) - 1) / 2;
w = 1/2 * erfc((sqrt(XiSqDivFresnel) - 1) / (sqrt(2) * sigma)); w = 1/2 * erfc((sqrt(XiSqDivFresnel) - 1) / (sqrt(2) * sigma));
regWeights = lim1 .* w + lim2 .* (1 - w); regWeights = lim1 .* w + lim2 .* (1 - w);
% Additional regularization level for spatial frequencies beyond the
% numerical aperture of the optical system. Only relevant for
% deeply holographic data
if nargin == 6
fresnelNumbersFOV = fresnelNumbers(:) .* sizeHolo(:).^2;
freqCutoffFOV = pi .* (fresnelNumbersFOV./sizeHolo(:));
wCutoffFOV = 1;
XiSqDivCutoffFOV = 0;
for jj = 1:ndim
XiSqDivCutoffFOV = XiSqDivCutoffFOV + (xi{jj}./freqCutoffFOV(jj)).^2;
wCutoffFOV = wCutoffFOV .* (0.5 .* erfc((abs(xi{jj})./freqCutoffFOV(jj)-1)./0.1));
end
regWeights = wCutoffFOV .* regWeights + (1-wCutoffFOV) .* lim3;
end
end end
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