-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfixed_mask_calc.m
191 lines (159 loc) · 7.54 KB
/
fixed_mask_calc.m
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% James Manton, 2020 %
% jmanton@mrc-lmb.cam.ac.uk %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SHOW_IMAGES = 1;
SAVE_IMAGES = 1;
ROOT_NAME = 'A12802-35-040_900nm';
wavelength = 900E-9; % Excitation wavelength (900 nm)
max_numerical_aperture = 0.70; % Numerical aperture of excitation lens
refractive_index = 1.33; % Refractive index of focal volume
multiphoton_order = 2; % 1 for single photon, 2 for two-photon, etc.
%% Calculate pupil zone NAs based on physical mask dimensions
% Hamamatsu A12802-35-040 mask
mask_dias = [1.8, 2.5, 3.1, 3.6, 4.0]; % Diameter of mask zones in mm
num_steps = length(mask_dias);
mask_to_pupil_mag = 0.6 * 4; % Magnification from mask to galvos (0.6) to pupil (4)
nas = [0, mask_to_pupil_mag * mask_dias ./ (2 * (180 / 25))]; % Uses focal length of 25x Olympus lens
dof = wavelength / (refractive_index * (1 - cosd(asind(nas(2) / refractive_index))));
%% Set up simulation parameters
M = 150; % PSF planes will have dimensions (2M + 1) by (2M + 1)
N = 1315; % Scaling factor -- increase for smaller voxel size
num_slices = 2 * M + 1; % Set so that PSF stack has same size in all dimensions
z_max = 43.9978E-06; % Set so that Z range and hence voxel size is same as in XY
% Sampling params
m = -M:M;
n = m;
[m_2D, n_2D] = meshgrid(m, n);
k0 = 2 * pi / wavelength;
delta_k = k0 * max_numerical_aperture / M;
R = sqrt(m_2D.^2 + n_2D.^2);
% Aperture stop
stop = zeros(2 * M + 1, 2 * M + 1);
stop(R < M) = 1;
% Coordinates
theta = asin(delta_k .* R / (k0 * refractive_index));
phi = atan2(n_2D, m_2D);
k_xy = delta_k * R;
k_z = sqrt((k0 * refractive_index)^2 - k_xy.^2);
% Distances
z_max = z_max / 2;
z_min = -z_max;
delta_x = M * wavelength / (max_numerical_aperture * N);
x_val = -delta_x * M : delta_x : delta_x * M;
y_val = x_val;
z_val = z_min : (z_max - z_min) / (num_slices - 1) : z_max;
disp(['refractive index: ', num2str(refractive_index)])
disp(['delta_x: ', num2str(delta_x)])
disp(['delta_z: ', num2str(z_val(2) - z_val(1))])
disp(['x FOV: ', num2str((x_val(end) - x_val(1)) * 1E6), ' microns'])
disp(['z FOV: ', num2str((z_max - z_min) * 1E6), ' microns'])
%% Pupil field
pupil_obliqueness = 1 ./ cos(theta);
pupil_polarisation(:, :, 1) = 1 / 2 * ((1 + cos(theta)) - (1 - cos(theta)) .* cos(2 * phi));
pupil_polarisation(:, :, 2) = -1 / 2 * (1 - cos(theta)) .* sin(2 * phi);
pupil_polarisation(:, :, 3) = -cos(phi) .* sin(theta);
pupil_phase = ones(size(pupil_obliqueness));
% Stop constraint
pupil_phase = pupil_phase .* stop;
pupil_polarisation = pupil_polarisation .* repmat(stop, [1, 1, 3]);
pupil_obliqueness = pupil_obliqueness .* stop;
%% Pupil amplitudes
r = R / max(m) * max_numerical_aperture;
intensity = zeros([size(pupil_obliqueness), num_slices, length(nas)]);
for i=1:(length(nas) - 1)
pupil_amplitude = zeros(size(pupil_obliqueness));
pupil_amplitude(r > nas(i) & r < nas(i+1)) = 1;
pupil_amplitude = pupil_amplitude .* stop;
% Focal field calculations
electric_field = zeros([size(pupil_polarisation), num_slices]);
progress = waitbar(0, ['Calculating for zone ', num2str(i), ' of ' num2str(num_steps), '...']);
tic;
for k_index = 1:num_slices
electric_field(:, :, :, k_index) = debye.propagate(M, N, k_z, z_val(k_index), pupil_obliqueness, pupil_polarisation, pupil_amplitude, pupil_phase);
intensity(:, :, k_index, i) = abs(electric_field(:, :, 1, k_index)).^2 + abs(electric_field(:, :, 2, k_index)).^2 + abs(electric_field(:, :, 3, k_index)).^2;
intensity(:, :, k_index, i) = intensity(:, :, k_index, i).^multiphoton_order;
waitbar(k_index / num_slices)
end
toc
disp(['For ' num2str(num_slices) ' slices.'])
close(progress)
end
%% Conventional focus
pupil_amplitude = zeros(size(pupil_obliqueness));
pupil_amplitude(r < nas(end)) = 1;
pupil_amplitude = pupil_amplitude .* stop;
electric_field = zeros([size(pupil_polarisation), num_slices]);
conventional_intensity = zeros([size(pupil_obliqueness), num_slices]);
progress = waitbar(0, 'Calculating conventional focus...');
tic;
for k_index = 1:num_slices
electric_field(:, :, :, k_index) = debye.propagate(M, N, k_z, z_val(k_index), pupil_obliqueness, pupil_polarisation, pupil_amplitude, pupil_phase);
conventional_intensity(:, :, k_index) = abs(electric_field(:, :, 1, k_index)).^2 + abs(electric_field(:, :, 2, k_index)).^2 + abs(electric_field(:, :, 3, k_index)).^2;
conventional_intensity(:, :, k_index) = conventional_intensity(:, :, k_index).^multiphoton_order;
waitbar(k_index / num_slices)
end
toc
disp(['For ' num2str(num_slices) ' slices.'])
close(progress)
%% Saving
if SAVE_IMAGES
intensity_16bit = uint16(intensity / max(intensity(:)) * 2^16 - 1);
for i = 1:(length(nas) - 1)
for k_index = 1:num_slices
imwrite(intensity_16bit(:, :, k_index, i), ['psfs/', ROOT_NAME, '-' num2str(i), '-', num2str(k_index, '%03d'), '-MP', num2str(multiphoton_order), '.tif'], 'tiff');
end
end
conventional_intensity_16bit = uint16(conventional_intensity / max(conventional_intensity(:)) * 2^16 - 1);
for k_index = 1:num_slices
imwrite(conventional_intensity_16bit(:, :, k_index), ['psfs/', ROOT_NAME, '-' num2str(i), '-', num2str(k_index, '%03d'), '-MP', num2str(multiphoton_order), '-conv.tif'], 'tiff');
end
end
%% Display
if SHOW_IMAGES
% Zone PSFs
for i=1:num_steps
subplot(4, num_steps + 2, i)
imshow(intensity(:, :, floor(size(intensity, 3) / 2), i), [0, max(intensity(:))])
title(['Band ', num2str(i), ' XY'])
subplot(4, num_steps + 2, num_steps + 2 + i)
imshow(intensity(:, :, floor(size(intensity, 3) / 2), i), [0, max(intensity(:))])
zoom(8)
title(['Band ', num2str(i), ' XY zoomed'])
subplot(4, num_steps + 2, 2 * (num_steps + 2) + i)
imshow(squeeze(intensity(:, floor(size(intensity, 2) / 2), :, i))', [0, max(intensity(:))])
title(['Band ', num2str(i), ' XZ'])
subplot(4, num_steps + 2, 3 * (num_steps + 2) + i)
imshow(squeeze(intensity(:, floor(size(intensity, 2) / 2), :, i))', [0, max(intensity(:))])
zoom(8)
title(['Band ', num2str(i), ' XZ zoomed'])
end
% Effective PSF
subplot(4, num_steps + 2, num_steps + 1)
imshow(sum(intensity(:, :, floor(size(intensity, 3) / 2), :), 4), [])
title('Overall XY')
subplot(4, num_steps + 2, 2 * (num_steps + 2) - 1)
imshow(sum(intensity(:, :, floor(size(intensity, 3) / 2), :), 4), [])
title('Overall XY')
subplot(4, num_steps + 2, 3 * (num_steps + 2) - 1)
imshow(sum(squeeze(intensity(:, floor(size(intensity, 2) / 2), :, :)), 3)', [])
title('Overall XZ')
subplot(4, num_steps + 2, 4 * (num_steps + 2) - 1)
imshow(sum(squeeze(intensity(:, floor(size(intensity, 2) / 2), :, :)), 3)', [])
zoom(8)
title('Overall XZ zoomed')
% Conventional PSF
subplot(4, num_steps + 2, num_steps + 2)
imshow(conventional_intensity(:, :, floor(size(conventional_intensity, 3) / 2)), [])
title('Conventional XY')
subplot(4, num_steps + 2, 2 * (num_steps + 2))
imshow(conventional_intensity(:, :, floor(size(conventional_intensity, 3) / 2)), [])
title('Conventional XY')
subplot(4, num_steps + 2, 3 * (num_steps + 2))
imshow(squeeze(conventional_intensity(:, floor(size(conventional_intensity, 2) / 2), :))', [])
title('Conventional XZ')
subplot(4, num_steps + 2, 4 * (num_steps + 2))
imshow(squeeze(conventional_intensity(:, floor(size(conventional_intensity, 2) / 2), :))', [])
zoom(8)
title('Conventional XZ zoomed')
end