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 |
|
---|