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,
  93.                  double radius,
  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.     {
  153.         return radius;
  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);
  184.     auto radius = (radius1 < radius2 ? radius1 : radius2);
  185.  
  186.     if( ! is_delaunay(pl, circle_intersect, radius, i, j))
  187.     {
  188.         return -1.0;
  189.     }
  190.  
  191.     return radius;
  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.  
  218.     if(std::isfinite(radius) &&
  219.        is_delaunay(pl, circumc, radius, i, j, k) &&
  220.        inside_semicircle(circumc))
  221.     {
  222.         return radius;
  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. OK, I Understand
 
Top