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

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

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

• 旋转后的坐标为：$(x\cos \theta -y\sin \theta ,x\sin \theta +y\cos \theta )$；
• 伸缩后的坐标为：$(x\times scale,y\times scale)$；
• 平移后的坐标为：(x+dx,y+dy)(x + d_x, y + d_y)(x+d​x​​,y+d​y​​)。

由于行人移动的特殊性，我们可以确保 $0<scale\le 10$。和简单版本不同的是，这道题处理的坐标为浮点数而非整数。

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

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

输入格式

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

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

坐标范围在 −109-10^9−10​9​​ 到 10910^910​9​​ 之间，输入的坐标都保留到 $6$ 位小数。

对于中等版本，$1\le n\le 500$；

对于困难版本，1≤n≤1051 \le n \le 10^{5}1≤n≤10​5​​。

输出格式

第一行输出一个浮点数 $\theta$，第二行输出一个浮点数 $scale$，第三行输出两个浮点数 dx,dyd_x,d_yd​x​​,d​y​​。

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

样例输入

5
0 0 -1 1
0 1 -2 1
1 0 -1 2
1 1 0 0
2 1 1 0

样例输出

1.5707963268
1
-1 1

～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～～

计算几何+随机化～

每个行人可以列出两个方程，所以两个行人就能解出四个参数。因为有一半的行人是正确的，所以我们枚举的正确率为0.5^2=0.25，枚举20次正确率0.997，几近正确。

至于四个方程怎么解出参数，我们不需要直接解方程，而是利用计算几何知识：

角度：转化后两点成的直线与转化前两点成的直线的夹角即为答案，我们利用atan2(x,y)求解。

倍数：转化后两点的距离/转化前两点的距离。

x,y：直接减。

需要注意的是这题答案不唯一，刚开始检查了半天结果是对的……

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<cmath>
using namespace std;
#define eps 1e-4

int n,x,y,num;

struct node{
	double x,y;
}a[100001],b[100001];

struct node1{
	double ang,k,x,y;
}tmp;

bool operator < (node1 u,node1 v)
{
	u.ang<v.ang;
}

node operator - (node u,node v)
{
	return (node){u.x-v.x,u.y-v.y};
}

bool operator == (node u,node v)
{
	return fabs(u.x-v.x)<eps && fabs(u.y-v.y)<eps;
}

double angle(node u)
{
	return atan2(u.y,u.x);
}

double dis(node u,node v)
{
	return sqrt((u.x-v.x)*(u.x-v.x)+(u.y-v.y)*(u.y-v.y));
}

node cal(double u,node v)
{
	double x=v.x,y=v.y;
	v.x=x*cos(u)-y*sin(u);
	v.y=x*sin(u)+y*cos(u);
	return v;
}

bool chec(node1 u)
{
	int now=0;
	node kkk;
	for(int i=1;i<=n;i++)
	{
		kkk=cal(u.ang,a[i]);
		kkk.x=kkk.x*u.k+u.x;
		kkk.y=kkk.y*u.k+u.y;
		if(kkk==b[i]) continue;
		now++;
		if(now*2>n) return 0;
	}
	return 1;
}

bool chec(int u,int v)
{
	tmp.ang=angle(b[v]-b[u])-angle(a[v]-a[u]);
	tmp.k=dis(b[u],b[v])/dis(a[u],a[v]);
	node now=cal(tmp.ang,a[u]);
	tmp.x=b[u].x-tmp.k*now.x;
	tmp.y=b[u].y-tmp.k*now.y;
	return chec(tmp);
}

int main()
{
	scanf("%d",&n);
	srand(20000619);
	for(int i=1;i<=n;i++) scanf("%lf%lf%lf%lf",&a[i].x,&a[i].y,&b[i].x,&b[i].y);
	if(n==1)
	{
		printf("0\n1\n%.10lf %.10lf",b[1].x-a[1].x,b[1].y-a[1].y);
		return 0;
	}
	for(int kk=20;kk;kk--)
	{
		x=rand()%n+1;y=rand()%n+1;
		while(x==y) y=rand()%n+1;
		if(chec(x,y)) break;
	}
	printf("%.10lf\n%.10lf\n%.10lf %.10lf",tmp.ang,tmp.k,tmp.x,tmp.y);
	return 0;
}