[37] | 1 | function X = vgg_intersect_quadrics(A, B, C);
|
---|
| 2 | %
|
---|
| 3 | % function X = vgg_intersect_quadrics(A, B, C);
|
---|
| 4 | %
|
---|
| 5 | % Purpose:
|
---|
| 6 | % Compute intersections of three quadrics in projective 3-space.
|
---|
| 7 | %
|
---|
| 8 | % Input:
|
---|
| 9 | % A, B, C are symmetric 4x4 matrices.
|
---|
| 10 | %
|
---|
| 11 | % Output:
|
---|
| 12 | % X is 4-by-8 complex. Each column is a root.
|
---|
| 13 | %
|
---|
| 14 | % The function works by converting the problem into an 8x8 generalized
|
---|
| 15 | % eigen-problem.
|
---|
| 16 | %
|
---|
| 17 | % fsm@robots.ox.ac.uk
|
---|
| 18 | %
|
---|
| 19 |
|
---|
| 20 | if any(size(A) ~= [4 4])
|
---|
| 21 | error('bad A')
|
---|
| 22 | end
|
---|
| 23 | if any(size(B) ~= [4 4])
|
---|
| 24 | error('bad B')
|
---|
| 25 | end
|
---|
| 26 | if any(size(C) ~= [4 4])
|
---|
| 27 | error('bad C')
|
---|
| 28 | end
|
---|
| 29 |
|
---|
| 30 | % symmetrize and normalize.
|
---|
| 31 | A = A + A.'; A = A / norm(A, 'fro'); a = [A(1, 1) 2*A(1, 2) 2*A(1, 3) 2*A(1, 4) A(2, 2) 2*A(2, 3) 2*A(2, 4) A(3, 3) 2*A(3, 4) A(4, 4)].';
|
---|
| 32 | B = B + B.'; B = B / norm(B, 'fro'); b = [B(1, 1) 2*B(1, 2) 2*B(1, 3) 2*B(1, 4) B(2, 2) 2*B(2, 3) 2*B(2, 4) B(3, 3) 2*B(3, 4) B(4, 4)].';
|
---|
| 33 | C = C + C.'; C = C / norm(C, 'fro'); c = [C(1, 1) 2*C(1, 2) 2*C(1, 3) 2*C(1, 4) C(2, 2) 2*C(2, 3) 2*C(2, 4) C(3, 3) 2*C(3, 4) C(4, 4)].';
|
---|
| 34 |
|
---|
| 35 | % make multiplication tables (for degrees 1+2=3, 2+2=4 and 3+1=4)
|
---|
| 36 | if 1
|
---|
| 37 | ijk = [
|
---|
| 38 | 1
|
---|
| 39 | 22
|
---|
| 40 | 43
|
---|
| 41 | 64
|
---|
| 42 | 82
|
---|
| 43 | 105
|
---|
| 44 | 126
|
---|
| 45 | 147
|
---|
| 46 | 163
|
---|
| 47 | 186
|
---|
| 48 | 208
|
---|
| 49 | 229
|
---|
| 50 | 244
|
---|
| 51 | 267
|
---|
| 52 | 289
|
---|
| 53 | 310
|
---|
| 54 | 325
|
---|
| 55 | 351
|
---|
| 56 | 372
|
---|
| 57 | 393
|
---|
| 58 | 406
|
---|
| 59 | 432
|
---|
| 60 | 454
|
---|
| 61 | 475
|
---|
| 62 | 487
|
---|
| 63 | 513
|
---|
| 64 | 535
|
---|
| 65 | 556
|
---|
| 66 | 568
|
---|
| 67 | 594
|
---|
| 68 | 617
|
---|
| 69 | 638
|
---|
| 70 | 649
|
---|
| 71 | 675
|
---|
| 72 | 698
|
---|
| 73 | 719
|
---|
| 74 | 730
|
---|
| 75 | 756
|
---|
| 76 | 779
|
---|
| 77 | 800
|
---|
| 78 | ];
|
---|
| 79 | X12 = zeros(20, 4, 10); X12(ijk) = 1;
|
---|
| 80 |
|
---|
| 81 | ijk = [
|
---|
| 82 | 1
|
---|
| 83 | 37
|
---|
| 84 | 73
|
---|
| 85 | 109
|
---|
| 86 | 145
|
---|
| 87 | 181
|
---|
| 88 | 217
|
---|
| 89 | 253
|
---|
| 90 | 289
|
---|
| 91 | 325
|
---|
| 92 | 352
|
---|
| 93 | 390
|
---|
| 94 | 426
|
---|
| 95 | 462
|
---|
| 96 | 501
|
---|
| 97 | 537
|
---|
| 98 | 573
|
---|
| 99 | 609
|
---|
| 100 | 645
|
---|
| 101 | 681
|
---|
| 102 | 703
|
---|
| 103 | 741
|
---|
| 104 | 778
|
---|
| 105 | 814
|
---|
| 106 | 852
|
---|
| 107 | 889
|
---|
| 108 | 925
|
---|
| 109 | 962
|
---|
| 110 | 998
|
---|
| 111 | 1034
|
---|
| 112 | 1054
|
---|
| 113 | 1092
|
---|
| 114 | 1129
|
---|
| 115 | 1165
|
---|
| 116 | 1203
|
---|
| 117 | 1240
|
---|
| 118 | 1276
|
---|
| 119 | 1313
|
---|
| 120 | 1349
|
---|
| 121 | 1385
|
---|
| 122 | 1405
|
---|
| 123 | 1446
|
---|
| 124 | 1482
|
---|
| 125 | 1518
|
---|
| 126 | 1561
|
---|
| 127 | 1597
|
---|
| 128 | 1633
|
---|
| 129 | 1669
|
---|
| 130 | 1705
|
---|
| 131 | 1741
|
---|
| 132 | 1756
|
---|
| 133 | 1797
|
---|
| 134 | 1834
|
---|
| 135 | 1870
|
---|
| 136 | 1912
|
---|
| 137 | 1949
|
---|
| 138 | 1985
|
---|
| 139 | 2022
|
---|
| 140 | 2058
|
---|
| 141 | 2094
|
---|
| 142 | 2107
|
---|
| 143 | 2148
|
---|
| 144 | 2185
|
---|
| 145 | 2221
|
---|
| 146 | 2263
|
---|
| 147 | 2300
|
---|
| 148 | 2336
|
---|
| 149 | 2373
|
---|
| 150 | 2409
|
---|
| 151 | 2445
|
---|
| 152 | 2458
|
---|
| 153 | 2499
|
---|
| 154 | 2537
|
---|
| 155 | 2573
|
---|
| 156 | 2614
|
---|
| 157 | 2652
|
---|
| 158 | 2688
|
---|
| 159 | 2726
|
---|
| 160 | 2762
|
---|
| 161 | 2798
|
---|
| 162 | 2809
|
---|
| 163 | 2850
|
---|
| 164 | 2888
|
---|
| 165 | 2924
|
---|
| 166 | 2965
|
---|
| 167 | 3003
|
---|
| 168 | 3039
|
---|
| 169 | 3077
|
---|
| 170 | 3113
|
---|
| 171 | 3149
|
---|
| 172 | 3160
|
---|
| 173 | 3201
|
---|
| 174 | 3239
|
---|
| 175 | 3275
|
---|
| 176 | 3316
|
---|
| 177 | 3354
|
---|
| 178 | 3390
|
---|
| 179 | 3428
|
---|
| 180 | 3464
|
---|
| 181 | 3500
|
---|
| 182 | ];
|
---|
| 183 | X22 = zeros(35, 10, 10); X22(ijk) = 1;
|
---|
| 184 |
|
---|
| 185 | ijk = [
|
---|
| 186 | 1
|
---|
| 187 | 37
|
---|
| 188 | 73
|
---|
| 189 | 109
|
---|
| 190 | 145
|
---|
| 191 | 181
|
---|
| 192 | 217
|
---|
| 193 | 253
|
---|
| 194 | 289
|
---|
| 195 | 325
|
---|
| 196 | 361
|
---|
| 197 | 397
|
---|
| 198 | 433
|
---|
| 199 | 469
|
---|
| 200 | 505
|
---|
| 201 | 541
|
---|
| 202 | 577
|
---|
| 203 | 613
|
---|
| 204 | 649
|
---|
| 205 | 685
|
---|
| 206 | 702
|
---|
| 207 | 740
|
---|
| 208 | 776
|
---|
| 209 | 812
|
---|
| 210 | 851
|
---|
| 211 | 887
|
---|
| 212 | 923
|
---|
| 213 | 959
|
---|
| 214 | 995
|
---|
| 215 | 1031
|
---|
| 216 | 1071
|
---|
| 217 | 1107
|
---|
| 218 | 1143
|
---|
| 219 | 1179
|
---|
| 220 | 1215
|
---|
| 221 | 1251
|
---|
| 222 | 1287
|
---|
| 223 | 1323
|
---|
| 224 | 1359
|
---|
| 225 | 1395
|
---|
| 226 | 1403
|
---|
| 227 | 1441
|
---|
| 228 | 1478
|
---|
| 229 | 1514
|
---|
| 230 | 1552
|
---|
| 231 | 1589
|
---|
| 232 | 1625
|
---|
| 233 | 1662
|
---|
| 234 | 1698
|
---|
| 235 | 1734
|
---|
| 236 | 1772
|
---|
| 237 | 1809
|
---|
| 238 | 1845
|
---|
| 239 | 1882
|
---|
| 240 | 1918
|
---|
| 241 | 1954
|
---|
| 242 | 1991
|
---|
| 243 | 2027
|
---|
| 244 | 2063
|
---|
| 245 | 2099
|
---|
| 246 | 2104
|
---|
| 247 | 2142
|
---|
| 248 | 2179
|
---|
| 249 | 2215
|
---|
| 250 | 2253
|
---|
| 251 | 2290
|
---|
| 252 | 2326
|
---|
| 253 | 2363
|
---|
| 254 | 2399
|
---|
| 255 | 2435
|
---|
| 256 | 2473
|
---|
| 257 | 2510
|
---|
| 258 | 2546
|
---|
| 259 | 2583
|
---|
| 260 | 2619
|
---|
| 261 | 2655
|
---|
| 262 | 2692
|
---|
| 263 | 2728
|
---|
| 264 | 2764
|
---|
| 265 | 2800
|
---|
| 266 | ];
|
---|
| 267 | X31 = zeros(35, 20, 4); X31(ijk) = 1;
|
---|
| 268 |
|
---|
| 269 | clear ijk
|
---|
| 270 | end
|
---|
| 271 |
|
---|
| 272 | % compose ideal in degree 2.
|
---|
| 273 | I2 = [a b c];
|
---|
| 274 |
|
---|
| 275 | % compute ideal in degrees 3 and 4.
|
---|
| 276 | I3 = reshape(reshape(X12, [20* 4 10]) * I2, [20 3* 4]);
|
---|
| 277 | I4 = reshape(reshape(X22, [35*10 10]) * I2, [35 3*10]);
|
---|
| 278 |
|
---|
| 279 | % the columns of C3 and C4 span complements of I3 and I4, respectively.
|
---|
| 280 | [U,S,V] = svd(I3.'); C3 = V(:, 13:20); %W3 = diag(S).'
|
---|
| 281 | [U,S,V] = svd(I4.'); C4 = V(:, 28:35); %W4 = diag(S).'
|
---|
| 282 |
|
---|
| 283 | % make the four multiplication operators.
|
---|
| 284 | M = cell(1, 4);
|
---|
| 285 | for k=1:4
|
---|
| 286 | M{k} = C4.' * X31(:, :, k) * C3;
|
---|
| 287 | end
|
---|
| 288 |
|
---|
| 289 | % FIXME: here we only use M{1} and M{2} to compute the Schur
|
---|
| 290 | % decomposition although we use all four multiplication operators to
|
---|
| 291 | % extract the root coordinates.
|
---|
| 292 | [AA, BB, Q, Z] = qz(M{1}, M{2});
|
---|
| 293 | for k=1:4
|
---|
| 294 | M{k} = Q*M{k}*Z;
|
---|
| 295 | end
|
---|
| 296 | %M{:}
|
---|
| 297 |
|
---|
| 298 | X = zeros(4, 8);
|
---|
| 299 | X(1, :) = diag(M{1}).';
|
---|
| 300 | X(2, :) = diag(M{2}).';
|
---|
| 301 | X(3, :) = diag(M{3}).';
|
---|
| 302 | X(4, :) = diag(M{4}).';
|
---|
| 303 |
|
---|
| 304 | % normalize
|
---|
| 305 | X = X ./ repmat(max(abs(X), [], 1), [4 1]);
|
---|
| 306 |
|
---|
| 307 | return;
|
---|
| 308 |
|
---|