-
Notifications
You must be signed in to change notification settings - Fork 0
/
Calculate_Calibration.m
182 lines (144 loc) · 6.16 KB
/
Calculate_Calibration.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
function [coefficients, fval, exitflag, output] = Calculate_Calibration( ...
phantom_img, spacing, phantom_radius, n_circle_steps, circle_width, ...
options)
%% Calculate_Calibration - calculates the FDC
% Algorithm accepts phantom image, its spacing, radius and an optimset and
% returns the coefficients for FDC after calculating them.
% Note: Execution might take up to 15 min, depending on image size and
% options.
%
% INPUTS:
% * phantom_img : [k x n x m] array
% OCT image data of spherical phantom. Dimensions are to be denoted
% as [Z x X x Y], with Z being the dimension along the optical axis
% (depth).
% * spacing : [1 x 3] array
% Spacing of the OCT imaging device ([Z x X x Y]).
% * phantom_radius : float
% REAL radius of spherical phantom.
% * options : [1x1] struct
% Optimization options for fminsearch-function. For more information,
% search documentation on 'fminsearch' and 'optimset'.
%
% OUTPUTS:
% * coefficients : struct
% Coefficients for usage in Apply_Coefficients_To_Surface.
% * fval : float
% Final loss function value in the end of optimization process.
% * exitflag : int
% Exitflag for abortion criteria. For more information,
% search documentation on 'fminsearch'.
% * output : [1x1] struct
% Information about the optimization process. For more information,
% search documentation on 'fminsearch'.
%
% DEPENDENCIES / TOOLBOXES
% * Apply_Coefficients_To_Surface.m
% * The MATLAB Optimization Toolbox is required for the execution of
% this file.
%
% Author: Maron Dolling
% Institute for Biomedical Optics - Universitaet zu Luebeck
% and
% Medical Laser Center Luebeck
% Email: [email protected]
% January 2023
%------------------------- START MAIN -------------------------------------
n = length(phantom_img);
%% Initialize start values for field correction coefficients
% Note: coefficients will later be saved in a struct with notation in
% paper-correspondence
coefficients_initial = [
0,1,0,0,... % in x direction
0,1,0,0,... % in y direction
1,0,0, ... % z-dependent scaling of x- and y-direction
0,0,0, ... % |
0,0,0, ... % |
0,0,0, ... % | => Zernike surface distortion approximation
0,0,0, ... % |
0,0,0, ... % |
];
% Add small random value
coefficients_initial = coefficients_initial + ...
(rand(1, numel(coefficients_initial))-0.5) .* 1e-3;
%% Prepare OCT data
phantom_img_surface = cell(1,n);
for i = 1:n
phantom_img_surface{i} = Surface_Detection_Phantom(phantom_img{i}, ...
spacing);
end
%% Start optimization process
[optim_out, fval, exitflag, output] = fminsearch( @(c) ...
loss_function(phantom_img_surface, phantom_radius, c, ...
n_circle_steps, circle_width), coefficients_initial, options);
%% Make coefficients a struct with corresponding notation
c.q10 = optim_out(1); c.q11 = optim_out(2); c.q12 = optim_out(3); c.q13 = optim_out(4);
c.q20 = optim_out(5); c.q21 = optim_out(6); c.q22 = optim_out(7); c.q23 = optim_out(8);
c.s01 = optim_out(9); c.s02 = optim_out(10); c.s03 = optim_out(11);
c.c0 = optim_out(12); c.c1 = optim_out(13); c.c2 = optim_out(14);
c.c3 = optim_out(15); c.c4 = optim_out(16); c.c5 = optim_out(17);
c.c6 = optim_out(18); c.c7 = optim_out(19); c.c8 = optim_out(20);
c.c9 = optim_out(21); c.c10 = optim_out(22); c.c11 = optim_out(23);
c.c12 = optim_out(24); c.c13 = optim_out(25); c.c14 = optim_out(26);
coefficients = c;
end
%------------------------- END MAIN ---------------------------------------
%% Loss function for Optimization prcedure
function loss = loss_function(phantom_surfaces, phantom_radius, ...
coefficients, n_circle_steps, circle_width)
% Calculates the mean deviation of radii at certain angles around the
% phantom from the expected.
n = length(phantom_surfaces);
radii = nan(n, n_circle_steps);
for j = 1:n
% Apply the coefficients to correct surface
phantom_surface_corrected = Apply_Coefficients_To_Surface( ...
phantom_surfaces{j}, coefficients);
x = phantom_surface_corrected(:,1);
y = phantom_surface_corrected(:,2);
z = phantom_surface_corrected(:,3);
% Remove any nan-values from data
ii = and(and(~isnan(x), ~isnan(y)), ~isnan(z));
x = x(ii);
y = y(ii);
z = z(ii);
% Find BFS and center the phantom in x-y-plane
[center, ~] = ellipsoid_fit([x,y,z], 'xyz');
x = x - center(1);
y = y - center(2);
z = z - center(3);
% Get all the angles the circle is fittet at
angles = linspace(-pi/2, pi/2 , n_circle_steps);
% add small random, so angles are not the same in every iteration
angles = angles + (rand * pi / n_circle_steps);
% Handle angles that exceed the limit
angles(angles > pi/2) = angles(angles > pi/2) - pi;
angles = sort(angles);
% Calculate the radius for every angle
for i = 1:n_circle_steps
% Calculate the two points l1, l2 that define the line that
% encapsules all points with a distance smaller than "width"
l1 = repmat([cos(angles(i)), sin(angles(i)), 0], length(x), 1);
l2 = -l1;
distances = point_to_line_distance([x,y,zeros(length(x),1)], l1, l2);
% Get indices of points within cirlce widt
ii = find(distances < circle_width / 2);
% Skip sphere fit for very sparse surfaces
if length(ii) < 4
continue
end
% Calculate radius of circle at angle
[~, r] = ellipsoid_fit([x(ii), y(ii), z(ii)], 'xyz');
radii(j,i) = r(1);
end % for semicircles
end % for surfaces
% Calculate loss as mean deviation of radii form expected
loss = mean((radii-phantom_radius).^2, 'all');
end
function d = point_to_line_distance(pt, l1, l2)
% Calculates the distance of a point pt to a line crossing points l1 and l2
a = l1 - l2;
b = pt - l2;
c = cross(a, b, 2);
d = sqrt(sum(c.^2, 2)) ./ sqrt(sum(a.^2, 2));
end