daily pastebin goal
4%
SHARE
TWEET

Sample code for closest pair of points

a guest Dec 6th, 2014 206 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <climits>
  2. #include <cmath>
  3. #include <algorithm>
  4. #include <chrono>
  5. #include <iostream>
  6. #include <numeric>
  7. #include <random>
  8. #include <utility>
  9.  
  10. template<typename CT>
  11. struct Point {
  12.     typedef CT CoordinateType;
  13.  
  14.     CoordinateType x;
  15.     CoordinateType y;
  16. };
  17.  
  18. template<typename CT>
  19. std::ostream &operator<<(std::ostream &os, const Point<CT> &p)
  20. {
  21.     return os << "x:" << p.x << " y:" << p.y;
  22. }
  23.  
  24. template<typename T1, typename T2>
  25. std::ostream &operator<<(std::ostream &os, const std::pair<T1, T2> &pair)
  26. {
  27.     return os << " (" << pair.first << "," << pair.second << ") ";
  28. }
  29.  
  30. template<typename CT>
  31. CT squareRoot(CT n) {
  32.     return static_cast<CT>(std::pow(n, 0.5));
  33. }
  34.  
  35. template<typename CT>
  36. CT distanceSquared(const Point<CT> &first, const Point<CT> &second) {
  37.     const CT dx = (first.x - second.x), dy = (first.y - second.y);
  38.     return (dx*dx) + (dy*dy);
  39. }
  40.  
  41. template<typename CT>
  42. CT distance(const Point<CT> &first, const Point<CT> &second) {
  43.     return static_cast<CT>(std::pow(distanceSquared(first, second), 0.5));
  44. }
  45.  
  46. template<typename PointType, typename PairType = std::pair<PointType, PointType>, typename CT = typename PointType::CoordinateType >
  47. std::pair<PairType, CT> exhaustiveSearch(const PointType p[], int begin, int end)
  48. {
  49.     PairType pair;
  50.     CT minSquared = std::numeric_limits<CT>::max();
  51.  
  52.     for (int i = begin; i < (end - 1); ++i) {
  53.         const PointType &first(p[i]);
  54.  
  55.         for (int j = i + 1; j < end; ++j) {
  56.             const PointType &second(p[j]);
  57.             const CT distSquared = distanceSquared(second, first);
  58.             if (distSquared < minSquared) {
  59.                 minSquared = distSquared;
  60.                 pair = std::make_pair(first, second);
  61.             }
  62.         }
  63.     }
  64.  
  65.     return std::make_pair(pair, squareRoot(minSquared));
  66. }
  67.  
  68. template<typename PointType, typename PairType = std::pair<PointType, PointType>, typename CT = typename PointType::CoordinateType >
  69. std::pair<PairType, CT> recursiveSearch(const PointType p[], int begin, int end)
  70. {
  71.     const int span = end - begin;
  72.     if (span < 64) {
  73.         return exhaustiveSearch(p, begin, end);
  74.     }
  75.  
  76.     const int mid = begin + (span + 1)/2;
  77.  
  78.     std::pair<PairType, CT> lresult = recursiveSearch(p, begin, mid);
  79.     std::pair<PairType, CT> rresult = recursiveSearch(p, mid, end);
  80.  
  81.     PairType &lpair = lresult.first;
  82.     CT &lmin = lresult.second;
  83.  
  84.     PairType &rpair = rresult.first;
  85.     CT &rmin = rresult.second;
  86.    
  87.     {
  88.         CT min = std::min(lmin, rmin);
  89.         CT minSquared = min * min;
  90.         CT upperBound = p[mid - 1].x + min;
  91.  
  92.         for (int i = mid; i < end; ++i) {
  93.             const PointType &rp = p[i];
  94.             if (rp.x >= upperBound) {
  95.                 break;
  96.             } else {
  97.                 CT lowerBound = rp.x - min;
  98.                 for (int j = mid - 1; j >= begin; --j) {
  99.                     const PointType &lp = p[j];
  100.                     if (lp.x <= lowerBound) {
  101.                         break;
  102.                     } else {
  103.                         const CT distSquared = distanceSquared(lp, rp);
  104.                         if (distSquared < minSquared) {
  105.                             minSquared = distSquared;
  106.                             min = squareRoot(distSquared);
  107.                             lmin = min;
  108.                             lpair = std::make_pair(lp, rp);
  109.  
  110.                             upperBound = p[mid - 1].x + min;
  111.                             lowerBound = rp.x - min;
  112.                         }
  113.                     }
  114.                 }
  115.             }
  116.         }
  117.     }
  118.     {
  119.         CT min = std::min(lmin, rmin);
  120.         CT minSquared = min * min;
  121.         CT lowerBound = p[mid].x - min;
  122.  
  123.         for (int i = mid - 1; i >= begin; --i) {
  124.             const PointType &lp = p[i];
  125.             if (lp.x <= lowerBound) {
  126.                 break;
  127.             } else {
  128.                 CT upperBound = lp.x + min;
  129.                 for (int j = mid; j < end; ++j) {
  130.                     const PointType &rp = p[j];
  131.                     if (rp.x >= upperBound) {
  132.                         break;
  133.                     } else {
  134.                         const CT distSquared = distanceSquared(lp, rp);
  135.                         if (distSquared < minSquared) {
  136.                             minSquared = distSquared;
  137.                             min = squareRoot(distSquared);
  138.                             rmin = min;
  139.                             rpair = std::make_pair(lp, rp);
  140.  
  141.                             lowerBound = p[mid].x - min;
  142.                             upperBound = lp.x + min;
  143.                         }
  144.                     }
  145.                 }
  146.             }
  147.         }
  148.     }
  149.  
  150.     return std::make_pair((rmin < lmin ? rpair : lpair), (rmin < lmin ? rmin : lmin));
  151. }
  152.  
  153. template<typename PointType, typename PairType = std::pair<PointType, PointType>, typename CT = typename PointType::CoordinateType >
  154. std::pair<PairType, CT> neighbourSearch(const PointType p[], int begin, int end)
  155. {
  156.     PairType pair;
  157.     CT minSquared = distanceSquared(p[begin], p[begin + 1]);
  158.     CT min = squareRoot(minSquared);
  159.  
  160.     for (int i = begin; i < (end - 1); ++i) {
  161.         const PointType &first = p[i];
  162.         CT upperBound = first.x + min;
  163.  
  164.         for (int j = i + 1; j < end; ++j) {
  165.             const PointType &second = p[j];
  166.             if (second.x >= upperBound) {
  167.                 break;
  168.             } else {
  169.                 const CT distSquared = distanceSquared(first, second);
  170.                 if (distSquared < minSquared) {
  171.                     minSquared = distSquared;
  172.                     min = squareRoot(distSquared);
  173.                     pair = std::make_pair(first, second);
  174.  
  175.                     upperBound = first.x + min;
  176.                 }
  177.             }
  178.         }
  179.     }
  180.  
  181.     return std::make_pair(pair, min);
  182. }
  183.  
  184. template<typename CT>
  185. void timeOperations()
  186. {
  187. #ifdef USE_EXHAUSTIVE
  188.     const int N = 40 * 1000;
  189. #else
  190.     const int N = 400 * 1000;
  191. #endif
  192.  
  193.     std::vector<Point<CT> > points(N);
  194.     std::vector<Point<CT> > sorted(N);
  195.  
  196.     std::default_random_engine rng(std::chrono::system_clock::to_time_t(std::chrono::system_clock::now()));
  197.     std::uniform_real_distribution<CT> dst(0.0, 1.0);
  198.  
  199.     for (int i = 0; i < N; ++i) { points[i] = Point<CT>{ dst(rng), dst(rng) }; }
  200.  
  201.     std::pair<std::pair<Point<CT>, Point<CT> >, CT> min;
  202.  
  203.     std::chrono::system_clock::time_point tp = std::chrono::system_clock::now();
  204.     std::chrono::system_clock::duration elapsed;
  205.  
  206. #ifdef USE_EXHAUSTIVE
  207.     min = exhaustiveSearch(points.data(), 0, N);
  208.     elapsed = std::chrono::system_clock::now() - tp;
  209.     std::cout << "exhaustive min:" << min.second << " distance:" << min.first << " elapsed:" << elapsed.count() / (1000 * 1000) << std::endl;
  210.  
  211.     tp = std::chrono::system_clock::now();
  212. #endif
  213.  
  214.     for (int i = 0; i < N; ++i) { sorted[i] = points[i]; }
  215.     std::sort(std::begin(sorted), std::end(sorted), [](const Point<CT> &l, const Point<CT> &r) { return l.x < r.x; });
  216.     min = recursiveSearch(sorted.data(), 0, N);
  217.     elapsed = std::chrono::system_clock::now() - tp;
  218.     std::cout << "recursive min:" << min.second << " distance:" << min.first << " elapsed:" << elapsed.count() / (1000 * 1000) << std::endl;
  219.  
  220.     tp = std::chrono::system_clock::now();
  221.  
  222.     for (int i = 0; i < N; ++i) { sorted[i] = points[i]; }
  223.     std::sort(std::begin(sorted), std::end(sorted), [](const Point<CT> &l, const Point<CT> &r) { return l.x < r.x; });
  224.     min = neighbourSearch(sorted.data(), 0, N);
  225.     elapsed = std::chrono::system_clock::now() - tp;
  226.     std::cout << "neighbour min:" << min.second << " distance:" << min.first << " elapsed:" << elapsed.count() / (1000 * 1000) << std::endl;
  227. }
  228.  
  229. int main(int,char**)
  230. {
  231.     timeOperations<float>();
  232.     return 0;
  233. }
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