phaserec_nonlinTikhonov.m 12.4 KB
Newer Older
1
function result = phaserec_nonlinTikhonov(holograms, fresnelNumbers, settings)
2
3
4
% PHASEREC_NONLINTIKHONOV reconstructs the phase based on a fully nonlinear
% phase-contrast model using Tikhonov regularization.
%
5
%    ``result = phaserec_nonlinTikhonov(holograms, fresnelNumbers, settings)``
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
%
% The method can be seen as a nonlinear generalization phaserec_ctf in the
% sense that both functions will yield exactly the same result in the limit
% of weak phase shifts (linear contrast regime). However, the present method
% may also yield accurate phase-reconstructions outside this regime, thus
% outperforming phaserec_ctf for some data sets.
% The minization of the Tikhonov-functional is performed by a gradient-descent
% method with adaptive (Barzilai-Borwein) stepsizes, using the result of
% phaserec_ctf as an initial guess.
%
%
% Parameters
% ----------
% holograms : cell-array or numerical array
%     hologram(s) to reconstruct from
21
22
23
24
25
% fresnelNumbers : numerical array
%     pixel size-based Fresnel numbers, at which the holograms were acquired: the columns
%     correspond to possibly multiple propagation distances (size(fresnelNumbers,2) must
%     match the number of holograms). If the array has two rows, the rows are interpreted
%     as possibly distinct Fresnel numbers along the two dimensions of the holograms (astigmatism).
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
% settings : structure, optional
%     contains additional settings, see *Other Parameters*
%
% Other Parameters
% ----------------
% lim1 : Default = 1e-3
%     Regularisation Parameter for low frequencies
% lim2 : Default = 1e-1
%     Regularisation Parameter for high frequencies
% betaDeltaRatio : Default = 0
%     Fixed ratio between absorption and phase shifts to be assumed in the
%     reconstruction. The value 0 corresponds to assuming a pure phase object.
% padx : Default = 0
%     Amount by which the holograms are symmetrically (replicate)-padded along the
%     x-dimension before reconstruction
% pady : Default = 0
%     Amount by which the holograms are symmetrically (replicate)-padded along the 
%     y-dimension before reconstruction
% support : Default = []
%     Array of true- and false-values of the same size as the input holograms, specifying
%     an optional support constraint as a binary mask. If support == [], no support 
%     constraint is imposed.
% minPhase : Default = -inf 
%     Lower bound / minimum-constraint for the reconstructed phases. Imposes that all values
%     in result are >= minPhase. If minPhase == -inf, no minimum-constraint is imposed.
% maxPhase : Default = inf
%     Upper bound / maximum-constraint for the reconstructed phases. Imposes that all values
%     in result are <= maxPhase. If maxPhase == inf, no maximum-constraint is imposed.
% useGPUIfPossible : Default = true
%     Run the iterative optimization on GPU if GPU-computations are supported.
%
% Returns
% -------
% result : numerical array
%     The reconstructed phase image of the same size as the input holograms
% stats : struct
%     Structure that contains statistics on the convergence of the gradient-descent
%     iterations performed to reconstruct the phase.
%
%
% See also
% --------
mtoeppe's avatar
mtoeppe committed
68
% functions.phaseRetrieval.holographic.phaserec_ctf
69
70
71
72
73
74
%
% Example
% -------
%
% .. code-block:: matlab
%
75
%     fresnelNumbers = [1.1,0.9,0.8,1.2]*1e-3;
76
%     phaseTrue = -getPhantom('dicty');
77
78
79
%    
%     simSettings.betaDeltaRatio = 0.01;
%     simSettings.simulateLinear = false;
80
%     holograms = simulateHologram(phaseTrue, fresnelNumbers, simSettings);
81
82
83
84
85
%    
%     settings = simSettings;
%     settings.lim1 = 1e-4;
%     settings.lim2 = 1e-2;
%     
86
87
%     phaseRecCTF = phaserec_ctf(holograms, fresnelNumbers, settings);
%     phaseRecNonlin = phaserec_nonlinTikhonov(holograms, fresnelNumbers, settings);
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
%     
%     subplot(1,3,1)
%     showImage(phaseTrue)
%     title('True phase-image')
%     
%     subplot(1,3,2)
%     showImage(phaseRecCTF)
%     title('Reconstructed phase-image (CTF = linear)');
%     
%     subplot(1,3,3)
%     showImage(phaseRecNonlin)
%     title('Reconstructed phase-image (Nonlinear Tikhonov = nonlinear)');
%     
% 
% 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/>.


121
% Default parameters (widely identical to phaserec_ctf)
122
defaults = phaserec_ctf;
123
124
125
126
127
128
129
%
% Control parameters for the iterative optimization by projected gradient-descent:
% Only modify if you know what you are doing!
defaults.optimization = struct;
defaults.optimization.maxIterations = 50;   % maximum number of iterations
defaults.optimization.tolerance = 1e-2;     % relative tolerance for determining convergence
defaults.optimization.verbose = false;      % display status messages during optimization?
130
131
132

% Return default parameters if called without argument
if (nargin == 0)
133
134
    result = defaults;
    return
135
136
137
138
end

% Complete settings with default parameters
if nargin < 3
139
    settings = struct;
140
141
142
end
settings = completeStruct(settings, defaults);

143
% If holograms are assigned as cell-arrays, convert to numerical arrays
144
if iscell(holograms)
145
    holograms = cat(3, holograms{:});
146
147
148
149
end

% Check sanity of input
numHolos = size(holograms,3);
150
if ~isnumeric(fresnelNumbers)
151
    error('fresnelNumbers must be a single number, vector or matrix. Cell-arrays are not supported.');
152
end
153
if size(fresnelNumbers,2) ~= numHolos
154
    error('size(fresnelNumbers,2) must be equal to the number of assigned holograms.');
155
156
end
if size(fresnelNumbers,1) ~= 1 && size(fresnelNumbers,1) ~= 2
157
158
    error(['size(fresnelNumbers,1) must be either 1 (no astigmatism) or 2', ...
           ' (possibly astigmatic setting).']);
159
160
161
end

% Same Fresnel number along both dimensions of the holograms if fresnelNumbers has only one row
162
fresnelNumbers = fresnelNumbers .* ones([2,numHolos]);
163
164
165
166
167
168



% SPECIAL CASE: if settings.betaDeltaRatio is sufficiently high,
% regularization of low frequencies is omitted
if settings.lim1 < 2*numHolos*settings.betaDeltaRatio^2
169
    settings.lim1 = 0;
170
171
172
173
174
175
end



% Run iterative reconstruction on GPU if desired and possible
useGPU = settings.useGPUIfPossible && checkGPUSupport();
176
if settings.useGPUIfPossible && ~useGPU
177
    warning(['Unable to initialize GPU. Performing optimization on CPU instead. ', ...
178
179
180
181
182
             'Set settings.useGPUIfPossible = false to suppress this warning.']);
end 



Simon Maretzke's avatar
Simon Maretzke committed
183
% Optional padding of holograms (and of optional support)
184
holograms = padarray(holograms, [settings.pady settings.padx], 'replicate');
Simon Maretzke's avatar
Simon Maretzke committed
185
if ~isempty(settings.support)
186
    settings.support = padarray(settings.support, [settings.pady settings.padx]);
Simon Maretzke's avatar
Simon Maretzke committed
187
end
188
189
190
191
192
193
194
195
196
197
N = size(holograms(:,:,1));




% ----------------------------------------------------------------------------------------------- %
% --------- Step 1: perform CTF-reconstruction to obtain initial guess (if not assigned) -------- %
% ----------------------------------------------------------------------------------------------- %

settingsCTF = settings;
198
199
settingsCTF.padx = 0;            % Holograms are already padded,
settingsCTF.pady = 0;            % so no more padding necessary
200
201
settingsCTF.useGPUIfPossible = useGPU;
settingsCTF.optimization = struct;
202
initialGuess = gpuArrayIf(phaserec_ctf(holograms, fresnelNumbers, settingsCTF), useGPU);
203

204
205
206
207
208
209
210
211
212
213
214
% Preliminary nonlinearity-correction in the reconstructed absorption, applied in the low-frequency
% regime. 
% The idea is the following: the CTF-reconstruction fits the mean value of the phase, phi_0, such
% that 1-(2*betaDeltaRatio)*phi_0 ~ I_0 where I_0 is the mean of the measured intensities. In the
% nonlinear model, however, it must hold that exp(-(2*betaDeltaRatio)*phi_0) ~ I_0. The correction
% below accounts for this difference between CTF- and nonlinear model by adjusting initialGuess
% accordingly (not only its zero Fourier-frequency (= mean value) but the whole low-frequency
% part of the image).
% For strongly absorbing objects, this correction greatly improves convergence of the iterative
% scheme applied below.
if settings.betaDeltaRatio > 0
215
216
217
218
    absorptionImage = (2*settings.betaDeltaRatio) * initialGuess;
    absorptionLowFreq = imgaussfilt(absorptionImage, mean(1./sqrt(2*pi*fresnelNumbers(:))));
    initialGuess = initialGuess +    (log(1 + absorptionLowFreq) - absorptionLowFreq) ...
                                  ./ (2*settings.betaDeltaRatio);
219
end
220
221
222
223
224
225
226
227
228
229
230
231



% ----------------------------------------------------------------------------------------------- %
% ------------------ Step 2: Iteratively improve result using nonlinear model ------------------- %
% ----------------------------------------------------------------------------------------------- %

% Assemble functional to minimize
parOp.gpuFlag = useGPU;
parOp.method = 'fourier';
parOp.betaDeltaRatio = settings.betaDeltaRatio;
%
232
% Forward operator F that maps a given phase-image onto the corresponding holograms
233
F = NonlinearPhaseContrastOperator(N, fresnelNumbers, parOp);
234
%
235
236
% 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
237
238
% (note that "*" denotes composition of Operators/functionals in the following statements)
holograms = gpuArrayIf(holograms, useGPU);
239
regWeightsTimes2 = 2*ctfRegWeights(N, mean(fresnelNumbers,2), settings.lim1, settings.lim2, useGPU);
240
TikhonovFunctional =   NormSqFunctional * (F - holograms) ...
241
                     + NormSqFunctional(@(f) real(ifft2(regWeightsTimes2 .* fft2(f))));
242
243
244
245
246



% Assemble projection onto optional support- and min/max-constraints
if settings.minPhase > -inf || settings.maxPhase < inf || ~isempty(settings.support)
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
    Id = @(f) f;
    if settings.minPhase > -inf
        projMin = @(f) max(f, settings.minPhase);
    else
        projMin = Id;
    end
    if settings.maxPhase < inf
        projMax = @(f) min(f, settings.maxPhase);
    else
        projMax = Id;
    end
    if numel(settings.support)
        % If padding is applied to the data, the support has to be padded as well
        support = gpuArrayIf(settings.support, useGPU);
        projSupport = @(f) support .* f;
    else
        projSupport = Id;
    end

    % Proximal operator of the total constraint functional G
    projection = @(f) projMax(projMin(projSupport(f)));
268
269

else
270
    projection = [];
271
272
273
274
end


% Perform reconstruction using projected gradient descent method with Barzilai-Borwein stepsizes
275
276
277
optimization = settings.optimization;
optimization.projection = projection;
optimization.initialStepsize = 1/(4*numHolos*(1+settings.betaDeltaRatio.^2) + max(settings.lim1, settings.lim2));
278
%
279
280
281
282
283
284
285
286
% Reference value for the norm of the gradient of TikhonovFunctional to measure convergence of
% the gradient descent setps. If betaDeltaRatio == 0, gradRef is exactly the gradient of
% TikhonovFunctional evaluated for a constant zero phase. For betaDeltaRatio > 0, the gradient
% is computed for zero-mean holograms to avoid overly large reference-gradients due to constant 
% background contributions (for betaDeltaRatio = 0, this is not necessary as the constant part of 
% holograms does not enter in the gradient)
F.evaluate(zeros(size(initialGuess), class(initialGuess)));
if settings.betaDeltaRatio == 0
287
    gradRef = F.adjoint(holograms);
288
else
289
    gradRef = F.adjoint(holograms - mean(mean(holograms,1),2));
290
end
291
optimization.resiGradRef = gather(norm(gradRef(:)));    
292
%
293
[result, stats] = projectedGradientDescent(TikhonovFunctional, initialGuess, optimization);
294
295
296
297


% Print warning if convergence has not been achiveded to the prescribed tolerance
if ~stats.isConverged
298
299
300
301
    warning(sprintf(['Optimization method has not converged to the required accuracy (tolerance = %.2e) ', ...
                     'within maxIterations = %i steps (final residual was %.2e > tolerance). ', ...
                     'Result may be inaccurate.'], ...
                    optimization.tolerance, optimization.maxIterations, stats.resiGradRel(end)));
302
303
304
305
end


% Undo optional padding
306
result = gather(croparray(result, [settings.pady, settings.padx]));
307
308
309
310


end