• API
• FAQ
• Tools
• Archive
daily pastebin goal
30%
SHARE
TWEET

# Sample code for closest pair of points

a guest Dec 6th, 2014 208 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.

Top