geometry/closestPair.cpp


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

double squaredDist(pt a, pt b) {
	return real(conj(a-b) * (a-b));
}

bool compY(pt a, pt b) {
	return (imag(a) == imag(b)) ? real(a) < real(b) 
															: imag(a) < imag(b);
}

bool compX(pt a, pt b) {
	return (real(a) == real(b)) ? imag(a) < imag(b)
															: real(a) < real(b);
}

double shortestDist(vector<pt>& pts) { // pts.size() > 1
	set<pt, bool(*)(pt, pt)> status(compY);
	sort(all(pts), compX);
	double opt = 1.0/0.0, sqrtOpt = 1.0/0.0;
	auto left = pts.begin(), right = pts.begin();
	status.insert(*right); right++;

	while (right != pts.end()) {
		if (left != right &&
				abs(real(*left - *right)) >= sqrtOpt) {
			status.erase(*left);
			left++;
		} else {
			auto lower = status.lower_bound({-1.0/0.0, imag(*right) - sqrtOpt});
			auto upper = status.upper_bound({-1.0/0.0, imag(*right) + sqrtOpt});
			for (;lower != upper; lower++) {
				double cand = squaredDist(*right, *lower);
				if (cand < opt) {
					opt = cand;
					sqrtOpt = sqrt(opt);
			}}
			status.insert(*right);
			right++;
	}}
	return sqrtOpt;
}