• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# semicircle points

a guest Sep 23rd, 2015 132 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. // randomly moves points in point_list while obeying symmetry constraints
34. void jostle_points_constrained(point_list& pl, double jostle_distance)
35. {
36.     size_t index;
37.     if(pl.size() == 6)
38.     {
39.         index = 2;
40.     }
41.     else
42.     {
43.         assert(false);
44.     }
45.
46.     // First point(s) stay on the y-axis
47.     for(auto i = 0; i < index; ++i)
48.     {
49.         while(true)
50.         {
51.             auto distance = random_uniform(-jostle_distance, jostle_distance);
52.             auto np = pl[i];
53.             np[0] = 0.0;
54.             np[1] += distance;
55.             if(inside_semicircle(np))
56.             {
57.                 pl[i] = np;
58.                 break;
59.             }
60.         }
61.     }
62.
63.     // Other points are mirrored across the y-axis
64.     for(size_t i = index; i < pl.size(); i += 2)
65.     {
66.         while(true)
67.         {
68.             auto distance = random_uniform(0, jostle_distance);
69.             auto angle = random_uniform(0, 2*pi);
70.
71.             auto dx = distance*std::cos(angle);
72.             auto dy = distance*std::sin(angle);
73.
74.             auto np = pl[i];
75.             np[0] += dx;
76.             np[1] += dy;
77.             if(inside_semicircle(np))
78.             {
79.                 pl[i] = np;
80.                 pl[i+1] = np;
81.                 pl[i+1][0] *= -1.0; // reflect in y-axis
82.                 break;
83.             }
84.         }
85.     }
86. }
87.
88. // raondomly move a single point
89. void jostle_point(vector& p, double jostle_distance)
90. {
91.     while(true)
92.     {
93.         auto distance = random_uniform(0, jostle_distance);
94.         auto angle = random_uniform(0, 2*pi);
95.
96.         auto dx = distance*std::cos(angle);
97.         auto dy = distance*std::sin(angle);
98.
99.         auto np = p;
100.         np[0] += dx;
101.         np[1] += dy;
102.         if(inside_semicircle(np))
103.         {
104.             p = np;
105.             break;
106.         }
107.     }
108. }
109.
110. // Move all points randomly and independantly
111. void jostle_points_free(point_list& pl, double jostle_distance)
112. {
113.     for(auto& p : pl)
114.     {
115.         jostle_point(p, jostle_distance);
116.     }
117. }
118.
119. // Write out point list to output stream
120. void write_list(const point_list& pl, std::ostream& os)
121. {
122.     os << "p = [\n";
123.     for(const auto& p : pl)
124.     {
125.         for(auto e : p)
126.         {
127.             os << std::setprecision(std::numeric_limits<double>::max_digits10) << e << ' ';
128.         }
129.         os << ";\n";
130.     }
131.     os << "];\n=====================================================" << std::endl;
132. }
133.
134. // write point list to console and file
135. std::ofstream ofs("result.txt");
136. void write_list(const point_list& pl)
137. {
138.     write_list(pl, std::cout);
139.     write_list(pl, ofs);
140. }
141.
142. // Given a radius associated with points i, j, k
143. // in the point list, check if any othe points
144. // fall within a radius of the given point (usually
145. // the circumcenter of i, j, k).
146. bool is_delaunay(const point_list& pl,
147.                  const vector& point,
149.                  size_t i, size_t j, size_t k = -1)
150. {
152.     {
153.         return false;
154.     }
155.
156.     for(size_t n = 0; n < pl.size(); ++n)
157.     {
158.         if((n == i) || (n == j) || (n == k))
159.         {
160.             continue;
161.         }
162.
164.         {
165.             return false;
166.         }
167.     }
168.
169.     return true;
170. }
171.
172. vector midpoint(const vector& a, const vector& b)
173. {
174.     return {(a[0] + b[0])/2.0, (a[1] + b[1])/2.0};
175. }
176.
177. // return a vector perpendicular to the line joining a and b
178. vector perpendicular(const vector& a, const vector& b)
179. {
180.     return {a[1] - b[1], b[0] - a[0]};
181. }
182.
183. // dot product
184. double dot(const vector& a, const vector& b)
185. {
186.     return (a[0]*b[0]) + (a[1]*b[1]);
187. }
188.
189. // Note: all maximum distance functions return -1.0
190. // to indicate that an invalid or irrelevant distance
191. // was the result. Since a maximum distance is being sought,
192. // this value will never be used.
193.
194. // local maximum distance from the line y = 0 (ye0)
195. double max_distance_to_ye0(size_t i, size_t j, const point_list& pl)
196. {
197.     auto mid = midpoint(pl[i], pl[j]);
198.     auto perp = perpendicular(pl[i], pl[j]);
199.
200.     // find point where mid + k*perp intersects y = 0
201.     auto t = -mid[1]/perp[1];
202.     auto x = mid[0] + t*perp[0];
203.     auto line_intersect = vector{x, 0.0};
204.     auto radius = distance(line_intersect, pl[i]);
205.     if(is_delaunay(pl, line_intersect, radius, i, j) &&
206.        inside_semicircle(line_intersect))
207.     {
209.     }
210.
211.     return -1.0;
212. }
213.
214. // Local maximum distance to the circle (r = 1)
215. double max_distance_to_re1(size_t i, size_t j, const point_list& pl)
216. {
217.     auto mid = midpoint(pl[i], pl[j]);
218.     auto perp = perpendicular(pl[i], pl[j]);
219.     auto dmp = dot(mid, perp);
220.     double dpp = dot(perp, perp);
221.     double dmm = dot(mid, mid);
222.     auto det = std::sqrt(dmp*dmp - dpp*(dmm-1.0));
223.     auto k1 = (-dmp + det)/dpp;
224.     auto k2 = (-dmp - det)/dpp;
225.     vector circle_intersect1 = {mid[0] + k1*perp[0], mid[1] + k1*perp[1]};
226.     vector circle_intersect2 = {mid[0] + k2*perp[0], mid[1] + k2*perp[1]};
227.     auto radius1 = distance(circle_intersect1, pl[i]);
228.     auto radius2 = distance(circle_intersect2, pl[i]);
229.     if(circle_intersect1[1] < 0.0)
230.     {
232.     }
233.     if(circle_intersect2[1] < 0.0)
234.     {
236.     } // these two conditions cannot be true at the same time
237.
240.
241.     if( ! is_delaunay(pl, circle_intersect, radius, i, j))
242.     {
243.         return -1.0;
244.     }
245.
247. }
248.
249. // Find point equidistant from three points
250. vector circumcenter(size_t i, size_t j, size_t k, const point_list& pl)
251. {
252.     auto m0 = midpoint(pl[i], pl[j]);
253.     auto m1 = midpoint(pl[j], pl[k]);
254.
255.     auto p0 = perpendicular(pl[i], pl[j]);
256.     auto p1 = perpendicular(pl[j], pl[k]);
257.
258.     auto t = ((m1[1] - m0[1])*p0[0] + (m0[0] - m1[0])*p0[1])/(p1[0]*p0[1] - p1[1]*p0[0]);
259.
260.     return {m1[0] + t*p1[0], m1[1] + t*p1[1]};
261. }
262.
263. // return the maximum distance from three points while still occupying
264. // the triangle formed by the three points
265. double circumcircle_radius(size_t i, size_t j, size_t k, const point_list& pl)
266. {
267.     auto circumc = circumcenter(i, j, k, pl);
268.     auto radius = distance(pl[i], circumc);
269.
270.     assert(distance(pl[i], circumc) - distance(pl[j], circumc) < 1e-4);
271.     assert(distance(pl[k], circumc) - distance(pl[j], circumc) < 1e-4);
272.
274.        is_delaunay(pl, circumc, radius, i, j, k) &&
275.        inside_semicircle(circumc))
276.     {
278.     }
279.
280.     return -1.0;
281. }
282.
283. // find the maximum distance to any point in the point list
284. // using the above maximum distance functions.
285. double maximum_minimum_distance(const point_list& pl)
286. {
287.     auto maximum_distance = 0.0;
288.
289.     // Find maximum distance inside convex hull of points
290.     for(size_t i = 0; i < pl.size(); ++i)
291.     {
292.         for(size_t j = i + 1; j < pl.size(); ++j)
293.         {
294.             for(size_t k = j + 1; k < pl.size(); ++k)
295.             {
297.
299.                 {
301.                 }
302.             }
303.         }
304.     }
305.
306.     // For each pair of points, find the maximum distance to the
307.     // semicircular border and straight edge
308.     for(size_t i = 0; i < pl.size(); ++i)
309.     {
310.         for(size_t j = i + 1; j < pl.size(); ++j)
311.         {
312.             // find point where mid + k*perp intersects y = 0
313.             auto max_distance_to_line = max_distance_to_ye0(i, j, pl);
314.             if(max_distance_to_line > maximum_distance)
315.             {
316.                 maximum_distance = max_distance_to_line;
317.             }
318.
319.             // find point where mid + k*perp intersects x^2 + y^2 = 1
320.             auto max_distance_to_circle = max_distance_to_re1(i, j, pl);
321.             if(max_distance_to_circle > maximum_distance)
322.             {
323.                 maximum_distance = max_distance_to_circle;
324.             }
325.         }
326.     }
327.
328.     // check each point for distance to corners (-1, 0) and (1, 0)
329.     auto minimum_distance_to_left_corner = 2.0;
330.     for(const auto& p : pl)
331.     {
332.         auto d = distance(p, {-1.0, 0.0});
333.         if(d < minimum_distance_to_left_corner)
334.         {
335.             minimum_distance_to_left_corner = d;
336.         }
337.     }
338.     if(minimum_distance_to_left_corner > maximum_distance)
339.     {
340.         maximum_distance = minimum_distance_to_left_corner;
341.     }
342.
343.     auto minimum_distance_to_right_corner = 2.0;
344.     for(const auto& p : pl)
345.     {
346.         auto d = distance(p, {1.0, 0.0});
347.         if(d < minimum_distance_to_right_corner)
348.         {
349.             minimum_distance_to_right_corner = d;
350.         }
351.     }
352.     if(minimum_distance_to_right_corner > maximum_distance)
353.     {
354.         maximum_distance = minimum_distance_to_right_corner;
355.     }
356.
357.     return maximum_distance;
358. }
359.
360. int main()
361. {
362.     point_list pl;
363.     for(size_t i = 0; i < pl.size(); ++i)
364.     {
365.         pl[i] = {0.0, 0.5};
366.     }
367.     auto jostle_distance = 1e-4;
368.     auto min_max_min_dis = 3.0;
369.     auto fail_count = 0;
370.     auto fail_threshold = 1000000;
371.     auto success = false;
372.
373.     while(jostle_distance < 4.0)
374.     {
375.         auto new_pl = pl;
376.         jostle_points_constrained(new_pl, jostle_distance);
377.         double max_min_dis = maximum_minimum_distance(new_pl);
378.         if(max_min_dis < min_max_min_dis)
379.         {
380.             pl = new_pl;
381.             min_max_min_dis = max_min_dis;
382.             fail_count = 0;
383.             success = true;
384.
385.             std::cout << "r = " << std::setprecision(std::numeric_limits<double>::max_digits10) << max_min_dis
386.                       << " % Jostle = " << jostle_distance << std::endl;
387.             ofs       << "r = " << std::setprecision(std::numeric_limits<double>::max_digits10) << max_min_dis
388.                       << " % Jostle = " << jostle_distance << std::endl;
389.             write_list(pl);
390.         }
391.         else
392.         {
393.             if(++fail_count >= fail_threshold)
394.             {
395.                 if(success)
396.                 {
397.                     // If new low found, take smaller steps
398.                     jostle_distance *= 0.5;
399.                 }
400.                 else
401.                 {
402.                     // if no new low found, take larger steps to escape local minimum
403.                     jostle_distance *= 2.0;
404.                 }
405.
406.                 // don't let the jostle_distance get so small that adding
407.                 // it to a position doesn't change anything
408.                 while(jostle_distance < std::numeric_limits<double>::epsilon())
409.                 {
410.                     jostle_distance *= 20.0;
411.                 }
412.
413.                 fail_count = 0;
414.                 success = false;
415.                 std::cout << "New jostle distance = " << jostle_distance << std::endl;
416.             }
417.         }
418.     }
419.
420.     return 0;
421. }
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