[37] | 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';
|
---|