1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
|
void build_hull(deque<pair<double,double>> raw_points,
deque<pair<double,double>> &output)
{
if(raw_points.size()<=4)
{
output = raw_points;
return;
}
output.clear();
//
// The initial array of points is stored in deque raw_points. I first
// sort it, which gives me the far left and far right points of the hull.
// These are special values, and they are stored off separately in the left
// and right members.
//
// I then go through the list of raw_points, and one by one determine whether
// each point is above or below the line formed by the right and left points.
// If it is above, the point is moved into the upper_partition_points sequence. If it
// is below, the point is moved into the lower_partition_points sequence. So the output
// of this routine is the left and right points, and the sorted points that are in the
// upper and lower partitions.
//
std::pair<double,double> left;
std::pair<double,double> right;
std::deque< std::pair<double,double> > upper_partition_points;
std::deque< std::pair<double,double> > lower_partition_points;
//
// Step one in partitioning the points is to sort the raw data
//
std::sort( raw_points.begin(), raw_points.end() );
//
// The the far left and far right points, remove them from the
// sorted sequence and store them in special members
//
left = raw_points.front();
raw_points.pop_front();
//raw_points.erase(raw_points.begin());
right = raw_points.back();
raw_points.pop_back();
//
// Now put the remaining points in one of the two output sequences
//
for ( size_t i = 0 ; i < raw_points.size() ; i++ )
{
int dir = direction( left, right, raw_points[ i ] );
if ( dir < 0 )
upper_partition_points.push_back( raw_points[ i ] );
else
lower_partition_points.push_back( raw_points[ i ] );
}
//std::deque< std::pair<double,double> > lower_hull;
std::deque< std::pair<double,double> > upper_hull;
build_half_hull(lower_partition_points, output, left, right, 1 );
build_half_hull(upper_partition_points, upper_hull, left, right, -1 );
for(deque< pair<double,double> >::reverse_iterator it = upper_hull.rbegin()+1;
it != upper_hull.rend()-1; it++)
{
output.push_back(*it);
}
if(output.size()<=4)
return;
simplifyHull(output,4);
}
//
// Building the hull consists of two procedures: building the lower and
// then the upper hull. The two procedures are nearly identical - the main
// difference between the two is the test for convexity. When building the upper
// hull, our rull is that the middle point must always be *above* the line formed
// by its two closest neighbors. When building the lower hull, the rule is that point
// must be *below* its two closest neighbors. We pass this information to the
// building routine as the last parameter, which is either -1 or 1.
//
// This is the method that builds either the upper or the lower half convex
// hull. It takes as its input a copy of the input array, which will be the
// sorted list of points in one of the two halfs. It produces as output a list
// of the points in the corresponding convex hull.
//
// The factor should be 1 for the lower hull, and -1 for the upper hull.
void build_half_hull(std::deque< std::pair<double,double> > input,
std::deque< std::pair<double,double> > &output,
std::pair<double,double> left,
std::pair<double,double> right,
int factor )
{
//
// The hull will always start with the left point, and end with the right
// point. According, we start by adding the left point as the first point
// in the output sequence, and make sure the right point is the last point
// in the input sequence.
//
input.push_back( right );
output.push_back( left );
//
// The construction loop runs until the input is exhausted
//
while ( input.size() != 0 ) {
//
// Repeatedly add the leftmost point to the null, then test to see
// if a convexity violation has occured. If it has, fix things up
// by removing the next-to-last point in the output suqeence until
// convexity is restored.
//
output.push_back( input.front() );
//plot_hull( f, "adding a new point" );
input.pop_front();
//input.erase( input.begin() );
while ( output.size() >= 3 ) {
size_t end = output.size() - 1;
if ( factor * direction( output[ end - 2 ],
output[ end ],
output[ end - 1 ] ) <= 0 ) {
output.erase( output.begin() + end - 1 );
//plot_hull( f, "backtracking" );
}
else
break;
}
}
}
// We can do this by by translating the points so that p1 is at the origin,
// then taking the cross product of p0 and p2. The result will be positive,
// negative, or 0, meaning respectively that p2 has turned right, left, or
// is on a straight line.
//
static int direction( std::pair<double,double> p0,
std::pair<double,double> p1,
std::pair<double,double> p2 )
{
return ( (p0.first - p1.first ) * (p2.second - p1.second ) )
- ( (p2.first - p1.first ) * (p0.second - p1.second ) );
}
void simplifyHull(std::deque< std::pair<double,double> > &input,int numberOfEdges)
{
int size = input.size();
int i=0;
int j;
double a1,a2,b1,b2;
double x, y, x1, y1, x2, y2;
map<double,int> air;
map<int,pair<double,double>> newPoint;
double surface;
while(size>numberOfEdges)
{
if(findIntersectionPoint(input[i],input[i+1],input[i+2],input[i+3],x,y,surface))
{
air[surface] = i;
newPoint[i] = pair<double,double>(x,y);
}
i++;
if(i == size-3)
{
if(findIntersectionPoint(input[i],input[i+1],input[i+2],input[0],x,y,surface))
{
air[surface] = i;
newPoint[i] = pair<double,double>(x,y);
}
i++;
if(findIntersectionPoint(input[i],input[i+1],input[0],input[1],x,y,surface))
{
air[surface] = i;
newPoint[i] = pair<double,double>(x,y);
}
i++;
if(findIntersectionPoint(input[i],input[0],input[1],input[2],x,y,surface))
{
air[surface] = i;
newPoint[i] = pair<double,double>(x,y);
}
map<double,int>::iterator it = air.begin();
j = it->second;
if(j > size-4)
{
if(j == size-1)
{
input.pop_front();
input.pop_front();
input.push_front(newPoint[j]);
}
else if(j == size-2)
{
input.pop_back();
input.pop_front();
input.push_back(newPoint[j]);
}
else if(j == size-3)
{
input.pop_back();
input.pop_back();
input.push_back(newPoint[j]);
}
}
else
{
input.erase(input.begin()+j+1, input.begin() + j+3);
input.insert(input.begin()+j+1, newPoint[j]);
}
size += -1;
i=0;
air.clear();
}
}
}
bool findIntersectionPoint( pair<double,double> p0,
pair<double,double> p1,
pair<double,double> p2,
pair<double,double> p3,
double &x, double &y, double &surface)
{
double deltax1;
double deltay1;
double deltax2;
double deltay2;
double x1,x2,y1,y2;
double a1,a2,b1,b2;
x1 = p1.first;
y1 = p1.second;
x2 = p2.first;
y2 = p2.second;
deltax1 = x1 - p0.first;
deltay1 = y1-p0.second;
deltax2 = x2 - p3.first;
deltay2 = y2-p3.second;
if(deltax1*deltay1*deltax2*deltay2 == 0)
{
if(deltax1 == 0 && deltax2 == 0) // the 2 segments are vertical
return false;
if(deltay1 == 0 && deltay2 == 0) // the 2 segments are horizontal
return false;
if(deltax1 == 0 && deltay2 == 0)
{
x = x1;
y = y2;
surface = 0.5*abs( (x - x2)*(y - y1) ); // right triangle
return true;
}
if(deltay1 == 0 && deltax2 == 0)
{
x = x2;
y = y1;
surface = 0.5*abs( (x - x1)*(y - y2) ); // right triangle
return true;
}
if(deltax1 == 0)
{
x = x1;
a2 = (y2 - p3.second)/(x2 - p3.first);
b2 = y2-(a2*x2);
y = a2*x + b2;
surface = 0.5*abs( (x2 - x)*(y1 - y) );
return true;
}
if(deltay1 == 0)
{
y = y1;
a2 = (y2 - p3.second)/(x2 - p3.first);
b2 = y2-(a2*x2);
x = (y - b2)/a2;
surface = 0.5*abs( (x1 - x)*(y2 - y) );
return true;
}
if(deltax2 == 0)
{
x = x2;
a1 = (y1 - p0.second)/(x1 - p0.first);
b1 = y1-(a1*x1);
y = a1*x + b1;
surface = 0.5*abs( (x1 - x)*(y2 - y) );
return true;
}
if(deltay2 == 0)
{
y = y2;
a1 = (y1 - p0.second)/(x1 - p0.first);
b1 = y1-(a1*x1);
x = (y - b1)/a1;
surface = 0.5*abs( (x2 - x)*(y1 - y) );
return true;
}
}
else
{
a1 = (y1-p0.second)/(x1 - p0.first);
a2 = (y2-p3.second)/(x2 - p3.first);
b1 = y1-(a1*x1);
b2 = y2-(a2*x2);
x = (b2-b1)/(a1-a2);
y = (a1*x) + b1;
surface = 0.5*abs( ( (x1 - x)*(y2 - y) ) - ( (x2 - x)*(y1 - y) ) );
return true;
}
} |
