Project

General

Profile

root / branches / compiler / cSharp / ooasCompiler / src / libs / c5 / UserGuideExamples / GConvexHull.cs @ 3

1
/*
2
 Copyright (c) 2003-2006 Niels Kokholm and Peter Sestoft
3
 Permission is hereby granted, free of charge, to any person obtaining a copy
4
 of this software and associated documentation files (the "Software"), to deal
5
 in the Software without restriction, including without limitation the rights
6
 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
7
 copies of the Software, and to permit persons to whom the Software is
8
 furnished to do so, subject to the following conditions:
9
 
10
 The above copyright notice and this permission notice shall be included in
11
 all copies or substantial portions of the Software.
12
 
13
 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
14
 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
15
 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
16
 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
17
 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
18
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
19
 SOFTWARE.
20
*/
21

    
22
// Compile with 
23
//    csc /r:C5.dll GConvexHull.cs
24

    
25
using System;
26
using C5;
27

    
28
namespace GConvexHull
29
{
30
// Find the convex hull of a point set in the plane
31

    
32
// An implementation of Graham's (1972) point elimination algorithm,
33
// as modified by Andrew (1979) to find lower and upper hull separately.
34

    
35
// This implementation correctly handle duplicate points, and
36
// multiple points with the same x-coordinate.
37

    
38
// 1. Sort the points lexicographically by increasing (x,y), thus 
39
//    finding also a leftmost point L and a rightmost point R.
40
// 2. Partition the point set into two lists, upper and lower, according as 
41
//    point is above or below the segment LR.  The upper list begins with 
42
//    L and ends with R; the lower list begins with R and ends with L.
43
// 3. Traverse the point lists clockwise, eliminating all but the extreme
44
//    points (thus eliminating also duplicate points).
45
// 4. Join the point lists (in clockwise order) in an array, 
46
//    leaving out L from lower and R from upper.
47

    
48
  public class Convexhull
49
  {
50
    public static Point[] ConvexHull(Point[] pts)
51
    {
52
      // 1. Sort points lexicographically by increasing (x, y)
53
      int N = pts.Length;
54
      Array.Sort(pts);
55
      Point left = pts[0], right = pts[N - 1];
56
      // 2. Partition into lower hull and upper hull
57
      IList<Point> lower = new LinkedList<Point>(),
58
        upper = new LinkedList<Point>();
59
      lower.InsertFirst(left); upper.InsertLast(left);
60
      for (int i = 0; i < N; i++)
61
      {
62
        double det = Point.Area2(left, right, pts[i]);
63
        if (det < 0)
64
          lower.InsertFirst(pts[i]);
65
        else if (det > 0)
66
          upper.InsertLast(pts[i]);
67
      }
68
      lower.InsertFirst(right);
69
      upper.InsertLast(right);
70
      // 3. Eliminate points not on the hull
71
      Eliminate(lower);
72
      Eliminate(upper);
73
      // 4. Join the lower and upper hull, leaving out lower.Last and upper.Last
74
      Point[] res = new Point[lower.Count + upper.Count - 2];
75
      lower[0, lower.Count - 1].CopyTo(res, 0);
76
      upper[0, upper.Count - 1].CopyTo(res, lower.Count - 1);
77
      return res;
78
    }
79

    
80
    // Graham's scan
81
    public static void Eliminate(IList<Point> lst)
82
    {
83
      IList<Point> view = lst.View(0, 0);
84
      int slide = 0;
85
      while (view.TrySlide(slide, 3))
86
        if (Point.Area2(view[0], view[1], view[2]) < 0)   // right turn
87
          slide = 1;
88
        else
89
        {                                                 // left or straight
90
          view.RemoveAt(1);
91
          slide = view.Offset != 0 ? -1 : 0;
92
        }
93
    }
94
  }
95

    
96
// ------------------------------------------------------------
97

    
98
// Points in the plane
99

    
100
  public class Point : IComparable<Point>
101
  {
102
    private static readonly C5Random rnd = new C5Random(42);
103

    
104
    public readonly double x, y;
105

    
106
    public Point(double x, double y)
107
    {
108
      this.x = x; this.y = y;
109
    }
110

    
111
    public override string ToString()
112
    {
113
      return "(" + x + ", " + y + ")";
114
    }
115

    
116
    public static Point Random(int w, int h)
117
    {
118
      return new Point(rnd.Next(w), rnd.Next(h));
119
    }
120

    
121
    public bool Equals(Point p2)
122
    {
123
      return x == p2.x && y == p2.y;
124
    }
125

    
126
    public int CompareTo(Point p2)
127
    {
128
      int major = x.CompareTo(p2.x);
129
      return major != 0 ? major : y.CompareTo(p2.y);
130
    }
131

    
132
    // Twice the signed area of the triangle (p0, p1, p2)
133
    public static double Area2(Point p0, Point p1, Point p2)
134
    {
135
      return p0.x * (p1.y - p2.y) + p1.x * (p2.y - p0.y) + p2.x * (p0.y - p1.y);
136
    }
137
  }
138

    
139
// ------------------------------------------------------------
140

    
141
  class GConvexHull
142
  {
143
    static void Main(String[] args)
144
    {
145
      if (args.Length == 1)
146
      {
147
        string arg = args[0];
148
        int N = int.Parse(arg);
149
        Point[] pts = new Point[N];
150
        for (int i = 0; i < N; i++)
151
          pts[i] = Point.Random(500, 500);
152
        Point[] chpts = Convexhull.ConvexHull(pts);
153
        Console.WriteLine("Area is " + Area(chpts));
154
        Print(chpts);
155
      }
156
      else
157
        Console.WriteLine("Usage: GConvexHull <pointcount>\n");
158
    }
159

    
160
    // The area of a polygon (represented by an array of ordered vertices)
161
    public static double Area(Point[] pts)
162
    {
163
      int N = pts.Length;
164
      Point origo = new Point(0, 0);
165
      double area2 = 0;
166
      for (int i = 0; i < N; i++)
167
        area2 += Point.Area2(origo, pts[i], pts[(i + 1) % N]);
168
      return Math.Abs(area2 / 2);
169
    }
170

    
171
    public static void Print(Point[] pts)
172
    {
173
      int N = pts.Length;
174
      for (int i = 0; i < N; i++)
175
        Console.WriteLine(pts[i]);
176
    }
177
  }
178
}