1 | % By Philip Torr 2002
|
---|
2 | % copyright Microsoft Corp.
|
---|
3 | %this function allows for self calibration, it implements the Sturm method given in CVPR 2001
|
---|
4 |
|
---|
5 | %following the steps of Sturm, first provide an f and CC
|
---|
6 | %CC is an estimate of the self calibration parameters other than the focal length
|
---|
7 |
|
---|
8 | %
|
---|
9 | % CC = [aspect_ratio * est_foc, 0 , pp_x;
|
---|
10 | % 0, est_foc , pp_y;
|
---|
11 | % 0, 0, 1/f_est]
|
---|
12 | %
|
---|
13 | % where est_foc is the estimate of the focal length
|
---|
14 | %
|
---|
15 | % G ~ C' F C
|
---|
16 |
|
---|
17 | %to get to the focal length we solve
|
---|
18 | % G ~ diag (1,1,f) E diag (1,1,f)
|
---|
19 |
|
---|
20 |
|
---|
21 |
|
---|
22 | %how far is the focal length from our initial estimate : 1/CC(3,3) ?
|
---|
23 |
|
---|
24 | % well note the focal length we have calculated here is relative to our initial guess
|
---|
25 | %therefore the true focal length is 1/CC(3,3) * focal_length
|
---|
26 | % is we are a factor 5 out then we should start to worry
|
---|
27 |
|
---|
28 | function [focal_length, E, CC_out] = torr_self_calib_f(F,CC)
|
---|
29 |
|
---|
30 |
|
---|
31 | temp_foc = CC(3,3);
|
---|
32 | CC(3,3) = 1;
|
---|
33 |
|
---|
34 | if nargin < 2
|
---|
35 | G = F / norm(F);
|
---|
36 | CC = diag(ones(3,1),0);
|
---|
37 | else
|
---|
38 | %step 3
|
---|
39 | G = CC' * F * CC;
|
---|
40 | G = G / norm(G);
|
---|
41 | end
|
---|
42 |
|
---|
43 | [U,S,V] = svd(G);
|
---|
44 |
|
---|
45 | if abs(S(3,3)) > 0.001
|
---|
46 | error('F must be rank 2 to self calibrate');
|
---|
47 | end
|
---|
48 |
|
---|
49 | %next construct the quadratic equation of Sturm
|
---|
50 | % c1 f^4 + c2 f^2 + c3 = 0
|
---|
51 |
|
---|
52 | a = S(1,1)^2;
|
---|
53 | b = S(2,2)^2;
|
---|
54 |
|
---|
55 |
|
---|
56 | c1 = a * (1 - U(3,1)^2) * (1 - V(3,1)^2) - b * (1 - U(3,2)^2) * (1 - V(3,2)^2);
|
---|
57 | c2 = a *( U(3,1)^2 + V(3,1)^2 - 2 * U(3,1)^2 * V(3,1)^2) - b * (U(3,2)^2 + V(3,2)^2 - 2 * U(3,2)^2 * V(3,2)^2);
|
---|
58 | c3 = a * U(3,1)^2 * V(3,1)^2 - b * U(3,2)^2 * V(3,2)^2;
|
---|
59 |
|
---|
60 |
|
---|
61 | temp = sqrt(c2^2 - 4 * c1 * c3);
|
---|
62 |
|
---|
63 |
|
---|
64 | %first need to check, have the equations degenerated to a linear, i.e is c1 small
|
---|
65 | if (c1^2/(c2^2 + c3^2)) < 0.001 %then linear
|
---|
66 | foc1 = -(U(3,2) * V(3,1) * (S(1,1) * U(3,1) * V(3,1) + S(2,2) * U(3,2) * V(3,2))) ...
|
---|
67 | /(S(1,1) * U(3,1) * U(3,2) * (1 - V(3,1)^2) + (S(2,2) * V(3,1) * V(3,2) * (1 - U(3,2)^2) ));
|
---|
68 | foc2 = -(V(3,2) * U(3,1) * (S(1,1) * U(3,1) * V(3,1) + S(2,2) * U(3,2) * V(3,2)) ) ...
|
---|
69 | /(S(1,1) * V(3,1) * V(3,2) * (1 - U(3,1)^2) + (S(2,2) * U(3,1) * U(3,2) * (1 - V(3,2)^2) ));
|
---|
70 | else
|
---|
71 | foc1 = sqrt((-c2 + temp)/(2 * c1));
|
---|
72 | foc2 = sqrt((-c2 - temp)/(2 * c1));
|
---|
73 | end
|
---|
74 |
|
---|
75 |
|
---|
76 | %we now have two solutions we need to eliminate one. To do this we resort to the linear equations,
|
---|
77 | %simply multiply the two linear equations together and take the minimum absolute value.
|
---|
78 |
|
---|
79 |
|
---|
80 | alg1 = abs( ...
|
---|
81 | foc1^2 * (S(1,1) * U(3,1) * U(3,2) * (1 - V(3,1)^2) + (S(2,2) * V(3,1) * V(3,2) * (1 - U(3,2)^2) ))...
|
---|
82 | + U(3,2) * V(3,1) * (S(1,1) * U(3,1) * V(3,1) + S(2,2) * U(3,2) * V(3,2)) );
|
---|
83 |
|
---|
84 | alg2 = abs( ...
|
---|
85 | foc1^2 * (S(1,1) * V(3,1) * V(3,2) * (1 - U(3,1)^2) + (S(2,2) * U(3,1) * U(3,2) * (1 - V(3,2)^2) ))...
|
---|
86 | + V(3,2) * U(3,1) * (S(1,1) * U(3,1) * V(3,1) + S(2,2) * U(3,2) * V(3,2)) );
|
---|
87 |
|
---|
88 |
|
---|
89 | alg_foc1 = alg1 * alg2;
|
---|
90 |
|
---|
91 |
|
---|
92 |
|
---|
93 | alg1 = abs( ...
|
---|
94 | foc1^2 * (S(1,1) * U(3,1) * U(3,2) * (1 - V(3,1)^2) + (S(2,2) * V(3,1) * V(3,2) * (1 - U(3,2)^2) ))...
|
---|
95 | + U(3,2) * V(3,1) * (S(1,1) * U(3,1) * V(3,1) + S(2,2) * U(3,2) * V(3,2)) );
|
---|
96 |
|
---|
97 | alg2 = abs( ...
|
---|
98 | foc2^2 * (S(1,1) * V(3,1) * V(3,2) * (1 - U(3,1)^2) + (S(2,2) * U(3,1) * U(3,2) * (1 - V(3,2)^2) ))...
|
---|
99 | + V(3,2) * U(3,1) * (S(1,1) * U(3,1) * V(3,1) + S(2,2) * U(3,2) * V(3,2)) );
|
---|
100 |
|
---|
101 |
|
---|
102 |
|
---|
103 |
|
---|
104 | alg_foc2 = alg1 * alg2;
|
---|
105 |
|
---|
106 | if ~isreal(foc1)
|
---|
107 | focal_length = foc2;
|
---|
108 | else
|
---|
109 | if ~isreal(foc2)
|
---|
110 | focal_length = foc1;
|
---|
111 | end
|
---|
112 | end
|
---|
113 |
|
---|
114 |
|
---|
115 |
|
---|
116 |
|
---|
117 | if isreal(foc1) & isreal(foc2)
|
---|
118 |
|
---|
119 | if (alg_foc1 < alg_foc2)
|
---|
120 | focal_length = foc1;
|
---|
121 | else
|
---|
122 | focal_length = foc2;
|
---|
123 | end
|
---|
124 | end
|
---|
125 |
|
---|
126 |
|
---|
127 | %how far is the focal length from our initial estimate : 1/CC(3,3) ?
|
---|
128 |
|
---|
129 | % well note the focal length we have calculated here is relative to our initial guess
|
---|
130 | %therefore the true focal length is 1/CC(3,3) * focal_length
|
---|
131 | % is we are a factor of 5 out then we should start to worry
|
---|
132 | if (~isreal(foc1) & ~isreal(foc2)) | (abs(focal_length ) > 5)| (abs(focal_length) < 0.2)
|
---|
133 | focal_length
|
---|
134 | focal_length = 1;
|
---|
135 | disp('either bad F or needs a stronger calibration method');
|
---|
136 | disp('ooo vicar big fat focal length setting it to original guess');
|
---|
137 | end
|
---|
138 |
|
---|
139 |
|
---|
140 | Cf = [1 0 0; 0 1 0; 0 0 1/focal_length];
|
---|
141 | CC = (Cf * CC);
|
---|
142 |
|
---|
143 | %now we need to get to the essential matrix which is C^t F C = E
|
---|
144 |
|
---|
145 | E = CC' * F * CC;
|
---|
146 |
|
---|
147 | %now remove estimate of the focal length effect
|
---|
148 | focal_length = 1/CC(3,3);
|
---|
149 | focal_length
|
---|
150 | CC_out = CC;
|
---|
151 |
|
---|
152 |
|
---|
153 | [U,S,V] = svd(E);
|
---|
154 | %note that there is a one p[arameter family of SVD's for E
|
---|
155 |
|
---|
156 | if abs(S(3,3)) > 0.00001
|
---|
157 | error('E must be rank 2 to self calibrate');
|
---|
158 | end
|
---|
159 |
|
---|
160 | % if abs(S(1,1) - S(2,2)) > 0.00001
|
---|
161 | % S
|
---|
162 | % % error('E must have two equal singular values');
|
---|
163 | % end
|
---|
164 |
|
---|
165 | %this a problem not pointed out in the Sturm paper, the essential matrix produced might have funny singular values
|
---|
166 | %fix the bugga to have equal singular values:
|
---|
167 |
|
---|
168 | S(3,3) = 0;
|
---|
169 | S(1,1) = S(2,2);
|
---|
170 | E = U*S*V';
|
---|