[计蒜客 商汤科技的行人检测（困难）]概率+几何

分类： Math probability

1. 题目链接

[计蒜客 商汤科技的行人检测（困难）]

2. 题意描述

商汤科技近日推出的 SenseVideo 能够对视频监控中的对象进行识别与分析，包括行人检测等。在行人检测问题中，最重要的就是对行人移动的检测。由于往往是在视频监控数据中检测行人，我们将图像上的行人抽象为二维平面上若干个的点。那么，行人的移动就相当于二维平面上的变换。

在这道题中，我们将行人的移动过程抽象为 旋转、伸缩、平移，有 4 个 移动参数：θ,scale,dx,dy
​​ 。每次行人的移动过程会将行人对应的 nn 个点全部依次应用旋转、伸缩、平移，对于平移前的点 (x, y)(x,y)，进行每种操作后的坐标如下：

旋转后的坐标为：(xcosθ−ysinθ,xsinθ+ycosθ)；
伸缩后的坐标为：(x×scale,y×scale)；
平移后的坐标为：(x+dx,y+dy)。
由于行人移动的特殊性，我们可以确保 0<scale≤10。和简单版本不同的是，这道题处理的坐标为浮点数而非整数。

很显然，通过变换前后的正确坐标，很容易算出行人的移动参数，但问题没有这么简单。由于行人实际的移动并不会完全按照我们预想的方式进行，因此，会有一部分变换后的坐标结果不正确，但可以确保 结果不正确的坐标数量严格不超过一半。

你现在作为商汤科技的实习生，接手了这个有趣的挑战：算出行人的移动参数。如果不存在一组合法的移动参数，则随意输出一组参数；如果有多种合法的移动参数，输出其中任意一组合法的即可。

输入格式

第一行输入一个整数 n，表示行人抽象出的点数。

接下来 n行，每行 4个 浮点数。前两个数表示平移前的坐标，后两个数表示平移后的坐标。

坐标范围在 −109​​ 到 109​ 之间，输入的坐标都保留到 6位小数。

对于中等版本，1≤n≤500；

对于困难版本，1≤n≤105。

输出格式

第一行输出一个浮点数θ，第二行输出一个浮点数 scale，第三行输出两个浮点数 dx,dy。

建议输出保留到 10 位小数或以上。我们会按照 10−3​的精度判断是否有超过一半的点变换后的坐标重合。

3. 解题思路

中等版本

在中等版本中，除了平移之外，还加入了旋转和拉伸。可以发现只需要枚举哪两对点是正确变换的，就可以计算出对应的拉伸、旋转、平移的量，从而验证是否有严格超过一半的点对满足这组变换。时间复杂度O(n3)。注意特判只有一个点对的情况。

困难版本

在困难版本中，点对数 n 从 100升级到了 100000。其算法本质并没有发生改变，依然是枚举两对点，然后验证。但是其枚举顺序，必须从按顺序枚举，改为随机枚举，以避免最坏复杂度。注意到错误的点对数严格不超过一半，因此我们有超过 1/4 的概率，枚举到的两对点就是正确的。对应的，枚举一次失败的概率就不足 3/4。这意味着：随机枚举 10 次，失败的概率不足 5.6%；随机枚举 20次，失败的概率不足 0.3%；随机枚举 50 次，失败的概率不足 0.00005%。所以，只需要常数次枚举，基本可以保证找到答案。时间复杂度O(n)。

注意：
这题精度问题比较严重，尽量不要用atan来求θ，精度误差会比较大。

4. 实现代码

#include <set>#include <map>#include <queue>#include <stack>#include <ctime>#include <cmath>#include <cctype>#include <cstdio>#include <string>#include <cstring>#include <cassert>#include <iomanip>#include <iostream>#include <algorithm>using namespace std;typedef long long LL;typedef long double LB;typedef unsigned long long ULL;typedef pair<int, int> PII;typedef pair<LL, LL> PLL;typedef pair<LB, LB> PLB;typedef vector<int> VI;const int INF = 0x3f3f3f3f;const LL INFL = 0x3f3f3f3f3f3f3f3fLL;const long double PI = acos(-1.0);const long double eps = 1e-4;void debug() { cout << endl; }template<typename T, typename ...R> void debug (T f, R ...r) { cout << "[" << f << "]"; debug (r...); }template<typename T> inline void umax(T &a, T b) { a = max(a, b); }template<typename T> inline void umin(T &a, T b) { a = min(a, b); }template <typename T> inline bool scan_d (T &ret) {    char c; int sgn;    if (c = getchar(), c == EOF) return 0; //EOF    while (c != '-' && (c < '0' || c > '9') ) if((c = getchar()) == EOF) return 0;    sgn = (c == '-') ? -1 : 1;    ret = (c == '-') ? 0 : (c - '0');    while (c = getchar(), c >= '0' && c <= '9') ret = ret * 10 + (c - '0');    ret *= sgn;    return 1;}template<typename T> void print(T x) {    static char s[33], *s1; s1 = s;    if (!x) *s1++ = '0';    if (x < 0) putchar('-'), x = -x;    while(x) *s1++ = (x % 10 + '0'), x /= 10;    while(s1-- != s) putchar(*s1);}template<typename T> void println(T x) { print(x); putchar('\n'); }const int MAXN = 1e5 + 5;int n;struct QNode {    LB x[2], y[2];    QNode() {}} cmd[MAXN];int sgn(const LB& a, const LB& b) {    if(fabs(b - a) <= eps) return 0;    return a < b ? -1 : 1;}PLB f(QNode qd, LB sita, LB scale, LB dx, LB dy) {    LB a = (qd.x[0] * cos(sita) - qd.y[0] * sin(sita)) * scale + dx;    LB b = (qd.x[0] * sin(sita) + qd.y[0] * cos(sita)) * scale + dy;    return make_pair(a, b);}int calc(LB sita, LB scale, LB dx, LB dy) {    int cnt = 0;    for(int i = 1; i <= n; ++i) {        PLB rs = f(cmd[i], sita, scale, dx, dy);        if(!sgn(rs.first, cmd[i].x[1]) && !sgn(rs.second, cmd[i].y[1])) ++ cnt;    }    return cnt;}inline LB dis(LB ax, LB ay, LB bx = 0.0, LB by = 0.0) {    LB a = bx - ax;    LB b = by - ay;    return sqrt(a * a + b * b);}int main() {#ifdef ___LOCAL_WONZY___    freopen ("input.txt", "r", stdin);#endif // ___LOCAL_WONZY___    while(~scanf("%d", &n)) {        double x[2], y[2];        for(int i = 1; i <= n; ++i) {            scanf("%lf %lf %lf %lf", &x[0], &y[0], &x[1], &y[1]);            cmd[i].x[0] = x[0]; cmd[i].x[1] = x[1];            cmd[i].y[0] = y[0]; cmd[i].y[1] = y[1];        }        LB sita, scale, dx, dy, delta_x[2], delta_y[2], up, dw;        LB ans[4];        if(n == 1) {            ans[0] = 0.0;            ans[1] = 1.0;            ans[2] = cmd[1].x[1] - cmd[1].x[0];            ans[3] = cmd[1].y[1] - cmd[1].y[0];            cout << fixed << setprecision(11) << ans[0] << endl;            cout << fixed << setprecision(11) << ans[1] << endl;            cout << fixed << setprecision(11) << ans[2] << " " << ans[3] << endl;            continue;        }        int cnt = 0;        for(int _ = 100; _--; ) {            int i = rand() % n + 1;            int j = rand() % n + 1;            if(i == j) continue;            up = dis(cmd[i].x[1], cmd[i].y[1], cmd[j].x[1], cmd[j].y[1]);            dw = dis(cmd[i].x[0], cmd[i].y[0], cmd[j].x[0], cmd[j].y[0]);            if(!sgn(dw, 0.0)) continue;            scale = up / dw;            if(sgn(scale, 10.0) == 1 || sgn(scale, 0) == -1) continue;            delta_x[0] = cmd[j].x[0] - cmd[i].x[0];            delta_y[0] = cmd[j].y[0] - cmd[i].y[0];            delta_x[1] = cmd[j].x[1] - cmd[i].x[1];            delta_y[1] = cmd[j].y[1] - cmd[i].y[1];            up = delta_x[0] * delta_x[1] + delta_y[0] * delta_y[1];            dw = dis(delta_x[0], delta_y[0]) * dis(delta_x[1], delta_y[1]);            sita = acos(up / dw);            dx = cmd[i].x[1] - (cmd[i].x[0] * cos(sita) - cmd[i].y[0] * sin(sita)) * scale;            dy = cmd[i].y[1] - (cmd[i].x[0] * sin(sita) + cmd[i].y[0] * cos(sita)) * scale;            int ret = calc(sita, scale, dx, dy);            if(ret <= cnt) continue;            ans[0] = sita;            ans[1] = scale;            ans[2] = dx;            ans[3] = dy;            cnt = ret;            if(cnt * 2 >= n) break;        }        cout << fixed << setprecision(11) << ans[0] << endl;        cout << fixed << setprecision(11) << ans[1] << endl;        cout << fixed << setprecision(11) << ans[2] << " " << ans[3] << endl;    }#ifdef ___LOCAL_WONZY___    cout << "Time elapsed: " << 1.0 * clock() / CLOCKS_PER_SEC * 1000 << " ms." << endl;#endif // ___LOCAL_WONZY___    return 0;}