Added function to test the Fresnel number of a hologram

functions/phaseRetrieval/utilities/checkFresnelNumber.m

+function checkFresnelNumber(hologram,fresnelNumberExpected,settings)
+% checkFresnelNumber generates a plot of the power-spectral-density (PSD) of a given hologram
+% together with the expected contrast-transfer-function (CTF). Thereby, it allows to check
+% correctness of the assumed Fresnel number.
+%
+%    ``checkFresnelNumber(hologram,fresnelNumberExpected,settings)``
+% 
+%
+% Parameters
+% ----------
+% hologram : numerical array
+%     hologram for which the Fresnel number is to be checked
+% fresnelNumberExpected : float
+%     pixel size-based Fresnel number at which the hologram was presumably acquired
+% settings : structure, optional
+%     contains additional settings, see *Other Parameters*
+%
+% Other Parameters 
+% ---------------- 
+% betaDeltaRatio : Default = 0
+%     Assumed ratio between absorption and phase shifts of the imaged object.
+% maxNumPlottedCTFPeriods : Default = 5
+%     Maximum number of periods of the contrast-transfer-function included in the plot
+%     of the hologram's PSD
+%
+% Returns 
+% ------- 
+%
+% See also
+% --------
+% functions.phaseRetrieval.holographic.phaserec_ctf
+%
+% Example
+% -------
+%
+% .. code-block:: matlab
+%
+%     data = getHologram('macrophage');
+%     checkFresnelNumber(data.holograms(:,:,1), data.fresnelNumbers(1,1));
+%
+%
+% See also PHASEREC_CTF
+    
+% 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/>.
+
+% Complete settings with defaults
+defaults.betaDeltaRatio = 0;
+defaults.maxNumPlottedCTFPeriods = 5;  
+if (nargin < 3)
+    settings = struct;
+end       
+settings = completeStruct(settings,defaults);
+
+
+% Check input
+if ndims(hologram) ~= 2 || any(size(hologram) == 1)
+	error('Input hologram must be a single 2D-image.');
+end
+
+
+% Ensure that hologram is square in order to allow radially averaging its Fourier transform
+hologram = cropToCenter(hologram, min(size(hologram))*[1,1]);
+
+% Apply fadeout-window to hologram in order to suppress boundary effects when computing fft
+hologram = fadeoutImage(hologram);
+
+% Compute angularly averaged, logarithmic power-spectral-density (PSD)
+hologramPSDAngularAvg = angularAverage(psd(double(hologram)));
+
+% Fourier-frequencies corresponding to the PSD-values, normalized with Fresnel number
+fourierFreqs = fftfreq(size(hologram,1));
+fourierFreqs = fourierFreqs(1:ceil(end/2)) ./ sqrt(4*pi*fresnelNumberExpected);
+
+% Restrict the frequency range to the specified number of CTF-periods
+maxFreq = sqrt(2*pi*settings.maxNumPlottedCTFPeriods);
+hologramPSDAngularAvg = hologramPSDAngularAvg(fourierFreqs < maxFreq);
+fourierFreqs = fourierFreqs(fourierFreqs < maxFreq);
+
+% Compute expected contrast-transfer-function CTF with fine sampling 
+fourierFreqsCTF = linspace(0, min(pi./sqrt(4*pi*fresnelNumberExpected), maxFreq), 10000).';
+ctfExpected = (sin(fourierFreqsCTF.^2) + settings.betaDeltaRatio * cos(fourierFreqsCTF.^2)).^2;
+
+% Scale CTF for better visibility in joint plot with PSD
+minPSD = min(hologramPSDAngularAvg(:));
+maxPSDHighFreq = max(hologramPSDAngularAvg(fourierFreqs > sqrt(pi/2)));
+ctfExpected = (maxPSDHighFreq./max(ctfExpected(:))) * ctfExpected;
+
+% Jointly plot computed PSD and expected CTF with logarithmic y-axis
+figure; plot(fourierFreqs, hologramPSDAngularAvg, '-b', fourierFreqsCTF, ctfExpected, '--r');
+set(gca, 'YScale', 'log');
+ylim([0.8*minPSD, 1.25*maxPSDHighFreq]);
+xlim([0, maxFreq]);
+legend({'PSD(hologram)', 'expected CTF'}, 'location', 'northeast');
+title(['Check that the location of the minima in the two curves coincide. ', ...
+       'If the minima differ systematically, the entered Fresnel number is probably wrong']);
+
+end
+
+