-
Notifications
You must be signed in to change notification settings - Fork 1
/
make_sine_plot_4.m
162 lines (137 loc) · 5.44 KB
/
make_sine_plot_4.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
% Make the sine plots. Takes one argument: the (sub?)panel into which to plot.
function make_sine_plot_4(panels, tiffAdj, methods, methodsValid, colours)
%Number of time samples in a full phase:
nsteps = 3333333 / 7980;
%Time
t = linspace(-pi, pi, nsteps);
SHOW_X_GRAPH = 1;
SHOW_Y_GRAPH = 0;
%Beam Velocity
if isempty(panels.children)
panels.pack(1,2 + SHOW_X_GRAPH + SHOW_Y_GRAPH);
end
panels(1,1).select();
x = cos(t);
ti = find(t >= -pi/2 & t <= pi/2);
plot(t(ti), x(ti), 'LineWidth', 2); hold on;
xlabel('Time (phase)');
ylabel('Speed');
%Format
xlim([-1.8, 1.8]);
ylim([0, 1]);
set(gca, 'box', 'off', 'TickDir', 'out');
set(gca, 'XTick', [-pi/2, 0, pi/2], 'XTickLabel', {'-\pi/2', 0, '\pi/2'})
set(gca, 'YTick', [0, 0.5, 1], 'YTickLabel', [0, 0.5, 1])
title('(a) Beam speed');
%Beam position
panels(1,2).select();
% labels = {'Resonant scanner beam position'};
x = sin(t);
ti = find(t >= -pi/2 & t <= pi/2);
plot(t(ti), x(ti), 'LineWidth', 2);
hold on;
% Imaging fraction of beam. Show voxel positions for slower control system for clarity...
D = 0.9;
ControlSystemFreq = 1e6;
ResonantFreq = 7980;
fov = 666; % um for whole FOV
zoomlevel = 1.6;
nsteps_show = ControlSystemFreq / 7980;
t_show = linspace(-pi, pi, nsteps_show);
x_show = sin(t_show);
tt = asin(D);
ti = find(t > -tt & t < tt);
ti_show = find(t_show > -tt & t_show < tt);
%plot(t_show(ti_show),x_show(ti_show), 'k.', 'MarkerSize', 10);
sz = 0.1;
lines_x = [t_show(ti_show)-sz; t_show(ti_show)+sz];
lines_y = [x_show(ti_show); x_show(ti_show)];
line(lines_x, lines_y, 'Color', [0 0 0], 'LineWidth', 0.1);
%hold off;
ylabel('X Position');
xlabel('Time (phase)');
%xlim([-pi,pi]);
legend({'Beam position', 'Voxels'}, 'Location', 'NorthWest', 'box', 'off');
%title('(A) Beam position through time');
title('(b) Voxel positions');
%Format
xlim([-1.8, 1.8]);
ylim([-1, 1]);
set(gca, 'box', 'off', 'TickDir', 'out');
set(gca, 'XTick', [-pi/2, 0, pi/2], 'XTickLabel', {'-\pi/2', 0, '\pi/2'})
set(gca, 'YTick', [-1, 0, 1], 'YTickLabel', {'-\xi', 0, '\xi'})
convert_phase_dist_to_microns = fov/(D*zoomlevel*2);
voxelsizes = diff(x_show(ti_show)) * convert_phase_dist_to_microns;
disp(sprintf('Voxel size: min %g, max %g um', min(voxelsizes), max(voxelsizes)));
if SHOW_X_GRAPH
panels(1,3).select();
% Power compensation
p = cos(t);
pstar = 0.5*(1+cos(t));
%% cos^4 vignetting compensation. Show it for zoom=2.2x, since that's what
%% pstar is calibrated to. That means that the usable part of the X axis has
%% size FOV/zoom, the lens working distance
positions_um = sin(t) * convert_phase_dist_to_microns;
%positions_um(2,:) = sqrt(positions_um(1,:).^2 + 200^2);
%disp(sprintf('Working FOV is %g um', D*(max(positions_um)-min(positions_um))));
%= lens_working_distance * convert_microns_to_phase_dist;
% Divide p by predicted vignetting compensation
%cos3 = cos(atan(positions_um./lens_working_distance)).^3;
cla;
hold on;
for i = find(methodsValid)
centreX = round(size(tiffAdj{i}.p, 1)/2);
centreY = round(size(tiffAdj{i}.p, 2)/2);
plot(asin(tiffAdj{i}.xc / convert_phase_dist_to_microns), ...
tiffAdj{i}.p(:, centreY) / tiffAdj{i}.p(centreX, centreY), ...
'Color', colours(i,:));
end
hold off;
% t(ti), pstar(ti), 'c', ... % Ad-hoc
ylabel('Power');
xlabel('Time (phase)');
% xlim([-pi,pi]);
% ylim([0 1.1]);
legend(methods(find(methodsValid)), 'Location', 'South', 'box', 'off');
%title('(B) Relative voxel size');
title(sprintf('(c) X power compensation, FOV = %d \\mu{}m', ...
round(D*(max(positions_um)-min(positions_um)))));
%Format
xlim([-1.8, 1.8]);
%yl = get(gca, 'YLim');
%ylim([0, 1.3]);
ylim([0.4 1.1]);
set(gca, 'box', 'off', 'TickDir', 'out');
set(gca, 'XTick', [-pi/2, 0,pi/2], 'XTickLabel', {'-\pi/2', 0, '\pi/2'})
set(gca, 'YTick', [0, 0.5, 1, 1.5, 2])
%set(gca, 'YLim', [0.35 2]);
end
if SHOW_Y_GRAPH
panels(1,3 + SHOW_X_GRAPH).select();
cla;
hold on;
for i = 1:length(tiffAdj)
centreX = round(size(tiffAdj{i}.p, 1)/2);
centreY = round(size(tiffAdj{i}.p, 2)/2);
plot(tiffAdj{i}.yc, ...
tiffAdj{i}.p(centreX, :) / tiffAdj{i}.p(centreX, centreY), ...
'Color', colours(i,:));
end
hold off;
% t(ti), pstar(ti), 'c', ... % Ad-hoc
ylabel('Power');
xlabel('Y position (\mu{}m)');
legend(methods, 'Location', 'NorthWest', 'box', 'off');
%title('(B) Relative voxel size');
title(sprintf('(d) Y power compensation, FOV = %d \\mu{}m', ...
round(D*(max(positions_um)-min(positions_um)))));
%Format
%xlim([-1.8, 1.8]);
%yl = get(gca, 'YLim');
%ylim([0, 1.3]);
ylim([0.7 1.8]);
xlim([-290 290]);
set(gca, 'box', 'off', 'TickDir', 'out');
set(gca, 'YTick', [0, 0.5, 1, 1.5, 2])
end
drawnow;