AOJ 2051 - Rotation Estimation

http://judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=2051

概要

n 個の星の位置を表す写真が2種類与えられるので，最小でどれくらい回転させればその2つの写真が一致するかどうかを答える．

制約

• n <= 1000
• |座標| <= 100
• 2つの写真が一致するような角度は必ず存在する

解法

回転させても重心の位置は変わらないので，まず2つの写真を重心が原点になるように平行移動させる．

あとはある2点を決めて，それが一致するように片方を回転させたとき，他の星もすべて一致するかどうかを調べる．

 1 #include <cstdio>
 2 #include <vector>
 3 #include <complex>
 4 #include <cmath>
 5 #include <algorithm>
 6 using namespace std;
 7 typedef complex<double> P;
 8 static const double EPS = 1e-6;
 9 static const double PI = acos(-1.0);
10 namespace std {
11   bool operator<(const P& a, const P& b)
12   {
13     if (fabs(a.real() - b.real()) < EPS) {
14       return a.imag() < b.imag() - EPS;
15     } else {
16       return a.real() < b.real() - EPS;
17     }
18   }
19 };
20 inline double dot(const P& a, const P& b) { return a.real()*b.real() + a.imag()*b.imag(); }
21 
22 void move_to_origin(vector<P>& ps)
23 {
24   P p(0, 0);
25   for (vector<P>::const_iterator it = ps.begin(); it != ps.end(); ++it) {
26     p += *it;
27   }
28   p /= ps.size();
29   for (vector<P>::iterator it = ps.begin(); it != ps.end(); ++it) {
30     *it -= p;
31   }
32 }
33 
34 vector<P> rot(const vector<P>& ps, double theta)
35 {
36   const P p = polar(1.0, theta);
37   vector<P> ans;
38   for (vector<P>::const_iterator it = ps.begin(); it != ps.end(); ++it) {
39     ans.push_back(*it * p);
40   }
41   sort(ans.begin(), ans.end());
42   return ans;
43 }
44 
45 bool equal(const vector<P>& ps, const vector<P>& qs)
46 {
47   const int N = ps.size();
48   for (int i = 0; i < N; i++) {
49     if (abs(ps[i] - qs[i]) > EPS) {
50       return false;
51     }
52   }
53   return true;
54 }
55 
56 int main()
57 {
58   int N;
59   while (scanf("%d", &N) != EOF && N != 0) {
60     vector<P> ps, qs;
61     for (int i = 0; i < N; i++) {
62       double x, y;
63       scanf("%lf %lf", &x, &y);
64       ps.push_back(P(x, y));
65     }
66     for (int i = 0; i < N; i++) {
67       double x, y;
68       scanf("%lf %lf", &x, &y);
69       qs.push_back(P(x, y));
70     }
71 
72     move_to_origin(ps);
73     move_to_origin(qs);
74     sort(ps.begin(), ps.end());
75     const double arg0 = fmod(arg(ps[0]) + 2*PI, 2*PI);
76     double ans = 100;
77     for (int i = 0; i < N; i++) {
78       const double theta = arg0 - fmod(arg(qs[i]) + 2*PI, 2*PI);
79       const vector<P> rs = rot(qs, theta);
80       if (equal(ps, rs)) {
81         const double t = fabs(theta);
82         ans = min(ans, min(t, 2*PI-t));
83       }
84     }
85     printf("%.7f\n", ans);
86   }
87   return 0;
88 }