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