• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Semi-circle covering

a guest Sep 23rd, 2015 114 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. #include <iostream>
2. #include <iomanip>
3. #include <fstream>
4. #include <array>
5. #include <random>
6. #include <cmath>
7. #include <cassert>
8. #include <limits>
9.
10. typedef std::array<double, 2> vector;
11. typedef std::array<vector, 6> point_list;
12.
13. const auto pi = 4*std::atan(1);
14.
15. double random_uniform(double min, double max)
16. {
17.     static std::default_random_engine generator;
18.     static std::uniform_real_distribution<double> unif_dist; // min = 0, max = 1
19.
20.     return min + (max - min)*unif_dist(generator);
21. }
22.
23. double distance(const vector& a, const vector& b)
24. {
25.     return std::sqrt(std::pow(a[0] - b[0], 2) + std::pow(a[1] - b[1], 2));
26. }
27.
28. bool inside_semicircle(const vector& p)
29. {
30.     return p[1] >= 0 && (distance(p, {0.0, 0.0}) <= 1);
31. }
32.
33. // raondomly move a single point
34. void jostle_point(vector& p, double jostle_distance)
35. {
36.     while(true)
37.     {
38.         auto distance = random_uniform(0, jostle_distance);
39.         auto angle = random_uniform(0, 2*pi);
40.
41.         auto dx = distance*std::cos(angle);
42.         auto dy = distance*std::sin(angle);
43.
44.         auto np = p;
45.         np[0] += dx;
46.         np[1] += dy;
47.         if(inside_semicircle(np))
48.         {
49.             p = np;
50.             break;
51.         }
52.     }
53. }
54.
55. // Move all points randomly and independantly
56. void jostle_points_free(point_list& pl, double jostle_distance)
57. {
58.     for(auto& p : pl)
59.     {
60.         jostle_point(p, jostle_distance);
61.     }
62. }
63.
64. // Write out point list to output stream
65. void write_list(const point_list& pl, std::ostream& os)
66. {
67.     os << "p = [\n";
68.     for(const auto& p : pl)
69.     {
70.         for(auto e : p)
71.         {
72.             os << std::setprecision(std::numeric_limits<double>::max_digits10) << e << ' ';
73.         }
74.         os << ";\n";
75.     }
76.     os << "];\n=====================================================" << std::endl;
77. }
78.
79. // write point list to console and file
80. std::ofstream ofs("result.txt");
81. void write_list(const point_list& pl)
82. {
83.     write_list(pl, std::cout);
84.     write_list(pl, ofs);
85. }
86.
87. // Given a radius associated with points i, j, k
88. // in the point list, check if any othe points
89. // fall within a radius of the given point (usually
90. // the circumcenter of i, j, k).
91. bool is_delaunay(const point_list& pl,
92.                  const vector& point,
94.                  size_t i, size_t j, size_t k = -1)
95. {
96.     if( ! std::isfinite(radius))
97.     {
98.         return false;
99.     }
100.
101.     for(size_t n = 0; n < pl.size(); ++n)
102.     {
103.         if((n == i) || (n == j) || (n == k))
104.         {
105.             continue;
106.         }
107.
108.         if(distance(point, pl[n]) < radius)
109.         {
110.             return false;
111.         }
112.     }
113.
114.     return true;
115. }
116.
117. vector midpoint(const vector& a, const vector& b)
118. {
119.     return {(a[0] + b[0])/2.0, (a[1] + b[1])/2.0};
120. }
121.
122. // return a vector perpendicular to the line joining a and b
123. vector perpendicular(const vector& a, const vector& b)
124. {
125.     return {a[1] - b[1], b[0] - a[0]};
126. }
127.
128. // dot product
129. double dot(const vector& a, const vector& b)
130. {
131.     return (a[0]*b[0]) + (a[1]*b[1]);
132. }
133.
134. // Note: all maximum distance functions return -1.0
135. // to indicate that an invalid or irrelevant distance
136. // was the result. Since a maximum distance is being sought,
137. // this value will never be used.
138.
139. // local maximum distance from the line y = 0 (ye0)
140. double max_distance_to_ye0(size_t i, size_t j, const point_list& pl)
141. {
142.     auto mid = midpoint(pl[i], pl[j]);
143.     auto perp = perpendicular(pl[i], pl[j]);
144.
145.     // find point where mid + k*perp intersects y = 0
146.     auto t = -mid[1]/perp[1];
147.     auto x = mid[0] + t*perp[0];
148.     auto line_intersect = vector{x, 0.0};
149.     auto radius = distance(line_intersect, pl[i]);
150.     if(is_delaunay(pl, line_intersect, radius, i, j) &&
151.        inside_semicircle(line_intersect))
152.     {
154.     }
155.
156.     return -1.0;
157. }
158.
159. // Local maximum distance to the circle (r = 1)
160. double max_distance_to_re1(size_t i, size_t j, const point_list& pl)
161. {
162.     auto mid = midpoint(pl[i], pl[j]);
163.     auto perp = perpendicular(pl[i], pl[j]);
164.     auto dmp = dot(mid, perp);
165.     double dpp = dot(perp, perp);
166.     double dmm = dot(mid, mid);
167.     auto det = std::sqrt(dmp*dmp - dpp*(dmm-1.0));
168.     auto k1 = (-dmp + det)/dpp;
169.     auto k2 = (-dmp - det)/dpp;
170.     vector circle_intersect1 = {mid[0] + k1*perp[0], mid[1] + k1*perp[1]};
171.     vector circle_intersect2 = {mid[0] + k2*perp[0], mid[1] + k2*perp[1]};
172.     auto radius1 = distance(circle_intersect1, pl[i]);
173.     auto radius2 = distance(circle_intersect2, pl[i]);
174.     if(circle_intersect1[1] < 0.0)
175.     {
176.         radius1 = 3.0; // larger than any valid radius
177.     }
178.     if(circle_intersect2[1] < 0.0)
179.     {
180.         radius2 = 3.0;
181.     } // these two conditions cannot be true at the same time
182.
183.     auto circle_intersect = (radius1 < radius2 ? circle_intersect1 : circle_intersect2);
185.
186.     if( ! is_delaunay(pl, circle_intersect, radius, i, j))
187.     {
188.         return -1.0;
189.     }
190.
192. }
193.
194. // Find point equidistant from three points
195. vector circumcenter(size_t i, size_t j, size_t k, const point_list& pl)
196. {
197.     auto m0 = midpoint(pl[i], pl[j]);
198.     auto m1 = midpoint(pl[j], pl[k]);
199.
200.     auto p0 = perpendicular(pl[i], pl[j]);
201.     auto p1 = perpendicular(pl[j], pl[k]);
202.
203.     auto t = ((m1[1] - m0[1])*p0[0] + (m0[0] - m1[0])*p0[1])/(p1[0]*p0[1] - p1[1]*p0[0]);
204.
205.     return {m1[0] + t*p1[0], m1[1] + t*p1[1]};
206. }
207.
208. // return the maximum distance from three points while still occupying
209. // the triangle formed by the three points
210. double circumcircle_radius(size_t i, size_t j, size_t k, const point_list& pl)
211. {
212.     auto circumc = circumcenter(i, j, k, pl);
213.     auto radius = distance(pl[i], circumc);
214.
215.     assert(distance(pl[i], circumc) - distance(pl[j], circumc) < 1e-4);
216.     assert(distance(pl[k], circumc) - distance(pl[j], circumc) < 1e-4);
217.
219.        is_delaunay(pl, circumc, radius, i, j, k) &&
220.        inside_semicircle(circumc))
221.     {
223.     }
224.
225.     return -1.0;
226. }
227.
228. // find the maximum distance to any point in the point list
229. // using the above maximum distance functions.
230. double maximum_minimum_distance(const point_list& pl)
231. {
232.     auto maximum_distance = 0.0;
233.
234.     // Find maximum distance inside convex hull of points
235.     for(size_t i = 0; i < pl.size(); ++i)
236.     {
237.         for(size_t j = i + 1; j < pl.size(); ++j)
238.         {
239.             for(size_t k = j + 1; k < pl.size(); ++k)
240.             {
241.                 auto del_radius = circumcircle_radius(i, j, k, pl);
242.
243.                 if(del_radius > maximum_distance)
244.                 {
245.                     maximum_distance = del_radius;
246.                 }
247.             }
248.         }
249.     }
250.
251.     // For each pair of points, find the maximum distance to the
252.     // semicircular border and straight edge
253.     for(size_t i = 0; i < pl.size(); ++i)
254.     {
255.         for(size_t j = i + 1; j < pl.size(); ++j)
256.         {
257.             // find point where mid + k*perp intersects y = 0
258.             auto max_distance_to_line = max_distance_to_ye0(i, j, pl);
259.             if(max_distance_to_line > maximum_distance)
260.             {
261.                 maximum_distance = max_distance_to_line;
262.             }
263.
264.             // find point where mid + k*perp intersects x^2 + y^2 = 1
265.             auto max_distance_to_circle = max_distance_to_re1(i, j, pl);
266.             if(max_distance_to_circle > maximum_distance)
267.             {
268.                 maximum_distance = max_distance_to_circle;
269.             }
270.         }
271.     }
272.
273.     // check each point for distance to corners (-1, 0) and (1, 0)
274.     auto minimum_distance_to_left_corner = 2.0;
275.     for(const auto& p : pl)
276.     {
277.         auto d = distance(p, {-1.0, 0.0});
278.         if(d < minimum_distance_to_left_corner)
279.         {
280.             minimum_distance_to_left_corner = d;
281.         }
282.     }
283.     if(minimum_distance_to_left_corner > maximum_distance)
284.     {
285.         maximum_distance = minimum_distance_to_left_corner;
286.     }
287.
288.     auto minimum_distance_to_right_corner = 2.0;
289.     for(const auto& p : pl)
290.     {
291.         auto d = distance(p, {1.0, 0.0});
292.         if(d < minimum_distance_to_right_corner)
293.         {
294.             minimum_distance_to_right_corner = d;
295.         }
296.     }
297.     if(minimum_distance_to_right_corner > maximum_distance)
298.     {
299.         maximum_distance = minimum_distance_to_right_corner;
300.     }
301.
302.     return maximum_distance;
303. }
304.
305. int main()
306. {
307.     point_list pl;
308.     for(size_t i = 0; i < pl.size(); ++i)
309.     {
310.         pl[i] = {0.0, 0.5};
311.     }
312.     auto jostle_distance = 2.0;
313.     auto min_max_min_dis = 3.0;
314.     auto fail_count = 0;
315.     auto fail_threshold = 1000000;
316.     auto success = false;
317.
318.     while(jostle_distance < 4.0)
319.     {
320.         auto new_pl = pl;
321.         jostle_points_free(new_pl, jostle_distance);
322.         double max_min_dis = maximum_minimum_distance(new_pl);
323.         if(max_min_dis < min_max_min_dis)
324.         {
325.             pl = new_pl;
326.             min_max_min_dis = max_min_dis;
327.             fail_count = 0;
328.             success = true;
329.
330.             std::cout << "r = " << std::setprecision(std::numeric_limits<double>::max_digits10) << max_min_dis
331.                       << " % Jostle = " << jostle_distance << std::endl;
332.             ofs       << "r = " << std::setprecision(std::numeric_limits<double>::max_digits10) << max_min_dis
333.                       << " % Jostle = " << jostle_distance << std::endl;
334.             write_list(pl);
335.         }
336.         else
337.         {
338.             if(++fail_count >= fail_threshold)
339.             {
340.                 if(success)
341.                 {
342.                     // If new low found, take smaller steps
343.                     jostle_distance *= 0.5;
344.                 }
345.                 else
346.                 {
347.                     // if no new low found, take larger steps to escape local minimum
348.                     jostle_distance *= 2.0;
349.                 }
350.
351.                 // don't let the jostle_distance get so small that adding
352.                 // it to a position doesn't change anything
353.                 while(jostle_distance < std::numeric_limits<double>::epsilon())
354.                 {
355.                     jostle_distance *= 20.0;
356.                 }
357.
358.                 fail_count = 0;
359.                 success = false;
360.                 std::cout << "New jostle distance = " << jostle_distance << std::endl;
361.             }
362.         }
363.     }
364.
365.     return 0;
366. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top