template_tomo_analysis_GINIX.m 16.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
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
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
%% Start from scratch
clear 'variables';  clc;
set(0,'DefaultFigureColor','w');



%% Parameters of the scan
year = '2016';
runname = 'run56_LTP';
date = '20161116';
detector = 'pirra';
filePrefix = '12816-6_Kleinhirn3_tomo7';



%% (modify if loading fails)
if ispc
    toolboxPath = 'S:/Projects_X-ray_Imaging/holotomotoolbox/functions';
    astraPath = 'S:/Projects_X-ray_Imaging/ASTRAToolbox/windows';
    prePath = 'S:/Messzeiten/';

elseif isunix
    toolboxPath = '/IRP/AG_Salditt/Projects_X-ray_Imaging/holotomotoolbox/functions';
    astraPath = '/IRP/AG_Salditt/Projects_X-ray_Imaging/ASTRAToolbox/linux';
    prePath = '/IRP/AG_Salditt/Messzeiten/';
end

dataPath = fullfile(prePath, year, 'GINIX', runname, '/data/current/raw', date, 'detectors', detector, filePrefix);

if ~exist(dataPath,'dir')
    error(['Assumed data path ''', dataPath, ''' does not exist. Scan parameters might be wrong. Please fix the path manually.']);
end

if ~exist(toolboxPath,'dir')
    error(['Assigned path to the holotomotoolbox ''', toolboxPath, ''' does not exist.']);
end

if ~exist(astraPath,'dir')
    error(['Assigned path to the ASTRA-toolbox ''', astraPath, ''' does not exist.']);
end  

%% Add holotomotoolbox to path
addpath(genpath(toolboxPath));

%% Add ASTRA-toolbox to path
addpath(genpath(astraPath));

%% Set up image-reader for the specified detector
imageReader = getImageReaderIRP(detector);

switch(lower(detector))
    case 'argos'
        fileEnding = '_%04i.tif'; dumpEnding = '_%04i.tif.txt';
        dx = 0.54e-6;
    case 'talos'
        fileEnding = '_%04i.tif'; dumpEnding = '_%04i.tif.txt';
        dx = 75e-6;
    case 'pirra'
        fileEnding = '_%i.tif'; dumpEnding = '_%04i.tif.txt';
        dx = 6.5e-6;
    case 'thyia'
        fileEnding = '_%i.tif'; dumpEnding = '_%04i.tif.txt';
        dx = 4.5e-6;
    case 'iris'
        fileEnding = '_%04i.img'; dumpEnding = '_%04i.img.txt';
        dx = 6.5e-6;
    case 'timepix'
        fileEnding = '_%04i.raw'; dumpEnding = '_%04i.raw.txt';
        dx = 55e-6;
    case 'zyla'
        fileEnding = '_%04i.tif'; dumpEnding = '_%04i.tif.txt';
        dx = 6.5e-6;
end

fileName = [filePrefix,fileEnding];
dumpName = [filePrefix,dumpEnding];

getFileName = @(number) fullfile(dataPath, sprintf(fileName,number));
getDumpName = @(number) fullfile(dataPath,'dump',sprintf(dumpName,number));

%% Read one image to check if path is correct
number = 100;
raw0 = imageReader(getFileName(number));

showImage(raw0);
colormap gray


%% setup the current run - edit here
fprintf('... add specifications of setup... \n' ) ;

angularRange = 180;

E = 8000;                   % energy in eV
lambda = 12.398/E*1e-10;    % wavelength in m

% detector distance
z02 = 5047;  % in mm

% calibration motor vs real distance to wg
xmotor = 'stx';
delta_x = 0;  % width of the waveguide (silicon channels)
corr_x = 0;     % width of the vertical 1dWG (crossed multilayer)
corr_y = 0;     % width of the horizontal 1dWG (crossed multilayer)


%% set up file numbers - edit here
fprintf('... set up filenumbers... \n' ) ;

numDistances = 4;       % number of defocus distances
numAngles = 1501;       % number of projections
numFlats = 50;          % number of flat fields
numDarks = 50;          % number of dark fields
offset=84;             % number of images prior to flatfields (wg_align)

% Compute file numbers of projections, flatfields and darkfields - do not edit

numImages=offset+2*numFlats+numAngles; %number of images per distance scan


flatFilenum = zeros([2,numDistances,numFlats]);
projFilenum = zeros([numDistances,numAngles]);
for i=1:numDistances
    flatFilenum(1,i,:) = offset+1+(i-1)*numImages:offset+numFlats+(i-1)*numImages;
    projFilenum(i,:) = flatFilenum(1,i,end)+1:flatFilenum(1,i,end)+numAngles;
    flatFilenum(2,i,:) = projFilenum(i,end)+1:projFilenum(i,end)+numFlats;
end

darknum = flatFilenum(2,end,end)+1:flatFilenum(2,end,end)+numDarks;


%% Check ranges of projection-, flat-field- and dark-field images
fprintf('... please check flat, projection and dark... \n' ) ;

for iii = 1 : numDistances
    tmp = imageReader(getFileName(flatFilenum(1,iii,1)));
    showImage(tmp);
    title('Is this a flatfield? (Y- enter) (N - Ctr+c)')
    colormap gray
    pause()
    
    tmp = imageReader(getFileName(projFilenum(iii,1)));
    showImage(tmp);
    title('Is this a projection? (Y- enter) (N - Ctr+c)')
    colormap gray
    pause()
    
end

tmp = imageReader(getFileName(darknum(1)));
showImage(tmp);
title('Is this a darkfield? (Y- enter) (N - Ctr+c)')
colormap gray
pause()


%% read flat fields - do not edit
fprintf('... read flatfields... \n' ) ;

empty = zeros([size(raw0,1),size(raw0,2),numDistances,2], 'single');

for indDistance = 1:numDistances
    for indEmpty = 1:size(flatFilenum,1)
        emptyTmp = zeros([size(raw0,1),size(raw0,2),numFlats], 'single');
        for indImg = 1:numFlats
            emptyTmp(:,:,indImg) = imageReader(getFileName(flatFilenum(indEmpty,indDistance,indImg)));
            fprintf('.');
        end
        % take median of all emptys
        empty(:,:,indDistance,indEmpty) = median(emptyTmp,3);
        fprintf('\n');
    end
end


%% read dark - do not edit
fprintf('... read darkfields... \n' ) ;

darkTmp = zeros(size(raw0,1),size(raw0,2),numDarks);
for indImg = 1:numDarks
    darkTmp(:,:,indImg) = imageReader(getFileName(darknum(indImg)));
    fprintf('.');
end
% take median of all darkfield images
dark = median(darkTmp,3);
fprintf('\n');


%% calculate geometry parameters - do not edit
fprintf('... calculate geometry parameters... \n' ) ;

% Read defocus distances from dump files
stx = zeros([1,numDistances]);
for distIdx = 1:numDistances
    stx(distIdx) = wwas(getDumpName(projFilenum(distIdx,2)),xmotor);
end

% calculate derived parameters: (suffixes 'x' and 'y' account for
% possible astigmatism)
%
% true defocus distances including corrections
z01x = stx+corr_x+delta_x;
z01y = stx+corr_y+delta_x;
%
% distances between object and detetor
z12x = z02-z01x;
z12y = z02-z01y;
%
% magnifications
Mx = z02./z01x;
My = z02./z01y;
%
% effective distance between object and detetor
z_effx = z12x./Mx;
z_effy = z12y./My;
%
% effective pixel size
dxeffx = dx./Mx;
dxeffy = dx./My;
%
% Fresnel numbers in the effective geometry
Fx = dxeffx.^2./(z_effx*lambda);
Fy = dxeffy.^2./(z_effy*lambda);


% Magnifications and Fresnel numbers for all defocus distances (columns), 
% possibly different along the x- and y-dimensions of the images (rows)
M = [My(:).'; Mx(:).'];     
F = [Fy(:).'; Fx(:).'];


%% Read, rescale and align holograms

% Which distances should be read?
distToUse = 1:numDistances;      

% Settings for alignment of the different defocused images. See
% documentation of rescaleDefocusSeries for details and additional options
alignSettings = rescaleDefocusSeries;
alignSettings.alignMethod = 'dftregistration';
alignSettings.medfilt = true;                   % median filter
alignSettings.sigmaHighpass = 50;               % high pass filter
alignSettings.alignReconstructedImages = true;  % should the holograms or single distance CTF reconstructed images be sued for registration?

% settings for the preliminary single-distance CTF-reconstruction used for aligning the images
ctfsettingsAlign = phaserec_ctf;
ctfsettingsAlign.betaDeltaRatio = 1/10;
ctfsettingsAlign.lim2 = 0.1;
alignSettings.ctfsettings = ctfsettingsAlign;

% no changes below this point
fprintf('... reading holograms... \n' ) ;
holograms = zeros([size(raw0,1), size(raw0,2), numel(distToUse), numAngles], 'single');
thetas = zeros([numAngles,1]);
parfor indProj = 1:numAngles
    
    fprintf('projection: %i \n',indProj)
    
    % read in all distances
    I = zeros([size(raw0,1), size(raw0,2), length(distToUse)]);
    idx = 0;
    for indDistance = distToUse
        idx = idx+1;
        % interpolate linearly between the two empty images
        a = (numAngles-indProj)/(numAngles-1);
        emptyInterpolated = a*empty(:,:,indDistance,1)+(1-a)*empty(:,:,indDistance,2);
        
        % read in data and perform empty beam correction
        tmp = imageReader(getFileName(projFilenum(indDistance,indProj)));
        tmpCorr=(tmp-dark)./(emptyInterpolated-dark);
        
        % remove hot pixels
        I(:,:,idx) = removeOutliers(tmpCorr);
        
    end
    
    % read rotation angle
    if indProj~=1
        thetas(indProj) = wwas(getDumpName(projFilenum(distToUse(1),indProj)),'stzrot');
    else
        thetas(indProj)=0;
    end
    
    % rescale and align the images recorded at different defocus distances
    IAligned = rescaleDefocusSeries(I, F(:,distToUse), M(:,distToUse), alignSettings);
    
    % Store aligned hologram
    holograms(:,:,:,indProj) = IAligned;
    
end

% Compute rescaled Fresnel-numbers (cannot be done in parfor-loop)
fresnelNumbers = F .* (M ./ max(M))).^2;

%% Show all holograms for a random tomographic angle to check results
projIdx = randi(numAngles);
for jj = 1:size(holograms,3)
    showImage(holograms(:,:,jj,projIdx));
    colormap gray
    pause(1)
end

%% Show full tomographic data set for first distance to check results
figure; showStack(squeeze(holograms(:,:,1,:)));

%% OPTIONAL: delete corrupted holograms due to missing wedge from data
missingWedge = [];      % Enter range of corrupted holograms here
holograms(:,:,:,missingWedge) = [];
thetas(missingWedge) = [];


%% OPTIONAL: high-pass filtering of holograms to reduce possible artifacts from empty-beam correction
sigmaHighpass = 100;

parfor holoIdx = 1:size(holograms,3)*size(holograms,4)
    holo = holograms(:,:,holoIdx);
    holo = mean(holo(:)) + (holo - imgaussfilt(holo, sigmaHighpass));
    holograms(:,:,holoIdx) = holo;
end


%% OPTIONAL: removal of vertical and/or horizontal stripes in the holograms
%  See documentation of removeStripes for details
stripeSettings = removeStripes;
stripeSettings.direction = 'vertical';      % direction of profile removal: 'vertical', 'horizontal', 'both'
stripeSettings.method = 'additive';         % method of stripe removal: 'additive' or 'multiplicative'
stripeSettings.rangeTop = 200;              % range to average over at top of the image
stripeSettings.rangeBottom = 200;           % range to average over at bottom of the image
stripeSettings.rangeLeft = 200;             % range to average over at left edge of the image
stripeSettings.rangeRight = 200;            % range to average over at right edge of the image

parfor holoIdx = 1:size(holograms,3)*size(holograms,4)
    holo = holograms(:,:,holoIdx);
    holo = removeStripes(holo, stripeSettings);
    holograms(:,:,holoIdx) = holo;
end




%% OPTIONAL: save preprocessed holograms as raw (+ save setup parameters as mat)
save('-v7', [filePrefix, '_paramaters.mat'], 'thetas', 'fresnelNumbers');
savePrefix = [filePrefix, '_preprocessedHolograms'];
saveraw(holograms, savePrefix);

%% OPTIONAL: Load holograms - insert correct filenames, data type and size
parameterFile = 'parameterFile.mat';
hologramsFile = 'hologramsFile.raw';
hologramsDataType = 'single';
hologramsSize = [2048,2048,4,1501];

load(parameterFile);
clear 'holograms';
holograms = readraw(hologramsFile, hologramsDataType, hologramsSize);




%% OPTIONAL: find correct Fresnelnumber by comparing radially averaged PSF of first hologram to expected CTF: change 'scale' until minima coincide
fprintf('... check Fresnel number... \n' ) ;

projIdx = 1;
distIdx = 1;
scale = 1;
fresnelNumberExpected = scale * fresnelNumbers(1,distIdx);

settings.betaDeltaRatio = 0;
checkFresnelNumber(holograms(:,:,distIdx,projIdx), fresnelNumberExpected, settings);




%% PHASE RECONSTRUCTION
%% First determine suitable method and parameters for a random projection
projIdx = randi(size(holograms,4));
holos = double(holograms(:,:,:,projIdx));
N = [size(holos,1),size(holos,2)];

settingsRecon = struct;
settingsRecon.lim1 = 1e-3;                          % Regularization parameter for low frequencies
settingsRecon.lim2 = 3e-2;                          % Regularization parameter for high frequencies
settingsRecon.betaDeltaRatio  = 1/100;              % Estimated ratio between absorption and phase shifts
settingsRecon.padx = 128;                           % Horizontal and...
settingsRecon.pady = 128;                           % vertical padding: use value > 0 if strong artifacts occur at the image boundaries
%settingsRecon.maxPhase = 0;                        % Uncomment to impose negative phases

% Uncomment to impose support constraint via a binary mask (choice below = 
% square support of half the aspect length of the full field-of-view)
%settingsRecon.support = padarray(true(N/2), N/4);  


% CTF-reconstruction (linear)
phaseRecCTF = phaserec_ctf(holos, fresnelNumbers, settingsRecon);

figure('name', 'CTF-reconstruction'); showImage(phaseRecCTF);       


% Nonlinear reconstruction: if the result differs significantly from
% the CTF-recon (if so, it looks hopefully better...), then it should be 
% considered to reconstruct all holograms with the nonlinear approach
phaseRecNL = phaserec_nonlinTikhonov(holos, fresnelNumbers, settingsRecon);

figure('name', 'Nonlinear reconstruction'); showImage(phaseRecNL); 



%% Reconstruct holograms for all angles with the chosen method and settings
method = 'phaserec_ctf';    % Replace by 'phaserec_nonlinTikhonov' to reconstruct nonlinearly

projs = zeros([size(holograms,1), size(holograms,2), size(holograms,4)], 'single');
for projIdx = 1:size(projs,3)
    disp(projIdx);
    projs(:,:,projIdx) = feval(method, double(holograms(:,:,:,projIdx)), fresnelNumbers, settingsRecon);
end



%% OPTIONAL: save projections as raw
savePrefix = [filePrefix, '_reconstructedProjections'];
saveraw(projs, savePrefix);

%% OPTIONAL: Load projections - insert correct filenames, data type and size
parameterFile = 'parameterFile.mat';
projectionsFile = 'projectionsFile.raw';
projectionsDataType = 'single';
hologramsSize = [2048,2048,1501];

load(parameterFile);
clear 'projs';
projs = readraw(projectionsFile, projectionsDataType, projectionsSize);



%% align rotation axis automatically
proj0Degrees = projs(:,:,1);
[~, idx180Degrees] = min(abs((thetas-thetas(1))-180));
proj180Degrees = projs(:,:,idx180Degrees);

% perform automatic alignment of shift and rotation (might fail for too large cone angles)
settings = alignImages;
settings.registrationMode = 'shiftRotation';
settings.registrationAccuracy = 20;
settings.sigmaLowpass = 0;
settings.medfilt = true;
settings.sigmaHighpass = 50;    % High-pass filtering before alignment to avoid
                                % errors due to low-frequency artifacts
[~,shift,tiltAngle,~] = alignImages(fadeoutImage(proj0Degrees), fadeoutImage(fliplr(proj180Degrees)), settings);
shiftAxis = shift(2)/2;
detectorTilt = tiltAngle/2;

% check result
figure; showImage(shiftRotateImage(proj0Degrees, [0,shiftAxis], detectorTilt) - fliplr(shiftRotateImage(proj180Degrees, [0,shiftAxis], detectorTilt)))
title('Check that images cancel each other out! If not, automatic rotation-axis alignment has possibly failed.')


%% correct for position of rotation axis and detector tilt (only if known!! otherwise find out manually)
parfor indProj = 1:size(projs,3)
    disp(indProj);
    projs(:,:,indProj) = shiftRotateImage(projs(:,:,indProj), [0,shiftAxis], detectorTilt);
end




%% TOMOGRAPHIC RECONSTRUCTION (WITH OPTIONAL RING-REMOVAL)
%% First determine suitable reconstruction parameters by trial-reconstructions of a single tomographic slice
slicenumber = ceil(size(projs,1)/2);

sinoCrop = 50;      % Number of pixels cropped from the edges of the sinogram before reconstruction
sinoPad = 500;      % Number of pixels padded to the edges of the sinogram before reconstruction
sliceSize = 1800;   % Size of the reconstructed slices
settingsTomoRec.outputSize = sliceSize;

% Ring-removal settings: see documentation of ringremove for details
doRingRemoval = false;
settingsRingremove = ringremove;

% Perform reconstruction with optional ring-removal
sino = squeeze(projs(slicenumber,:,:));
if doRingRemoval
    sino = ringremove(sino, settingsRingremove);
end
slice = astraFBP(padarray(sino(sinoCrop+1:end-sinoCrop,:), [sinoPad,0], 'replicate'), thetas, settingsTomoRec);

% View results to see whether parameters have to be adjusted
figure; showImage(slice); colormap bone;


%% Reconstruct all slices with determined parameters

% optional ring removal for all slices
if doRingRemoval
    parfor sliceIdx = 1:size(projs,1)
        projs(sliceIdx,:,:) = ringremove(squeeze(projs(sliceIdx,:,:)), settingsRingremove);
    end
end

slices = astraFBP(padarray(projs(:,sinoCrop+1:end-sinoCrop,:), [0, sinoPad], 'replicate'), thetas, settingsTomoRec);



%% show the reconstructed slices
fprintf('... show reconstructed slices... \n' ) ;
showStack(slices);


%% save results as raw
savePrefix = [filePrefix, '_reconstructedSlices'];
saveraw(slices, savePrefix);



%% OPTIONAL: filtering with Gaussian for better SNR
fprintf('... filter slices with a Gaussian function... \n' ) ;
filterSigma = 1; % standard deviation in pixels
slices = imgaussfilt3(slices,filterSigma);


%% save filtered results as raw
savePrefix = [filePrefix, '_reconstructedSlices_filtered_sigma=',num2str(filterSigma)];
saveraw(slices, savePrefix);