Classroom Examples of Robustness Problems in Geometric Computations

1
Classroom Examples of Robustness Problems in
Geometric Computations

• Sylvain Pion
• INRIA, Sophia Antipolis
• Chee Yap
• New York University

Lutz Kettner, Kurt Mehlhorn MPI für Informatik,
Saarbrücken Stefan Schirra Otto-von-Guericke-Univ
ersität, Magdeburg
2
Overview

• Motivation
• 2D Convex Hull Algorithm
• Constructing Bad Input Examples
• Analysis

3
Motivation

• The algorithms of computational geometry are
  designed for a machine model with exact real
  arithmetic.
• Substituting floating point arithmetic for the
  assumed real arithmetic may cause implementations
  to fail.
• Although this is well known, it is not common
  knowledge.
• There is no good material for the classroom to
  show the inadequacy of floating point arithmetic.

4
Example

• Arrangements of circular arcs
• Seen in Antibes, France, during SGP04

5
Example

• Arrangements of circular arcs, with double
  arithmetic

6
Example

• Arrangements of circular arcs, with float
  arithmetic

7
Motivation

• Lets be more serious

8
Intersection of Four Simple Solids

• Rhino3D
• ((( s1 n s2) n c2) n c1) ? successful
• ((( c1 n c2) n s1) n s2) ? Boolean operation
  failed''

9
Intersection of Four Simple Solids
output is a combi- natorial object
plus coordinates (not a point set)

• Rhino3D
• ((( s1 n s2) n c2) n c1) ? successful
• ((( c1 n c2) n s1) n s2) ? Boolean operation
  failed''

10
Motivation

• Lets be more serious
• , but that is hard to explain in all detail in
  class.
• Lets do a really simple 2D convex hull algorithm.

11
2D-Orientation of Three Points

• Orientation( p, q, r) sign
• Orientation( p, q, r) sign((qx-px)(ry-py)
  (qy-py)(rx-px))
• Implemented with IEEE 784 double precision
• ? float_orient (p, q, r)

12
Convex Hulls in the Plane
v3
v1
r
v2

• maintain current hull as a circular list
  L(v0,v1,...,vk-1) of its extreme points in
  counter-clockwise order
• start with three non-collinear points in S.
• consider the remaining points r one by one.

13
Convex Hulls in the Plane
v4
v1
v3
v2

• maintain current hull as a circular list
  L(v0,v1,...,vk-1) of its extreme points in
  counter-clockwise order
• start with three non-collinear points in S.
• consider the remaining points r one by one.

14
Convex Hulls in the Plane
v4
v1
v3
v2

• maintain current hull as a circular list
  L(v0,v1,...,vk-1) of its extreme points in
  counter-clockwise order
• start with three non-collinear points in S.
• consider the remaining points r one by one.

15
Convex Hulls in the Plane
v4
r
v1
v3
v2

• maintain current hull as a circular list
  L(v0,v1,...,vk-1) of its extreme points in
  counter-clockwise order
• start with three non-collinear points in S.
• consider the remaining points r one by one.

16
Convex Hulls in the Plane
v4
r
v1
v3
v2

• Property A A point r is outside CH iff there is
  an i such that the edge (vi, vi1) is visible for
  r. (orientation(vi,vi1,r) gt 0)

17
Convex Hulls in the Plane
v4
r
v1
v3
v2

• Property A A point r is outside CH iff there is
  an i such that the edge (vi, vi1) is visible for
  r. (orientation(vi,vi1,r) gt 0)

18
Convex Hulls in the Plane
v4
r
v1
v3
v2

• Property A A point r is outside CH iff there is
  an i such that the edge (vi, vi1) is visible for
  r. (orientation(vi,vi1,r) gt 0)

19
Convex Hulls in the Plane
v4
r
v1
v3
v2

• Property A A point r is outside CH iff there is
  an i such that the edge (vi, vi1) is visible for
  r. (orientation(vi,vi1,r) gt 0)
• Property B If r is outside CH, then the set of
  edges that are weakly visible ( orientation is
  non-negative) from r forms a non-empty
  consecutive subchain so does the set of edges
  that are not weakly visible from r.

20
Convex Hulls in the Plane
v4
r
v1
v3
v2

• Property A A point r is outside CH iff there is
  an i such that the edge (vi, vi1) is visible for
  r. (orientation(vi,vi1,r) gt 0)
• Property B If r is outside CH, then the set of
  edges that are weakly visible ( orientation is
  non-negative) from r forms a non-empty
  consecutive subchain so does the set of edges
  that are not weakly visible from r.

21
Convex Hulls in the Plane
v4
r
v1
v3
v2

• Property A A point r is outside CH iff there is
  an i such that the edge (vi, vi1) is visible for
  r. (orientation(vi,vi1,r) gt 0)
• Property B If r is outside CH, then the set of
  edges that are weakly visible ( orientation is
  non-negative) from r forms a non-empty
  consecutive subchain so does the set of edges
  that are not weakly visible from r.

22
Convex Hulls in the Plane
v4
r
v1
v3
v2

• Property A A point r is outside CH iff there is
  an i such that the edge (vi, vi1) is visible for
  r. (orientation(vi,vi1,r) gt 0)
• Property B If r is outside CH, then the set of
  edges that are weakly visible ( orientation is
  non-negative) from r forms a non-empty
  consecutive subchain so does the set of edges
  that are not weakly visible from r.

23
Convex Hulls in the Plane
v4
r
v1
v3
v2

• Property A A point r is outside CH iff there is
  an i such that the edge (vi, vi1) is visible for
  r. (orientation(vi,vi1,r) gt 0)
• Property B If r is outside CH, then the set of
  edges that are weakly visible ( orientation is
  non-negative) from r forms a non-empty
  consecutive subchain so does the set of edges
  that are not weakly visible from r.

24
Convex Hulls in the Plane
v4
r
v1
v3
v2

• Property A A point r is outside CH iff there is
  an i such that the edge (vi, vi1) is visible for
  r. (orientation(vi,vi1,r) gt 0)
• Property B If r is outside CH, then the set of
  edges that are weakly visible ( orientation is
  non-negative) from r forms a non-empty
  consecutive subchain so does the set of edges
  that are not weakly visible from r.

25
Convex Hulls in the Plane
v4
v3
v1
v2

• Property A A point r is outside CH iff there is
  an i such that the edge (vi, vi1) is visible for
  r. (orientation(vi,vi1,r) gt 0)
• Property B If r is outside CH, then the set of
  edges that are weakly visible ( orientation is
  non-negative) from r forms a non-empty
  consecutive subchain so does the set of edges
  that are not weakly visible from r.

26
Devils Advocate

• Systematic construction of instances leading to
  violations of properties (A) and (B) when
  executed with doubles
• and in all possible ways
• a point outside sees no edge
• a point inside sees an edge
• a point outside sees all edges
• a point outside sees a non-contiguous set of
  edges
• examples involve nearly collinear points, of
  course
• examples are realistic as many real-life
  instances contain collinear points (which become
  nearly collinear by conversion to doubles)

27
Global Consequences

• A point outside sees no edge of the current hull.
• p1(7.3000000000000194, 7.3000000000000167)
• p2(24.000000000000068, 24.000000000000071)
• p3(24.00000000000005, 24.000000000000053)
• p4(0.50000000000001621, 0.50000000000001243)
• p5(8, 4) p6(4, 9) p7(15, 27)
• p8(26, 25) p9(19, 11)
• float_orient(p1,p2,p3) gt 0
• float_orient(p1,p2,p4) gt 0
• float_orient(p2,p3,p4) gt 0
• float_orient(p3,p1,p4) gt 0

28
Global Consequences

• A point outside sees all edges of the current
  hull.
• p1(200, 49.200000000000003)
• p2(100, 49.600000000000001)
• p3(-233.33333333333334, 50.93333333333333)
• p4(166.66666666666669, 49.333333333333336)
• float_orient(p1,p2,p3) gt 0
• float_orient(p1,p2,p4) lt 0
• float_orient(p2,p3,p4) lt 0
• float_orient(p3,p1,p4) lt 0

29
Global Consequences

• A point outside sees all edges of the current
  hull.
• p1(200, 49.200000000000003)
• p2(100, 49.600000000000001)
• p3(-233.33333333333334, 50.93333333333333)
• p4(166.66666666666669, 49.333333333333336)
• float_orient(p1,p2,p3) gt 0
• float_orient(p1,p2,p4) lt 0
• float_orient(p2,p3,p4) lt 0
• float_orient(p3,p1,p4) lt 0

• Depending on the implementation
• Algorithm does not terminate!
• Algorithm crashes!

30
Global Consequences

• A point inside sees an edge of the current hull.
• p1(24.00000000000005, 24.000000000000053)
• p2(24.0, 6.0)
• p3(54.85, 6.0)
• p4 (54.850000000000357, 61.000000000000121)
• p5(24.000000000000068, 24.000000000000071)
• (p1,p2,p3,p4) form a convex quadrilateral
• p4 is truly inside this quadrilateral, but
• float_orient(p4,p1,p5) gt 0

p4
p5
p1
p2
p3
31
Global Consequences
p4

• A point inside sees an edge of the current hull.
• p1(24.00000000000005, 24.000000000000053)
• p2(24.0, 6.0)
• p3(54.85, 6.0)
• p4 (54.850000000000357, 61.000000000000121)
• p5(24.000000000000068, 24.000000000000071)
• (p1,p2,p3,p4) form a convex quadrilateral
• p4 is truly inside this quadrilateral, but
• float_orient(p4,p1,p5) gt 0

p5
p1
p4
p2
p5
p1
p2
p3
32
Global Consequences

• A point outside sees a non-contiguous set of
  edges
• p1(24.00000000000005, 24.000000000000053)
• p2(24.0, 6.0)
• p3(54.85, 6.0)
• p4 (54.850000000000357, 61.000000000000121)
• p5(24.000000000000068, 24.000000000000071)
• p6(6, 6)

p4
p5
p1
p2
p3
p6
33
Global Consequences
p4

• A point outside sees a non-contiguous set of
  edges
• p1(24.00000000000005, 24.000000000000053)
• p2(24.0, 6.0)
• p3(54.85, 6.0)
• p4 (54.850000000000357, 61.000000000000121)
• p5(24.000000000000068, 24.000000000000071)
• p6(6, 6)

p5
p1
p4
p2
p5
p1
p2
p3
p6
34
Global Consequences
p4

• A point outside sees a non-contiguous set of
  edges
• p1(24.00000000000005, 24.000000000000053)
• p2(24.0, 6.0)
• p3(54.85, 6.0)
• p4 (54.850000000000357, 61.000000000000121)
• p5(24.000000000000068, 24.000000000000071)
• p6(6, 6)
• float_orient(p4,p5,p6) lt 0
• float_orient(p5,p1,p6) gt 0
• float_orient(p1,p2,p6) lt 0

p5
p1
p4
p2
p5
p1
p2
p3
p6
35
Global Consequences
p4

• A point outside sees a non-contiguous set of
  edges
• p1(24.00000000000005, 24.000000000000053)
• p2(24.0, 10.0)
• p3(54.85, 6.0)
• p4 (54.850000000000357, 61.000000000000121)
• p5(24.000000000000068, 24.000000000000071)
• p6(6, 6)

p5
p1
p4
p2
p5
p1
p2
p3
p6
36
Global Consequences

• A point outside sees a non-contiguous set of
  edges
• p1(24.00000000000005, 24.000000000000053)
• p2(24.0, 10.0)
• p3(54.85, 6.0)
• p4 (54.850000000000357, 61.000000000000121)
• p5(24.000000000000068, 24.000000000000071)
• p6(6, 6)

p4
p5
p1
p2
p3
p6
37
Global Consequences
p4

• A point outside sees a non-contiguous set of
  edges
• p1(24.00000000000005, 24.000000000000053)
• p2(24.0, 10.0)
• p3(54.85, 6.0)
• p4 (54.850000000000357, 61.000000000000121)
• p5(24.000000000000068, 24.000000000000071)
• p6(6, 6)

p5
p1
p6
p4
p5
p1
p2
p3
p6
38
Global Consequences
39
Global Consequences
p4

• A point outside sees a non-contiguous set of
  edges
• p1(24.00000000000005, 24.000000000000053)
• p2(24.0, 6.0)
• p3(54.85, 6.0)
• p4 (54.850000000000357, 61.000000000000121)
• p5(24.000000000000068, 24.000000000000071)
• p6(6, 6)

p5
p1
p4
p2
p5
p1
p2
p3
p6
40
Global Consequences

• A point outside sees a non-contiguous set of
  edges
• p1(24.00000000000005, 24.000000000000053)
• p2(24.0, 6.0)
• p3(54.85, 6.0)
• p4 (54.850000000000357, 61.000000000000121)
• p5(24.000000000000068, 24.000000000000071)
• p6(6, 6)

p4
p5
p1
p2
p3
p6
41
Global Consequences
p4

• A point outside sees a non-contiguous set of
  edges
• p1(24.00000000000005, 24.000000000000053)
• p2(24.0, 6.0)
• p3(54.85, 6.0)
• p4 (54.850000000000357, 61.000000000000121)
• p5(24.000000000000068, 24.000000000000071)
• p6(6, 6)

p6
p5
p1
p4
p6
p5
p1
p2
p3
p6
42
Global Consequences
43
2D-Orientation of Three Points

• Orientation( p, q, r) sign((qx-px)(ry-py)
  (qy-py)(rx-px))

44
2D-Orientation of Three Points

• Orientation( p, q, r) sign((qx-px)(ry-py)
  (qy-py)(rx-px))
• 256x256 pixel image
• redpos., yellow0, blueneg.
• orientation evaluated with double

45
2D-Orientation of Three Points

• Orientation( p, q, r) sign((qx-px)(ry-py)
  (qy-py)(rx-px))
• 256x256 pixel image
• redpos., yellow0, blueneg.
• orientation evaluated with double

46
2D-Orientation of Three Points

• Orientation( p, q, r) sign((qx-px)(ry-py)
  (qy-py)(rx-px))
• Keep y-coordinate fix

47
2D-Orientation of Three Points

• Orientation( p, q, r) sign((qx-px)(ry-py)
  (qy-py)(rx-px))
• Keep y-coordinate fix
• block sizes of 25 and 24

48
2D-Orientation of Three Points

• Orientation( p, q, r) sign((qx-px)(ry-py)
  (qy-py)(rx-px))
• 256x256 pixel image
• redpos., yellow0, blueneg.
• orientation evaluated with double

49
2D-Orientation of Three Points

• Orientation( p, q, r) sign((qx-px)(ry-py)
  (qy-py)(rx-px))
• 256x256 pixel image
• redpos., yellow0, blueneg.
• orientation evaluated with double

50
2D-Orientation of Three Points

• Orientation( p, q, r) sign((qx-px)(ry-py)
  (qy-py)(rx-px))
• 256x256 pixel image
• redpos., yellow0, blueneg.
• orientation evaluated with ext double

51
(Semi-) Systematic Search

• a point outside sees no edge
• p1 ( 0.5, 0.5)
• p2 ( 7.3000000000000194,
• 7.3000000000000167)
• p3 ( 24.000000000000068,
• 24.000000000000071)
• p4 ( 24.00000000000005,
• 24.000000000000053)
• (p2, p3, p4) form a
• counter-clockwise triangle

• Classification of (x(p1) xu, y(p1) yu) with
  respect to the edges (p2, p3) and (p4,p2).
• red sees no edge, ocher collinear with one,
  yellow collinear
• with both, blue sees one but not the other,
  green sees both

52
Choice of the Pivot?
Orientation( p, q, r) sign((qx-px)(ry-py)
(qy-py)(rx-px))
p
53
Choice of the Pivot?
Orientation( p, q, r) sign((qx-px)(ry-py)
(qy-py)(rx-px))
p
q
r
54
Conclusion

• Classroom examples
• Data sets and C programs available online
• http//www.mpi-sb.mpg.de/kettner/proj/NonRobust/
• Long version (available soon)
• More analysis
• Grahams scan algorithm
• 3D Delaunay triangulation