rescaleDefocusSeries.m 12 KB
Newer Older
1
function [imAligned,fresnelNumbersRescaled,imFiltered] = rescaleDefocusSeries(images,fresnelNumbers,magnifications,settings)
mtoeppe's avatar
mtoeppe committed
2
3
% RESCALEDEFOCUSSERIES transforms set of images in cone geometry to parallel-beam
% geometry.
4
%
5
%    ``[imAligned,fresnelNumbersRescaled,imFiltered] = rescaleDefocusSeries(images,fresnelNumbers,magnifications,settings)``
mtoeppe's avatar
mtoeppe committed
6
7
8
9
10
11
12
13
14
15
% 
% To account for a varying magnification in a defocus scan obtained in divergent
% geometry, all images are scaled to the effective pixel size of the projection with
% the highest magnification. Subsequently, the images are either shifted by a
% specified amount, aligned to each other in Fourier space or by manual
% identification of a feature of interest in all images to identify overlapping
% regions. By subsequently cutting all projections accordingly, projection images
% with the same field of view and effective pixel size but varying Fresnel numbers
% can be obtained, well suited for the application of multi-distance phase-retrieval
% approaches.
16
%
mtoeppe's avatar
mtoeppe committed
17
18
% Parameters
% ----------
mtoeppe's avatar
mtoeppe committed
19
% images : numerical array or cell-array
mtoeppe's avatar
mtoeppe committed
20
%     set of projections to be aligned
21
22
23
24
25
26
27
28
29
30
% fresnelNumbers : numerical array
%     pixel size-based Fresnel numbers, at which the holograms were acquired: the columns
%     correspond to the different defocus distances (hence, 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).
% magnifications : numerical
%     geometrical magnifications, at which the holograms were acquired: the columns
%     correspond to the different defocus distances (hence, size(magnifications,2) must
%     match the number of holograms). If the array has two rows, the rows are interpreted
%     as possibly distinct magnifications along the two dimensions of the holograms (astigmatism).
mtoeppe's avatar
mtoeppe committed
31
32
33
34
35
36
% settings : structure, optional
%     contains additional settings, see *Other Parameters*
%
% Other Parameters 
% ---------------- 
% alignMethod : Default = 'dftreg'
37
38
%     Method for aligning the images. 'dftreg' uses dftregistration, 'manualFeature' uses
%     a manually selected feature of interest.
39
40
41
42
43
44
% registrationMode : Default = 'shift'
%     Determines whether alignment is performed w.r.t. shifts (case: 'shift'), rotation
%     (case: 'rotation') or both (case: 'shiftRotation')
% registrationAccuracy : Default = 10
%     Accuracy-factor for the image-registration. Shifts and/or rotations between the images
%     are detected to an accuracy of 1/registrationAccuracy lengths of a pixel.
45
% alignReconstructedImages : Default = false
mtoeppe's avatar
mtoeppe committed
46
%     Should alignment be performed for single-distance CTF reconstructed
47
%     projections? 
mtoeppe's avatar
mtoeppe committed
48
% ctfsettings : Default = []
49
%     Settings for the CTF reconstruction, if alignReconstructedImages is true. If it is
mtoeppe's avatar
mtoeppe committed
50
51
52
53
54
55
56
57
%     empty, the settings betaDeltaRatio = 1/100, lim1 = 0.01 and lim2 = 0.1 are
%     chosen.
% sigmaHighpass : Default = 0
%     Standard deviation of the Gaussian function used for highpass filtering. If
%     sigmaHighpass is 0, no filtering will be performed.
% sigmaLowpass : Default = 0
%     Standard deviation of the Gaussian function used for lowpass filtering. If
%     sigmaLowpass is 0, no filtering will be performed.
58
59
% medfilt : Default = false
%     Should median filtering be applied prior to the alignment? 
mtoeppe's avatar
mtoeppe committed
60
61
62
63
64
65
66
67
% cutLeft : Default = 0
%     Amount of cutting at the left edge in pixels.
% cutRight : Default = 0
%     Amount of cutting at the right edge in pixels.
% cutTop : Default = 0
%     Amount of cutting at the top edge in pixels.
% cutBottom : Default = 0
%     Amount of cutting at the bottom edge in pixels.
68
69
% givestatus : Default = false
%     Should status of the alignment progress be printed in Matlab terminal? 
mtoeppe's avatar
mtoeppe committed
70
71
72
%
% Returns 
% ------- 
mtoeppe's avatar
mtoeppe committed
73
74
% imAligned : numerical array
%     numerical array containing the scaled, aligned and cropped projections
75
% fresnelNumbersRescaled : numerical array
mtoeppe's avatar
mtoeppe committed
76
%     array of rescaled Fresnel numbers
mtoeppe's avatar
mtoeppe committed
77
78
% imFiltered : numerical array 
%     exemplary image for which all image processing steps prior to the alignment were applied
mtoeppe's avatar
mtoeppe committed
79
%
mtoeppe's avatar
mtoeppe committed
80
81
% See also
% --------
mtoeppe's avatar
mtoeppe committed
82
% functions.imageProcessing.alignment.alignImages
mtoeppe's avatar
mtoeppe committed
83
%
mtoeppe's avatar
mtoeppe committed
84
85
% Example
% -------
mtoeppe's avatar
mtoeppe committed
86
87
88
%
% .. code-block:: matlab
%
mtoeppe's avatar
mtoeppe committed
89
%     projs = zeros(512,512,4);
90
91
92
%     magnifications = [1,2,3,4];
%     fresnelNumbers = 0.01./magnifications; 
%     shifts = rand([2,4])*10;
Simon Maretzke's avatar
Simon Maretzke committed
93
%     for ii = 1:4
94
%         projs(:,:,ii) = shiftImage(magnifyImage(phantom(512), magnifications(:,ii)), shifts(:,ii));
mtoeppe's avatar
mtoeppe committed
95
96
97
98
99
%     end
% 
%     settings = rescaleDefocusSeries;
%     settings.alignMethod = 'dftreg';
%     settings.sigmaHighpass = 50;
100
%     settings.medfilt = true;
mtoeppe's avatar
mtoeppe committed
101
% 
102
%     [projsAligned, fresnelNumbersRescaled] = rescaleDefocusSeries(projs,fresnelNumbers,magnifications,settings);
mtoeppe's avatar
mtoeppe committed
103
104
%
%     figure(1)
105
106
107
108
%     for jj = 1:4
%         subplot(2,2,jj); showImage(projs(:,:,jj)); title(['Original image ', num2str(jj)]);
%     end
%
mtoeppe's avatar
mtoeppe committed
109
%     figure(2)
110
111
112
113
%     for jj = 1:4
%         subplot(2,2,jj); showImage(projsAligned(:,:,jj)); title(['Rescaled image ', num2str(jj)]);
%     end
%
mtoeppe's avatar
mtoeppe committed
114
115
%
% See also ALIGNIMAGES
mtoeppe's avatar
mtoeppe committed
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
    
% 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/>.
132
133
134
135
136
137

if (nargin < 4)
    settings = struct;
end

% default settings for function alignImages
mtoeppe's avatar
mtoeppe committed
138
139
defaults.sigmaHighpass = 0;        
defaults.sigmaLowpass = 0;         
140
defaults.medfilt = false;              
mtoeppe's avatar
mtoeppe committed
141
142
143
144
145
defaults.cutLeft = 0;             
defaults.cutRight = 0;            
defaults.cutTop = 0;              
defaults.cutBottom = 0;           
defaults.alignMethod = 'dftreg';   
146
147
defaults.registrationMode = 'shift';            
defaults.registrationAccuracy = 10;            
148

mtoeppe's avatar
mtoeppe committed
149
defaults.givestatus = 0;
150
defaults.alignReconstructedImages = false;
151
152
153
154
defaults.ctfsettings = [];

% return default settings
if (nargin == 0)
mtoeppe's avatar
mtoeppe committed
155
    imAligned = defaults;
156
157
158
    return
end

mtoeppe's avatar
mtoeppe committed
159
settings = completeStruct(settings,defaults);
160
161
162

% default CTF settings
if isempty(settings.ctfsettings)
163
164
165
    settings.ctfsettings = phaserec_ctf;
    settings.ctfsettings.betaDeltaRatio = 1/100;
    settings.ctfsettings.lim1 = 0.01;
166
167
168
169
170
    settings.ctfsettings.lim2 = 0.1;
end


%% Verification of correct usage
mtoeppe's avatar
mtoeppe committed
171
172
if  iscell(images)
    images = cat(3, images{:});
173
end
174
175
176
177
numHolos = size(images,3);
%
% If only one hologram is assigned, then output = input
if numHolos<2
178
    imAligned = images;
179
180
    fresnelNumbersRescaled = fresnelNumbers;
    imFiltered = [];
181
182
    return
end
183
184
185
%
% Check that Fresnel numbers are assigned in the right format
if ~isnumeric(fresnelNumbers)
186
    error('fresnelNumbers must be a single number, vector or matrix. Cell-arrays are not supported.');
187
end
p.jhagema's avatar
p.jhagema committed
188
if length(fresnelNumbers) ~= numHolos
189
    error('size(fresnelNumbers,2) must be equal to the number of assigned holograms.');
190
end
191
if size(fresnelNumbers,1) ~= 1 && size(fresnelNumbers,1) ~= 2
192
193
    error(['size(fresnelNumbers,1) must be either 1 (no astigmatism) or 2', ...
           ' (possibly astigmatic setting).']);
194
end
195
196
%
% Same Fresnel number along both dimensions of the holograms if fresnelNumbers has only one row
197
fresnelNumbers = fresnelNumbers .* ones([2,numHolos]);
198
199
200
%
% Check that magnifications are assigned in the right format
if ~isnumeric(magnifications)
201
    error('magnifications must be a single number, vector or matrix. Cell-arrays are not supported.');
202
end
203
if size(magnifications,2) ~= numHolos
204
    error('size(magnifications,2) must be equal to the number of assigned holograms.');
205
206
end
if size(magnifications,1) ~= 1 && size(magnifications,1) ~= 2
207
208
    error(['size(magnifications,1) must be either 1 (no astigmatism) or 2', ...
           ' (possibly astigmatic setting).']);
209
210
211
end
%
% Same magnifications along both dimensions of the holograms if magnifications has only one row
212
magnifications = magnifications .* ones([2,numHolos]);
213
214
215
216
217



%% Perform the rescaling and alignment
% initialize
218
219
220
221
222
imAligned = zeros([size(images,1),size(images,2),numHolos], class(images));

% determine hologram with maximum magnification
[~,idxMaxMagnification] = max(magnifications(1,:));
maxMagnification = magnifications(:,idxMaxMagnification);
223
224

% direction of defocus series
Simon Maretzke's avatar
Bug fix    
Simon Maretzke committed
225
if idxMaxMagnification == numHolos
226
    direction = -1; % from smallest to largest magnification
Simon Maretzke's avatar
Bug fix    
Simon Maretzke committed
227
elseif idxMaxMagnification == 1
228
    direction = 1;  % from largest to smallest magnification
Simon Maretzke's avatar
Bug fix    
Simon Maretzke committed
229
else
230
    error('The magnifications must be either monotonically increasing or decreasing.');
231
232
233
end

% rescale Fresnel numbers
234
fresnelNumbersRescaled = fresnelNumbers .* (magnifications ./ maxMagnification).^2;
235
236

% determine scaling factors with respect to image with largest magnification
237
scalings = maxMagnification ./ magnifications;
238
239

if (settings.givestatus)
240
    disp (['Aligning image 1 / ' num2str(numHolos)])
241
242
243
end

% identify feature for manually determining shifts between images
244
featureCoords = zeros([numHolos,2]);
245
if strcmp(settings.alignMethod,'manualFeature')
246
    showImage(images(:,:,idxMaxMagnification))
247
248
    colormap gray
    % get coordinates (y, x)
249
    featureCoords(idxMaxMagnification,:) = fliplr(ginput(1));
250
251
end

252
% image with largest magnification is reference --> no alignment needed
253
imAligned(:,:,idxMaxMagnification) = images(:,:,idxMaxMagnification);
254
255

% all other images are aligned to the predecessor
mtoeppe's avatar
mtoeppe committed
256
shift = 0; number = 1;
257
for imageIdx = idxMaxMagnification+direction:direction:idxMaxMagnification+direction*(numHolos-1)
258
259
    
    if (settings.givestatus)
mtoeppe's avatar
mtoeppe committed
260
        number = number+1;
261
        disp (['Aligning image ' num2str(number) ' / ' num2str(numHolos)])
262
263
264
265
    end
    
    % scale all images to the pixelsize of the one with the largest
    % magnification
266
    imScaled = magnifyImage(images(:,:,imageIdx), scalings(:,imageIdx));  
267
    
268
    % determine reference image
mtoeppe's avatar
mtoeppe committed
269
    imReference = imAligned(:,:,imageIdx-direction);
270
271
272
    
    % optional CTF reconstruction to better align images to each other as
    % holograms might look too different
273
    if settings.alignReconstructedImages
274
275
        imReference = phaserec_ctf(imReference, fresnelNumbersRescaled(:,imageIdx-direction), settings.ctfsettings);
        imScaled = phaserec_ctf(imScaled, fresnelNumbersRescaled(:,imageIdx), settings.ctfsettings);
276
277
278
    end
           
    % identify feature for manually determining shifts between images
279
    if strcmp(settings.alignMethod,'manualFeature')
280
281
282
283
        imagesc(imScaled)
        colormap gray
        axis equal tight off
        % get coordinates (y, x)
mtoeppe's avatar
mtoeppe committed
284
        featureCoords(imageIdx,:) = fliplr(ginput(1)) - ...
285
286
            [(size(imScaled,1)-size(images(:,:,idxMaxMagnification),1))/2 ...
            (size(imScaled,2)-size(images(:,:,idxMaxMagnification),2))/2];
mtoeppe's avatar
mtoeppe committed
287
        shift_diff = featureCoords(imageIdx-direction,:) - featureCoords(imageIdx,:);
288
289
        shift = shift+shift_diff;
        
290
        imFiltered = imScaled;
291
    else
292
        % align the cropped image to the reference using alignImages and dftregistration
293
294
295
296
297
298
299
300
        keyboard
        [~, shifts, rotAngleDegree, imFiltered] = alignImages(imScaled(...
            settings.cutTop : size(imScaled,1)-settings.cutBottom ...
            , settings.cutLeft : size(imScaled,2)-settings.cutRight) ...
            , imReference(...
            settings.cutTop : size(imReference,1)-settings.cutBottom ...
            , settings.cutLeft : size(imReference,2)-settings.cutRight) ...
            ,settings);
301
302
    end
    
303
    % Apply scaling and registered shifts/rotations to original image (performed in a single step
304
305
306
    % to minimize the degradation of image-quality due to interpolation)
    imAligned(:,:,imageIdx) = shiftRotateMagnifyImage(images(:,:,imageIdx), ...
                                    shifts, rotAngleDegree, scalings(:,imageIdx));
307
308
309
310
311
end

end