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

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

分类：

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;
}