#include<stdio.h>
#include<algorithm>
#include<string.h>
#include<math.h>
#include<stdlib.h>
using namespace std;
const int MAXN=550;
const double eps=1e-8;
struct Point
{
double x,y,z;
Point(){}
Point(double xx,double yy,double zz):x(xx),y(yy),z(zz){}
//两向量之差
Point operator -(const Point p1)
{
return Point(x-p1.x,y-p1.y,z-p1.z);
}
//两向量之和
Point operator +(const Point p1)
{
return Point(x+p1.x,y+p1.y,z+p1.z);
}
//叉乘
Point operator *(const Point p)
{
return Point(y*p.z-z*p.y,z*p.x-x*p.z,x*p.y-y*p.x);
}
Point operator *(double d)
{
return Point(x*d,y*d,z*d);
}
Point operator / (double d)
{
return Point(x/d,y/d,z/d);
}
//点乘
double operator ^(Point p)
{
return (x*p.x+y*p.y+z*p.z);
}
};
struct CH3D
{
struct face
{
//表示凸包一个面上的三个点的编号
int a,b,c;
//表示该面是否属于最终凸包上的面
bool ok;
};
//初始顶点数
int n;
//初始顶点
Point P[MAXN];
//凸包表面的三角形数
int num;
//凸包表面的三角形
face F[8*MAXN];
//凸包表面的三角形
int g[MAXN][MAXN];
//向量长度
double vlen(Point a)
{
return sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
}
//叉乘
Point cross(const Point &a,const Point &b,const Point &c)
{
return Point((b.y-a.y)*(c.z-a.z)-(b.z-a.z)*(c.y-a.y),
(b.z-a.z)*(c.x-a.x)-(b.x-a.x)*(c.z-a.z),
(b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x)
);
}
//三角形面积*2
double area(Point a,Point b,Point c)
{
return vlen((b-a)*(c-a));
}
//四面体有向体积*6
double volume(Point a,Point b,Point c,Point d)
{
return (b-a)*(c-a)^(d-a);
}
//正:点在面同向
double dblcmp(Point &p,face &f)
{
Point m=P[f.b]-P[f.a];
Point n=P[f.c]-P[f.a];
Point t=p-P[f.a];
return (m*n)^t;
}
void deal(int p,int a,int b)
{
int f=g[a][b];//搜索与该边相邻的另一个平面
face add;
if(F[f].ok)
{
if(dblcmp(P[p],F[f])>eps)
dfs(p,f);
else
{
add.a=b;
add.b=a;
add.c=p;//这里注意顺序,要成右手系
add.ok=true;
g[p][b]=g[a][p]=g[b][a]=num;
F[num++]=add;
}
}
}
void dfs(int p,int now)//递归搜索所有应该从凸包内删除的面
{
F[now].ok=0;
deal(p,F[now].b,F[now].a);
deal(p,F[now].c,F[now].b);
deal(p,F[now].a,F[now].c);
}
bool same(int s,int t)
{
Point &a=P[F[s].a];
Point &b=P[F[s].b];
Point &c=P[F[s].c];
return fabs(volume(a,b,c,P[F[t].a]))<eps &&
fabs(volume(a,b,c,P[F[t].b]))<eps &&
fabs(volume(a,b,c,P[F[t].c]))<eps;
}
//构建三维凸包
void create()
{
int i,j,tmp;
face add;
num=0;
if(n<4)return;
//**********************************************
//此段是为了保证前四个点不共面
bool flag=true;
for(i=1;i<n;i++)
{
if(vlen(P[0]-P[i])>eps)
{
swap(P[1],P[i]);
flag=false;
break;
}
}
if(flag)return;
flag=true;
//使前三个点不共线
for(i=2;i<n;i++)
{
if(vlen((P[0]-P[1])*(P[1]-P[i]))>eps)
{
swap(P[2],P[i]);
flag=false;
break;
}
}
if(flag)return;
flag=true;
//使前四个点不共面
for(int i=3;i<n;i++)
{
if(fabs((P[0]-P[1])*(P[1]-P[2])^(P[0]-P[i]))>eps)
{
swap(P[3],P[i]);
flag=false;
break;
}
}
if(flag)return;
//*****************************************
for(i=0;i<4;i++)
{
add.a=(i+1)%4;
add.b=(i+2)%4;
add.c=(i+3)%4;
add.ok=true;
if(dblcmp(P[i],add)>0)swap(add.b,add.c);
g[add.a][add.b]=g[add.b][add.c]=g[add.c][add.a]=num;
F[num++]=add;
}
for(i=4;i<n;i++)
{
for(j=0;j<num;j++)
{
if(F[j].ok&&dblcmp(P[i],F[j])>eps)
{
dfs(i,j);
break;
}
}
}
tmp=num;
for(i=num=0;i<tmp;i++)
if(F[i].ok)
F[num++]=F[i];
}
//表面积
double area()
{
double res=0;
if(n==3)
{
Point p=cross(P[0],P[1],P[2]);
res=vlen(p)/2.0;
return res;
}
for(int i=0;i<num;i++)
res+=area(P[F[i].a],P[F[i].b],P[F[i].c]);
return res/2.0;
}
double volume()
{
double res=0;
Point tmp(0,0,0);
for(int i=0;i<num;i++)
res+=volume(tmp,P[F[i].a],P[F[i].b],P[F[i].c]);
return fabs(res/6.0);
}
//表面三角形个数
int triangle()
{
return num;
}
//表面多边形个数
int polygon()
{
int i,j,res,flag;
for(i=res=0;i<num;i++)
{
flag=1;
for(j=0;j<i;j++)
if(same(i,j))
{
flag=0;
break;
}
res+=flag;
}
return res;
}
//三维凸包重心
Point barycenter()
{
Point ans(0,0,0),o(0,0,0);
double all=0;
for(int i=0;i<num;i++)
{
double vol=volume(o,P[F[i].a],P[F[i].b],P[F[i].c]);
ans=ans+(o+P[F[i].a]+P[F[i].b]+P[F[i].c])/4.0*vol;
all+=vol;
}
ans=ans/all;
return ans;
}
//点到面的距离
double ptoface(Point p,int i)
{
return fabs(volume(P[F[i].a],P[F[i].b],P[F[i].c],p)/vlen((P[F[i].b]-P[F[i].a])*(P[F[i].c]-P[F[i].a])));
}
};
CH3D hull;
int main()
{
// freopen("in.txt","r",stdin);
// freopen("out.txt","w",stdout);
while(scanf("%d",&hull.n)==1) //输入点
{
for(int i=0;i<hull.n;i++)
{
scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z) //输入每个点;
}
hull.create(); //创建凸包
printf("%d\n",hull.polygon());//求凸包表面多边形个数
}
//给一个三维凸包,求重心到表面的最短距离
//模板题:三维凸包+多边形重心+点面距离
while(scanf("%d",&hull.n)==1)
{
for(int i=0;i<hull.n;i++)
{
scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z);
}
hull.create(); //建立凸包
Point p=hull.barycenter() //凸包重心;
double minn=1e20;
for(int i=0;i<hull.num;i++) //凸包面的个数
{
minn=min(minn,hull.ptoface(p,i)); //p点到第i个面的距离
}
printf("%.3lf\n",minn);
}
while(scanf("%d",&hull.n)==1)
{
for(int i=0;i<hull.n;i++)
{
scanf("%lf%lf%lf",&hull.P[i].x,&hull.P[i].y,&hull.P[i].z);
}
hull.create();
printf("%.3f\n",hull.area());//凸包表面积
}
return 0;
}
给出n个点,用这些点构造三角形,然后找出相似三角形个数的最大值。
/*
数相似三角形个数
用map记录角度
注意:
1、去除重点
2、排除共线的情况
3、精度问题。
WA了8次了
*/
#include<stdio.h>
#include<math.h>
#include<map>
#include<iostream>
#include<string.h>
using namespace std;
const double eps=1e-6;
map<long long,int>mp;
struct Node
{
int x,y;
}node[50];
double dis(Node a,Node b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
void solve(Node a,Node b,Node c)
{
double la=dis(b,c);
double lb=dis(a,c);
double lc=dis(a,b);
if(la+lb<=lc+eps)return;
if(lb+lc<=la+eps)return;
if(la+lc<=lb+eps)return;
double A=acos((lb*lb+lc*lc-la*la)/(2*lb*lc));
double B=acos((la*la+lc*lc-lb*lb)/(2*la*lc));
double C=acos((la*la+lb*lb-lc*lc)/(2*la*lb));
if(A<eps||B<eps||C<eps)return;
int t1=(int)(A*10000);
int t2=(int)(B*10000);
int t3=(int)(C*10000);
if(t1>t2)swap(t1,t2);
if(t1>t3)swap(t1,t3);
if(t2>t3)swap(t2,t3);
if(t1==0)return;
long long t=t1*1000000*1000000+t2*1000000+t3;
mp[t]++;
}
int hole[220][220];
int main ()
{
int n;
while(scanf("%d",&n),n)
{
int t=0;
memset(hole,0,sizeof(hole));
for(int i=0;i<n;i++)
{
scanf("%d%d",&node[t].x,&node[t].y);
if(hole[node[t].x+100][node[t].y+100]==0)
{
hole[node[t].x+100][node[t].y+100]=1;
t++;
}
}
n=t;
mp.clear();
for(int i=0;i<n;i++)
for(int j=i+1;j<n;j++)
for(int k=j+1;k<n;k++)
solve(node[i],node[j],node[k]);
int ans=0;
map<long long,int>::iterator it;
for(it=mp.begin();it!=mp.end();it++)
{
int t=it->second;
if(t>ans)ans=t;
}
printf("%d\n",ans);
}
return 0;
}
题意:
给定n个点,现在要求找一个点作为一个半径为1的圆的圆心,使得这个圆能覆盖的点个数最多,输出最多能覆盖多少个点。
题解:
我们枚举两个点,认为他们是在这个圆上,那么就能确定圆心。之后可以从所有点到圆心的距离判断覆盖情况。
证:假设现在确定了一个能覆盖最多点的圆,那么我们可以移动这个圆使得覆盖的点数不变,且至少有两个点在这个圆上。
每两个点可以确定两个圆心,但只要用一个就够了(都用两点构成向量的上方圆心或者下方圆心)。
证明:若圆可以覆盖三个以上的点(两个点不需要找圆心。。),那么任意取覆盖的最外围三个点连线得到三角形ABC。我们枚举的时候枚举AB,AC,无论是先向量上方还是下方,都必定存在一个枚举的圆心,在三角形ABC内部。也就是说只要枚举一个圆心就能确定答案了(可以减少一半的时间消耗)。
注意:两点距离大于2的可以不枚举;初始化ans=1,否则一个点的时候会出错;距离判断时可以用距离的平方判断,因为sqrt很耗时间。
#include<stdio.h>
#include<string.h>
#include<algorithm>
#include<iostream>
#include<math.h>
using namespace std;
const int MAXN=330;
const double eps=1e-4;
double p[MAXN][2];
double xx1,yy1,xx2,yy2;
double dis(int i,int j)
{
return sqrt((p[i][0]-p[j][0])*(p[i][0]-p[j][0])+(p[i][1]-p[j][1])*(p[i][1]-p[j][1]));
}
void get_center_point(int a,int b)
{
double x0=(p[a][0]+p[b][0])/2;
double y0=(p[a][1]+p[b][1])/2;
double d=sqrt(1-((p[a][0]-p[b][0])*(p[a][0]-p[b][0])+(p[a][1]-p[b][1])*(p[a][1]-p[b][1]))/4);
if(fabs(p[a][1]-p[b][1])<1e-6)
{
xx1=x0;
yy1=y0+d;
xx2=x0;
yy2=y0-d;
}
else
{
double tmp=atan(-(p[a][0]-p[b][0])/(p[a][1]-p[b][1]));
double dx=d*cos(tmp);
double dy=d*sin(tmp);
xx1=x0+dx;
yy1=y0+dy;
xx2=x0-dx;
yy2=y0-dy;
}
}
int main()
{
int T;
int n;
scanf("%d",&T);
while(T--)
{
scanf("%d",&n);
for(int i=0;i<n;i++)
scanf("%lf%lf",&p[i][0],&p[i][1]);
int ans=1;
for(int i=0;i<n;i++)
for(int j=i+1;j<n;j++)
{
if(dis(i,j)>2.0)continue;
get_center_point(i,j);
int cnt=0;
for(int i=0;i<n;i++)
if(sqrt((p[i][0]-xx1)*(p[i][0]-xx1)+(p[i][1]-yy1)*(p[i][1]-yy1))<1+eps)
cnt++;
if(cnt>ans)ans=cnt;
cnt=0;
for(int i=0;i<n;i++)
if(sqrt((p[i][0]-xx2)*(p[i][0]-xx2)+(p[i][1]-yy2)*(p[i][1]-yy2))<1+eps)
cnt++;
if(cnt>ans)ans=cnt;
}
printf("%d\n",ans);
}
return 0;
}
题目:在三维空间中,给出太阳的位置(看成点光源),空间中存在一个凸物体,问物体在一个给定的平面内的影子大小。
分析:计算几何、凸包、三维坐标转换。求解分为四个步骤:1.物体顶点的位置判断;2.求物体顶点在平面上的投影;3.将投影点转化到x-y平面上;4.构造凸包、求面积。
1.物体顶点的位置判断。过光源做平行于已知平面的平面。有三种情况:
(1)所有点都在光源平面的上方或在光源平面内,此时没有影子;
(2)所有点都在光源下方,此时有一个可以计算的影子;
(3)其他情况,影子的面积为无限大。
至于判断方法,把点带入平面方程,根据结果的符号可以确定点在平面的哪一侧。注意给定的法向量方向,需要先转化到指向点的方向(如果点带入已知平面结果为负,则把abcd均乘以-1即可,完成平面翻转)。
2.求解顶点在平面上的投影。根据直线的上的已知点(x,y,z)和法向量(dx,dy,dz)可以得到直线点集:(x,y,z)+t(dx,dy,dz)。然后将点带入平面ax+by+cz=d即可求出t,确定直线与平面交点(x+tdx,y+tdy,z+tdz)。在求解之前先进点的位置判断,可以避免通过t的正负来确定交点是否存在(因为实际上是射线,t>=0)。
3.三维坐标转换。先绕z轴旋转(让y轴落在(a,b,0)方向上),然后绕x轴旋转即可(让z周落在(a,b,c)方向上)。下面给出坐标转换公式:
绕z轴旋转:x`=xcosC-ysinC;y`=xsinC+ycosC;z`=z;
绕x轴旋转:x`=x;y`=ycosA-zsinA;z`=ysinA+zcosA;
绕z轴旋转:x`=zsinB+xcosB;y`=y;z`=zcosB-xsinB;
推导过程并不复杂,因为有一维是不变的,可以直接看成平面的推导。
4.构造凸包求面积。这里直接看成是二维的即可(z全相等,可以忽略)。然后将凸包分割成三角形求利用叉乘面积即可。
#include <iostream>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <algorithm>
#include <vector>
using namespace std;
const int MAXN=110;
const double eps=1e-8;
const double PI=acos(-1.0);
inline int sign(double d)
{
if(d > eps)return 1;
else if(d < -eps)return -1;
else return 0;
}
struct Point //三维点坐标
{
double x,y,z;
Point(double xx=0,double yy=0,double zz=0):x(xx),y(yy),z(zz){}
Point operator +(const Point p1)
{
return Point(x+p1.x,y+p1.y,z+p1.z);
}
Point operator -(const Point p1)
{
return Point(x-p1.x,y-p1.y,z-p1.z);
}
Point operator *(const Point p1)
{
return Point(y*p1.z-z*p1.y,z*p1.x-x*p1.z,x*p1.y-y*p1.x);
}
Point operator *(double d)
{
return Point(x*d,y*d,z*d);
}
Point operator /(double d)
{
return Point(x/d,y/d,z/d);
}
double operator ^(const Point p1)
{
return x*p1.x+y*p1.y+z*p1.z;
}
double getLen()
{
return sqrt(x*x+y*y+z*z);
}
void read()
{
scanf("%lf%lf%lf",&x,&y,&z);
}
};
Point ps[MAXN];
inline Point get_point(Point st,Point ed,Point tp)//求点tp在直线st-ed上的垂足
{
double t1=(tp-st)^(ed-st);
double t2=(ed-st)^(ed-st);
double t=t1/t2;
Point ans=st + (ed-st)*t;
return ans;
}
inline double dist(Point st,Point ed)
{
return sqrt((ed-st)^(ed-st));
}
Point rotate(Point st,Point ed,Point tp,double A)//沿着直线st-ed旋转角度A,从ed往st看是逆时针
{
Point root=get_point(st,ed,tp);
Point e=(ed-st)/dist(ed,st);//单位向量
Point r=tp-root;
Point tmp=e*r;
Point ans=r*cos(A)+tmp*sin(A)+root;
return ans;
}
double a,b,c,d;
int n;
bool input()
{
scanf("%lf%lf%lf%lf",&a,&b,&c,&d);
if(a==0&&b==0&&c==0&&d==0)return false;
scanf("%d",&n);
for(int i=0;i<=n;i++)
ps[i].read();
return true;
}
//****************************************************
//二维凸包部分
struct point
{
double x,y;
};
point list[MAXN];
int stack[MAXN],top;
double cross(point p0,point p1,point p2)//p0p1 X p0p2
{
return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}
double dis(point p1,point p2)
{
return sqrt((p2.x-p1.x)*(p2.x-p1.x)+(p2.y-p1.y)*(p2.y-p1.y));
}
bool cmp(point p1,point p2)
{
double tmp=cross(list[0],p1,p2);
if(sign(tmp)>0)return true;
else if(sign(tmp)==0&&dis(list[0],p1)<dis(list[0],p2))return true;
else return false;
}
void init(int n)
{
point p0;
int k;
p0.x=list[0].x;
p0.y=list[0].y;
k=0;
for(int i=1;i<n;i++)
{
if(p0.y>list[i].y || ((p0.y==list[i].y)&&(p0.x>list[i].x)))
{
p0.x=list[i].x;
p0.y=list[i].y;
k=i;
}
}
list[k]=list[0];
list[0]=p0;
sort(list+1,list+n,cmp);
}
void graham(int n)
{
int i;
if(n==1)
{
top=0;
stack[0]=0;
return;
}
if(n==2)
{
top=1;
stack[0]=0;
stack[1]=1;
return;
}
for(int i=0;i<n;i++)stack[i]=i;
top=1;
for(i=2;i<n;i++)
{
while(top>0 && cross(list[stack[top-1]],list[stack[top]],list[i])<=0)top--;
top++;
stack[top]=i;
}
}
//****************************************************
void solve()
{
if(n<=2)
{
printf("0.00\n");
return;
}
Point p1(a,b,c);//法向量
Point tp(0,0,0);
Point e(0,0,1);
if(sign(a)!=0)tp.x=d/a;
else if(sign(b)!=0)tp.y=d/b;
else if(sign(c)!=0)tp.x=d/c;
ps[n+1]=tp;
Point vec=p1*e;
double A=(e^p1)/p1.getLen();
A=acos(A);
if(sign(A)!=0 && sign(A-PI)!=0)
{
Point p0(0,0,0);
for(int i=0;i<=n+1;i++)
ps[i]=rotate(p0,vec,ps[i],A);
}
for(int i=0;i<=n+1;i++)
ps[i].z-=ps[n+1].z;
if(sign(ps[n].z)<0)
{
for(int i=0;i<=n;i++)
ps[i].z=-ps[i].z;
}
int cnt1=0;//记录ps[n]上方的点的个数
int cnt2=0;//记录和ps[n]高度一样的点的个数
for(int i=0;i<n;i++)
{
int k=sign(ps[i].z-ps[n].z);
if(k>0)cnt1++;
else if(k==0)cnt2++;
}
if(cnt1+cnt2==n)
{
printf("0.00\n");
return;
}
if(cnt1>0 || cnt2>0)
{
printf("Infi\n");
return;
}
for(int i=0;i<n;i++)
{
double u=ps[n].z/(ps[i].z-ps[n].z);
list[i].x=ps[n].x+u*(ps[i].x-ps[n].x);
list[i].y=ps[n].y+u*(ps[i].y-ps[n].y);
}
init(n);
graham(n);
double ans=0;
for(int i=0;i<top;i++)
ans+=(list[stack[i]].x * list[stack[i+1]].y - list[stack[i+1]].x*list[stack[i]].y);
ans+=(list[stack[top]].x * list[stack[0]].y - list[stack[0]].x*list[stack[top]].y);
printf("%.2lf\n",ans/2);
}
int main()
{
while(input())solve();
return 0;
}
两园面积交
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <string>
#include <math.h>
#include <stdlib.h>
#include <time.h>
using namespace std;
const double eps = 1e-8;
const double PI = acos(-1.0);
struct Point
{
double x,y;
void input()
{
scanf("%lf%lf",&x,&y);
}
};
double dist(Point a,Point b)
{
return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}
//两个圆的公共部分面积
double Area_of_overlap(Point c1,double r1,Point c2,double r2)
{
double d = dist(c1,c2);
if(r1 + r2 < d + eps)return 0;
if(d < fabs(r1 - r2) + eps)
{
double r = min(r1,r2);
return PI*r*r;
}
double x = (d*d + r1*r1 - r2*r2)/(2*d);
double t1 = acos(x / r1);
double t2 = acos((d - x)/r2);
return r1*r1*t1 + r2*r2*t2 - d*r1*sin(t1);
}
int main()
{
//freopen("in.txt","r",stdin);
//freopen("out.txt","w",stdout);
Point c1,c2;
double r1,r2;
int T;
scanf("%d",&T);
int iCase = 0;
while(T--)
{
iCase ++;
c1.input();
scanf("%lf",&r1);
c2.input();
scanf("%lf",&r2);
printf("Case %d: %.6lf\n",iCase,Area_of_overlap(c1,r1,c2,r2));
}
return 0;
}
给你两条直线,问你最短距离是多少,并且求出最短距离的点,需要求得点是唯一的。
#pragma comment(linker, "/STACK:1024000000,1024000000")
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <string>
#include <math.h>
#include <stdlib.h>
#include <time.h>
using namespace std;
struct Point
{
double x,y,z;
Point(double _x = 0,double _y = 0,double _z = 0)
{
x = _x;
y = _y;
z = _z;
}
Point operator +(const Point &b)const
{
return Point(x + b.x, y + b.y,z + b.z);
}
Point operator -(const Point &b)const
{
return Point(x - b.x, y - b.y,z - b.z);
}
Point operator /(double k)
{
return Point(x/k,y/k,z/k);
}
Point operator *(double k)
{
return Point(x*k,y*k,z*k);
}
double operator *(const Point &b)const
{
return x*b.x + y*b.y + z*b.z;
}
Point operator ^(const Point &b)const
{
return Point(y*b.z - z *b.y,z*b.x - x*b.z, x*b.y - y * b.x);
}
void input()
{
scanf("%lf%lf%lf",&x,&y,&z);
}
void output()
{
printf("%.6lf %.6lf %.6lf",x,y,z);
}
};
double dis(Point a,Point b)
{
return sqrt((a.x - b.x) * (a.x- b.x) + (a.y - b.y) *(a.y - b.y) + (a.z - b.z)*(a.z - b.z));
}
double norm(Point a)
{
return sqrt(a.x *a.x + a.y *a.y + a.z * a.z);
}
Point A,B,C,D;
Point mid1,mid2;
int main()
{
//freopen("in.txt","r",stdin);
//freopen("out.txt","w",stdout);
int T;
scanf("%d",&T);
while(T--)
{
A.input();
B.input();
C.input();
D.input();
Point p1 = (B-A);
Point p2 = (D-C);
Point p = p1^p2;
double dd = (p*(A-C))/norm(p);
dd = fabs(dd);
printf("%.6lf\n",dd);
double t1 = ( (C-A)^p2 )*(p1^p2);
t1 /= norm(p1^p2)*norm(p1^p2);
double t2 = ( (C-A)^p1 )*(p1^p2);
t2 /= norm(p1^p2)*norm(p1^p2);
mid1 = A + (p1 * t1);
mid2 = C + (p2 * t2);
mid1.output();
printf(" ");
mid2.output();
printf("\n");
}
return 0;
} /*
题目意思:
给了平面上的n条线段:
让你在x轴上[0,L]的范围内找一个点作为圆心,画一个圆,这个圆不能把线段包含在里面去。
求最大的半径。
求解:
二分最终的答案,求解。
对于二分的半径值,对每条线段找出对于x位置上满足半径要求的段。
看n条线段有没有交集就可以了
*/
#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <string>
#include <math.h>
#include <stdlib.h>
#include <time.h>
using namespace std;
const double eps = 1e-8;
const double inf = 1e5;
const double pi = acos(-1.0);
int sgn(double x)
{
if(fabs(x) < eps)return 0;
else if(x < 0)return -1;
return 1;
}
struct Point
{
double x,y;
Point(){}
Point(double _x,double _y)
{
x = _x;
y = _y;
}
void input()
{
scanf("%lf%lf",&x,&y);
}
Point operator -(const Point &b)const
{
return Point(x-b.x,y-b.y);
}
double operator *(const Point &b)const
{
return x*b.x + y*b.y;
}
double operator ^(const Point &b)const
{
return x*b.y - y*b.x;
}
double len()
{
return hypot(x,y);
}
double distance(Point p)
{
return hypot(x-p.x,y-p.y);
}
};
struct Line
{
Point s,e;
Line(){}
Line(Point _s,Point _e)
{
s = _s;
e = _e;
}
void input()
{
s.input();
e.input();
}
double length()
{
return s.distance(e);
}
double dispointtoline(Point p)
{
return fabs((p-s)^(e-s))/length();
}
double dispointtoseg(Point p)
{
if(sgn((p-s)*(e-s)) < 0 || sgn((p-e)*(s-e)) < 0)
return min(p.distance(s),p.distance(e));
return dispointtoline(p);
}
};
const int MAXN = 2020;
Line line[MAXN];
double get(int i,double L)
{
double l = 0, r = L;
while(r - l >= eps)
{
double mid = (l + r)/2;
double midmid = (r + mid)/2;
if(line[i].dispointtoseg(Point(mid,0)) < line[i].dispointtoseg(Point(midmid,0)))
r = midmid - eps;
else l = mid + eps;
}
return l;
}
int n;
double L;
struct Node
{
double a;
int c;
Node(double _a = 0,int _c = 0)
{
a = _a;
c = _c;
}
};
bool cmp(Node a,Node b)
{
if(sgn(a.a - b.a) != 0)return a.a < b.a;
else return a.c > b.c;
}
double bin1(int i,double l,double r,double R)
{
while(r-l >= eps)
{
double mid = (l+r)/2;
if(line[i].dispointtoseg(Point(mid,0)) <= R)
r = mid - eps;
else l = mid + eps;
}
return l;
}
double bin2(int i,double l,double r,double R)
{
while(r-l >= eps)
{
double mid = (l+r)/2;
if(line[i].dispointtoseg(Point(mid,0)) <= R)
l = mid + eps;
else r = mid - eps;
}
return l;
}
bool gao(double r)
{
vector<Node>vec;
for(int i = 0;i < n;i++)
{
double tmp = get(i,L);
if(sgn(line[i].dispointtoseg(Point(tmp,0)) - r) >= 0)
{
vec.push_back(Node(0,1));
vec.push_back(Node(L,-1));
continue;
}
if(sgn(line[i].dispointtoseg(Point(0,0)) - r) >= 0)
{
double tt = bin1(i,0,tmp,r);
vec.push_back(Node(0,1));
vec.push_back(Node(tt,-1));
}
if(sgn(line[i].dispointtoseg(Point(L,0)) - r) >= 0)
{
double tt = bin2(i,tmp,L,r);
vec.push_back(Node(tt,1));
vec.push_back(Node(L,-1));
}
}
sort(vec.begin(),vec.end(),cmp);
int cnt = 0;
for(int i = 0;i < vec.size();i++)
{
cnt += vec[i].c;
if(cnt >= n)return true;
}
return false;
}
double solve()
{
double l = 0, r = inf;
while(r-l >= eps)
{
double mid = (l+r)/2;
if(gao(mid))l = mid+eps;
else r = mid - eps;
}
return l;
}
int main()
{
//freopen("in.txt","r",stdin);
//freopen("out.txt","w",stdout);
int T;
scanf("%d",&T);
while(T--)
{
scanf("%d%lf",&n,&L);
for(int i = 0;i < n;i++)
line[i].input();
printf("%.3lf\n",solve());
}
return 0;
}
计算几何模板
最新推荐文章于 2018-02-16 22:50:13 发布