...
 
Commits (2)
......@@ -67,6 +67,10 @@ function imTransformed = rotateImage(im, rotAngleDegree, settings)
% 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 == 0
imTransformed = rmfield(shiftRotateMagnifyImage, 'invertTransform');
return
end
if nargin < 3
settings = struct;
end
......
......@@ -66,6 +66,10 @@ function imTransformed = shiftImage(im, shifts, settings)
% 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 == 0
imTransformed = rmfield(shiftRotateMagnifyImage, 'invertTransform');
return
end
if nargin < 3
settings = struct;
end
......
......@@ -87,6 +87,10 @@ function imTransformed = shiftRotateImage(im, shifts, rotAngleDegree, settings)
% 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 == 0
imTransformed = shiftRotateMagnifyImage;
return
end
if nargin < 4
settings = struct;
end
......
......@@ -236,7 +236,8 @@ F = NonlinearPhaseContrastOperator(N, fresnelNumbers, parOp);
% 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
% (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) ...
+ 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
% 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)
......@@ -32,11 +32,27 @@ xi = cell([ndim,1]);
[xi{:}] = fftfreq(sizeHolo, 1);
XiSqDivFresnel = 0;
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
sigma = (sqrt(2) - 1) / 2;
w = 1/2 * erfc((sqrt(XiSqDivFresnel) - 1) / (sqrt(2) * sigma));
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