19
|
1 |
== Predicates/Constructions
|
|
2 |
* circumsphere, side of circumsphere
|
|
3 |
* barycentric coordinates
|
|
4 |
* degenerate simplex
|
|
5 |
* projection on a k-flat, squared distance to the k-flat
|
|
6 |
* signed square distance to an oriented (d-1)-flat (?)
|
|
7 |
|
|
8 |
|
|
9 |
== Computing circumcenter and circumradius of a k-simplex in d dimensions.
|
|
10 |
|
|
11 |
* Ken Clarkson's Gram-Schmidt-based algorithm works perfectly over any field.
|
|
12 |
The idea is as follows.
|
|
13 |
|
|
14 |
Given points p_i, the circumcenter c satisfies (p_i-c)^2 = r^2. By
|
|
15 |
translation, we can assume that some p_i is at the origin, so that c^c=r^2,
|
|
16 |
and that equation becomes p_i \cdot c - p_i^2/2 = 0. That is, [c -1/2] is
|
|
17 |
normal to [p_i p_i^2]. Also, c must be an affine combination of the p_i,
|
|
18 |
so if we start with a point [c' 0] with c' in the affine hull of the p_i,
|
|
19 |
and remove all multiples of all [p_i p_i^2] from it, via Gram-Schmidt for
|
|
20 |
example, and then scale the result so that the last entry is -1/2, we get
|
|
21 |
c (alternatively we can start with [0 1]). Note that if the simplex is
|
|
22 |
degenerate (its affine hull is less than k-dimensional), the last entry of
|
|
23 |
c is 0.
|
|
24 |
|
|
25 |
P = {p_0, ..., p_k}
|
|
26 |
p'_i = p_i - p_0 # make p_0 the origin
|
|
27 |
p''_i = [p'_i p'_i^2] # lift the points
|
|
28 |
c' = 0 # c is in the affine hull of p'_i
|
|
29 |
c'' = [c' 1] # c'' is not in the affine hull of p'_i^2
|
|
30 |
|
|
31 |
for i = 1 to k: # create an orthogonal basis for aff(p''_i)
|
|
32 |
a''_i = p''_i
|
|
33 |
for j = k-1 downto 1:
|
|
34 |
a''_i -= a''_j * (a''_i \cdot c''_j)/c''_j^2
|
|
35 |
|
|
36 |
for i = 1 to k:
|
|
37 |
c'' -= a''_i * (c'' \cdot a''_i) / (a''_i^2)
|
|
38 |
c' = c''[1..d]/(-2 * c''[d+1]) # scale c'' down
|
|
39 |
|
|
40 |
* Ulfar Erlingsson's, Erich Kaltofen's and David Musser's "Generic
|
|
41 |
Gram-Schmidt Orthogonalization by Exact Division" (ISSAC'96) gives an
|
|
42 |
algorithm to compute Gram-Schmidt for vector space with integral domain
|
|
43 |
coefficients. I need to understand the paper better to see the advantages
|
|
44 |
over normalizing rationals at every step of the computation above
|
|
45 |
(the paper's introduction says that it avoids "costly GCD computations").
|
|
46 |
|
|
47 |
* Using QR decomposition (or any other orthogonal decomposition), and
|
|
48 |
then using the formula of the circumsphere on the matrix R (which is
|
|
49 |
just the simplex rotated into a k-dimensional space):
|
|
50 |
|
|
51 |
* Only (d-2) determinants of a k by k matrix are required which is nice.
|
|
52 |
|
|
53 |
* All (implemented) computation of the QR decomposition that I have found
|
|
54 |
so far requires square roots (for norms), divisions, and often
|
|
55 |
comparisons. While divisions are acceptable, square roots and
|
|
56 |
comparisons aren't.
|
|
57 |
|
|
58 |
* There is a string of papers in late 80s, early 90s on square root free
|
|
59 |
and division free algorithms for QR decomposition out of VLSI community.
|
|
60 |
There is a paper that claims to put all such techniques into one
|
|
61 |
framework. There is also a paper that says that all such techniques don't
|
|
62 |
work, and square root implementation needs to be reconsidered. (authors
|
|
63 |
to look up if searching for these papers are K.J.R. Liu and E.
|
|
64 |
Frantzeskakis)
|
|
65 |
|
|
66 |
* Using the formulas from Herbert's Weighted Alpha Shapes paper which are
|
|
67 |
derived from formulas of the planes that are lifted Voronoi cells. This
|
|
68 |
approach requires (d-k)*(d+2) determinants of k by k matrices, which is not
|
|
69 |
great.
|