Commit 01cd201d authored by Simon Maretzke's avatar Simon Maretzke
Browse files

rudimentary docu + improved naming of ctfRegWeights/ctfRegMatrix

parent ef6edc27
...@@ -196,9 +196,9 @@ end ...@@ -196,9 +196,9 @@ end
% Assemble CTF-inversion formula: % Assemble CTF-inversion formula:
% %
% volPhase = argmin_f { sum_j ||2*ctf{j}.* FT(f) - FT(holograms{j}-1)||_2^2 % volPhase = argmin_f { sum_j ||2*ctf{j}.* FT(f) - FT(holograms{j}-1)||_2^2
% + 2*||sqrt(regmatrix) .* FT(f)||_2^2 } % + 2*||sqrt(regWeights) .* FT(f)||_2^2 }
% %
% = iFT( (sum_j ctf{j} .* FT(holograms{j}-1)) / (sum_j ctf{j}.^2 + regmatrix) ) % = iFT( (sum_j ctf{j} .* FT(holograms{j}-1)) / (sum_j ctf{j}.^2 + regWeights) )
% %
% ctf{j} = sin(xi^2/(4*pi*fresnel(j))) + settings.betaDeltaRatio * cos(xi^2/(4*pi*fresnel(j))); % ctf{j} = sin(xi^2/(4*pi*fresnel(j))) + settings.betaDeltaRatio * cos(xi^2/(4*pi*fresnel(j)));
% %
...@@ -236,15 +236,15 @@ end ...@@ -236,15 +236,15 @@ end
sumCTFSq = 2*sumCTFSq; sumCTFSq = 2*sumCTFSq;
% Assemble regularization matrix in Fourier-space that smoothly transitions from the value % Assemble regularization weights in Fourier-space that smoothly transitions from the value
% lim1 in the low-frequency regime (around the central CTF-minimum) to lim2 at larger % lim1 in the low-frequency regime (around the central CTF-minimum) to lim2 at larger
% larger spatial frequencies beyond the first CTF-maximum % larger spatial frequencies beyond the first CTF-maximum
meanFresnelNumbers = [mean(fresnelNumbers(2,:))*[1;1]; mean(fresnelNumbers(1,:))]; meanFresnelNumbers = [mean(fresnelNumbers(2,:))*[1;1]; mean(fresnelNumbers(1,:))];
regmatrix = ctfRegMatrix(N, meanFresnelNumbers, settings.lim1, settings.lim2); regWeights = ctfRegWeights(N, meanFresnelNumbers, settings.lim1, settings.lim2);
% Add regularization term to sum of squared CTFs % Add regularization term to sum of squared CTFs
sumCTFSqReg = sumCTFSq + regmatrix; clear 'regmatrix'; clear 'sumCTFSq'; sumCTFSqReg = sumCTFSq + regWeights; clear 'regWeights'; clear 'sumCTFSq';
...@@ -258,7 +258,7 @@ if ~isempty(settings.support) || settings.minPhasePerVoxel > -inf || settings.ma ...@@ -258,7 +258,7 @@ if ~isempty(settings.support) || settings.minPhasePerVoxel > -inf || settings.ma
% Assemble proximal map of the regularized CTF-functional % Assemble proximal map of the regularized CTF-functional
% %
% F(f) := ||2*ctf{j}.*FT(f) - FT(holograms{j}-1)||_2^2 % F(f) := ||2*ctf{j}.*FT(f) - FT(holograms{j}-1)||_2^2
% + 2*||sqrt(regmatrix) .* FT(f)||_2^2 % + 2*||sqrt(regWeights) .* FT(f)||_2^2
% %
proxCTF = @(f, sigma) real(ifftn( (sumCTFHolograms + fftn((1./sigma)*f)) ... proxCTF = @(f, sigma) real(ifftn( (sumCTFHolograms + fftn((1./sigma)*f)) ...
./ (sumCTFSqReg + 1./sigma) )); ./ (sumCTFSqReg + 1./sigma) ));
......
...@@ -56,8 +56,6 @@ function result = phaserec_ctf(holograms, fresnelNumbers, settings) ...@@ -56,8 +56,6 @@ function result = phaserec_ctf(holograms, fresnelNumbers, settings)
% ------- % -------
% result : numerical array % result : numerical array
% The reconstructed phase image of the same size as the input holograms % The reconstructed phase image of the same size as the input holograms
% regmatrix : numerical array
% The regularization weights applied in Fourier-space
% %
% %
% See also % See also
...@@ -187,9 +185,9 @@ N = size(holograms(:,:,1)); ...@@ -187,9 +185,9 @@ N = size(holograms(:,:,1));
% Assemble CTF-inversion formula: % Assemble CTF-inversion formula:
% %
% result = argmin_f { sum_j ||2*ctf{j}.* fft2(f) - fft2(holograms{j}-1)||_2^2 % result = argmin_f { sum_j ||2*ctf{j}.* fft2(f) - fft2(holograms{j}-1)||_2^2
% + 2*||sqrt(regmatrix) .* fft2(f)||_2^2 } % + 2*||sqrt(regWeights) .* fft2(f)||_2^2 }
% %
% = ifft2( (sum_j ctf{j} .* fft2(holograms{j}-1)) / (sum_j ctf{j}.^2 + regmatrix) ) % = ifft2( (sum_j ctf{j} .* fft2(holograms{j}-1)) / (sum_j ctf{j}.^2 + regWeights) )
% %
% ctf{j} = sin(xi^2/(4*pi*F(j))) + settings.betaDeltaRatio * cos(xi^2/(4*pi*F(j))); % ctf{j} = sin(xi^2/(4*pi*F(j))) + settings.betaDeltaRatio * cos(xi^2/(4*pi*F(j)));
% %
...@@ -212,10 +210,10 @@ sumCTFSq = 2*sumCTFSq; ...@@ -212,10 +210,10 @@ sumCTFSq = 2*sumCTFSq;
sumCTFHolograms(1,1) = sumCTFHolograms(1,1) - prod(N) * numHolos * settings.betaDeltaRatio; sumCTFHolograms(1,1) = sumCTFHolograms(1,1) - prod(N) * numHolos * settings.betaDeltaRatio;
% Assemble regularization matrix in Fourier-space that smoothly transitions from the value % Assemble regularization weights in Fourier-space that smoothly transitions from the value
% lim1 in the low-frequency regime (around the central CTF-minimum) to lim2 at larger % lim1 in the low-frequency regime (around the central CTF-minimum) to lim2 at larger
% larger spatial frequencies beyond the first CTF-maximum % larger spatial frequencies beyond the first CTF-maximum
regmatrix = ctfRegMatrix(N, mean(fresnelNumbers,2), settings.lim1, settings.lim2); regWeights = ctfRegWeights(N, mean(fresnelNumbers,2), settings.lim1, settings.lim2);
...@@ -236,9 +234,9 @@ if numel(settings.support) || settings.minPhase > -inf || settings.maxPhase < i ...@@ -236,9 +234,9 @@ if numel(settings.support) || settings.minPhase > -inf || settings.maxPhase < i
% Assemble proximal map of the regularized CTF-functional % Assemble proximal map of the regularized CTF-functional
% %
% F(f) := ||2*ctf{j}.*fft2(f) - fft2(holograms{j}-1)||_2^2 % F(f) := ||2*ctf{j}.*fft2(f) - fft2(holograms{j}-1)||_2^2
% + 2*||sqrt(regmatrix) .* fft2(f)||_2^2 % + 2*||sqrt(regWeights) .* fft2(f)||_2^2
% %
sumCTFSqReg = gpuArrayIf(sumCTFSq + regmatrix, useGPU); sumCTFSqReg = gpuArrayIf(sumCTFSq + regWeights, useGPU);
sumCTFHolograms = gpuArrayIf(sumCTFHolograms, useGPU); sumCTFHolograms = gpuArrayIf(sumCTFHolograms, useGPU);
proxCTF = @(f, sigma) real(ifft2( (sumCTFHolograms + fft2((1./sigma)*f)) ... proxCTF = @(f, sigma) real(ifft2( (sumCTFHolograms + fft2((1./sigma)*f)) ...
./ (sumCTFSqReg + 1./sigma) )); ./ (sumCTFSqReg + 1./sigma) ));
...@@ -276,7 +274,7 @@ if numel(settings.support) || settings.minPhase > -inf || settings.maxPhase < i ...@@ -276,7 +274,7 @@ if numel(settings.support) || settings.minPhase > -inf || settings.maxPhase < i
% ---------- Default case: Minimization without constraints. Can be performed directly ---------- % % ---------- Default case: Minimization without constraints. Can be performed directly ---------- %
% ----------------------------------------------------------------------------------------------- % % ----------------------------------------------------------------------------------------------- %
else else
result = real(ifft2(sumCTFHolograms ./ (sumCTFSq + regmatrix))); result = real(ifft2(sumCTFHolograms ./ (sumCTFSq + regWeights)));
end end
......
...@@ -227,13 +227,13 @@ parOp.betaDeltaRatio = settings.betaDeltaRatio; ...@@ -227,13 +227,13 @@ parOp.betaDeltaRatio = settings.betaDeltaRatio;
% Forward operator F that maps a given phase-image onto the corresponding holograms % Forward operator F that maps a given phase-image onto the corresponding holograms
F = NonlinearPhaseContrastOperator(N, fresnelNumbers, parOp); F = NonlinearPhaseContrastOperator(N, fresnelNumbers, parOp);
% %
% Tikhonov functional(f) = || F(f) - holograms ||_2^2 + 2*|| regmatrix * fft2(f) ||_2^2, % Tikhonov functional(f) = || F(f) - holograms ||_2^2 + 2*|| regWeights * fft2(f) ||_2^2,
% where regmatrix is the same regularization matrix 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)
holograms = gpuArrayIf(holograms, useGPU); holograms = gpuArrayIf(holograms, useGPU);
regmatrixTimes2 = 2*ctfRegMatrix(N, mean(fresnelNumbers,2), settings.lim1, settings.lim2, useGPU); regWeightsTimes2 = 2*ctfRegWeights(N, mean(fresnelNumbers,2), settings.lim1, settings.lim2, useGPU);
TikhonovFunctional = NormSqFunctional * (F - holograms) ... TikhonovFunctional = NormSqFunctional * (F - holograms) ...
+ NormSqFunctional(@(f) real(ifft2(regmatrixTimes2 .* fft2(f)))); + NormSqFunctional(@(f) real(ifft2(regWeightsTimes2 .* fft2(f))));
......
function regmatrix = ctfRegMatrix(sizeHolo, fresnelNumbers, lim1, lim2, gpuFlag)
if nargin < 5
gpuFlag = false;
end
ndim = numel(sizeHolo);
xi = cell([ndim,1]);
[xi{:}] = fftfreq(sizeHolo, 1, false);
XiSqDivFresnel = 0;
for jj = 1:ndim
XiSqDivFresnel = XiSqDivFresnel + gpuArrayIf(xi{jj}, gpuFlag).^2 ./ (2*pi^2*fresnelNumbers(jj));
end
sigma = (sqrt(2) - 1) / 2;
w = 1/2 * erfc((sqrt(XiSqDivFresnel) - 1) / (sqrt(2) * sigma));
regmatrix = lim1 .* w + lim2 .* (1 - w);
end
function regWeights = ctfRegWeights(sizeHolo, fresnelNumbers, lim1, lim2, gpuFlag)
% Computes a regularization weights in Fourier space, suitable for regularized inversion
% 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)
% regularization parameter lim1 at low Fourier frequencies around the central CTF-mimimum
% to a second (typically value) value lim2 that determines the regularization in the
% higher Fourier frequencies beyond the first CTF-maximum.
% See phaserec_ctf, phaserec_nonlinTikhonov and phaserec3D_ctf for typical usage.
% HoloTomoToolbox
% Copyright (C) 2019 Institut fuer Roentgenphysik, Universitaet Goettingen
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
if nargin < 5
gpuFlag = false;
end
ndim = numel(sizeHolo);
xi = cell([ndim,1]);
[xi{:}] = fftfreq(sizeHolo, 1, false);
XiSqDivFresnel = 0;
for jj = 1:ndim
XiSqDivFresnel = XiSqDivFresnel + gpuArrayIf(xi{jj}, gpuFlag).^2 ./ (2*pi^2*fresnelNumbers(jj));
end
sigma = (sqrt(2) - 1) / 2;
w = 1/2 * erfc((sqrt(XiSqDivFresnel) - 1) / (sqrt(2) * sigma));
regWeights = lim1 .* w + lim2 .* (1 - w);
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