#include<iostream>
#include<math.h>
#include<vector>
#include<algorithm>
#include<stdio.h>
using namespace std;
const double eps = 1e-10;
const double pi = acos( -1.0 );
const double inf = 1e100;
inline int sig( double x ){
return (x > eps) - (x < - eps);
}
struct point{
double x, y, dis, angle;
point(){}
point( double x0, double y0 ):x(x0),y(y0){}
//逻辑运算
inline bool operator < ( point t )const {
return sig ( y - t.y ) < 0 || sig ( y - t.y ) == 0 && sig( x - t.x ) < 0 ;
}
inline bool operator == ( point t )const{
return sig ( x - t.x ) == 0 && sig ( y - t.y ) == 0;
}
//算术运算
inline point operator - ( point b){
return point( x - b.x , y - b.y );
}
inline point operator + ( point b ){
return point( x + b.x, y + b.y );
}
inline point operator * ( double t ){
return point ( x * t , y * t );
}
inline point operator / ( double t ){
return point ( x / t , y / t );
}
//求长度
inline double len(){
return sqrt( x * x + y * y );
}
//长度的平方
inline double len2(){
return x * x + y * y ;
}
inline double dist( point t ){
return sqrt( (x - t.x)*(x - t.x) + (y - t.y)*(y - t.y) );
}
//输入输出
int read(){
return scanf("%lf%lf",&x,&y);
}
void out(){
printf("(%.2lf %.2lf) ",x,y);
}
};
inline bool _cmp( point a, point b ){
return sig( a.angle - b.angle ) < 0 || sig( a.angle - b.angle ) == 0 && sig( a.dis - b.dis ) < 0;
}
// 点积
inline double dot( point a, point b ){
return a.x * b.x + a.y * b.y ;
}
inline double dot( point p , point a, point b ){//这种形式更常用
return ( a.x - p.x ) * (b.x - p.x) + (a.y - p.y) * (b.y - p.y);
}
//叉积 ,遵循 左手定则（左旋为正）
inline double cross( point a, point b){
return a.y * b.x - a.x * b.y ;
}
inline double cross( point p , point a, point b){//这种形式更常用
return (a.y - p.y) *(b.x - p.x) - (a.x - p.x) * (b.y - p.y);
}
//angle > 0 表示逆时针旋转，angle < 0 为顺时针
//矢量pv以p为顶点逆时针旋转angle并放大scale倍
inline point rotate( point v, point p, double angle, double scale ){
point ret = p;
v.x -= p.x, v.y -= p.y;
p.x = scale * cos(angle); p.y = scale * sin(angle);
ret.x += v.x * p.x - v.y * p.y;
ret.y += v.x * p.y + v.y * p.x;
return ret;
}
//向量 p 的法向量 ，p逆时针旋转90度
inline point faxiang( point p){
return point( - p.y , p.x );
}
//向量 ab 左移 r 长度
inline point left( point a, point b, double r){
point ta, tb, tt;
tt.x = b.y - a.y; tt.y = a.x - b.x;
double k = r / sqrt( tt.x * tt.x + tt.y * tt.y );
tt.x = - tt.x * k; tt.y = - tt.y * k;
return tt;
}
//两向量之间的夹角[0,pi] , 乘以 cross( a, b ) 后便成了有向夹角。
inline double point_angle(point a,point b){
return acos( dot( a, b ) / sqrt( a.len2() * b.len2() ));
return cross( a, b ) * acos( dot( a, b ) / sqrt( a.len2() * b.len2() ));
}
//点到直线的距离
inline double point_to_line(point p, point a, point b){
return fabs( cross( a, p, b ) / ( b - a ).len());
}
//返回点p与线段ab的距离,
inline double point_to_seg(point p , point a,point b){
if ( sig ( dot ( a, p, b ) ) < 0 )//p在a端外
return ( p - a ).len();
if ( sig ( dot ( b, p, a ) ) < 0 )//p在b端外
return ( p - b ).len();
return point_to_line( p , a , b );//p正对于ab
}
//返回线段ab与线段cd之间的距离，相交（有一个公共点）返回0,交点为p。自己通过<,<=,来处理边界
inline double seg_to_seg ( point a , point b ,point c ,point d,point &p ){
if( ! ( (sig ( max( a.x , b.x ) - min( c.x , d.x ) ) < 0 ) || ( sig ( max( d.x , c.x ) - min( a.x , b.x ) ) < 0 )
||( sig ( max ( a.y , b.y ) - min( c.y , d.y ) ) < 0 ) || ( sig ( max( d.y , c.y ) - min ( a.y , b.y ) ) < 0 ) ) )
{
if ( sig ( cross( a, b, c ) * cross ( a, b, d) ) <= 0 && sig ( cross( c, d, a ) * cross ( c, d, b ) ) <= 0 ){
p = a;
double t=((a.x-c.x)*(c.y-d.y)-(a.y-c.y)*(c.x-d.x))
/((a.x-b.x)*(c.y-d.y)-(a.y-b.y)*(c.x-d.x));
p.x+=(b.x-a.x)*t;
p.y+=(b.y-a.y)*t;
return 0;
}
}
//不相交
return min( min( point_to_seg( a, c, d ) , point_to_seg( b, c, d ) ), min( point_to_seg( c, a, b), point_to_seg( d, a, b ) ));
}
// 严格相交
inline int seg_insert_line ( point a, point b, point c, point d, point & p){
if ( sig ( cross( c, d, a ) * cross ( c, d, b ) ) < 0 ){
p = a;
double t=((a.x-c.x)*(c.y-d.y)-(a.y-c.y)*(c.x-d.x))
/((a.x-b.x)*(c.y-d.y)-(a.y-b.y)*(c.x-d.x));
p.x+=(b.x-a.x)*t;
p.y+=(b.y-a.y)*t;
return 1;
}
return 0;
}
// 求向量 a 的极角，角度范围为 [0,pi*2)
double angle( point a ){
double t = atan2( a.y, a.x );
return sig( t ) >= 0 ? t : t + pi * 2;
}
// 角度修正函数，根据不同的题目有不同写法和功能
void fix_angle( double & rad ){
if( sig( rad - pi ) > 0 ) rad -= pi * 2;
if( sig( rad + pi ) < 0 ) rad += pi * 2;
}
//求直线和园相交。0为相离，1为相切，2为相交，交点v1，v2 ,且 v2到v1是逆时针指向
int cir_insert_line( point c, double r, point a, point b, point &v1, point &v2){
double h = fabs( cross( a, c, b )/(b - a).len());
if( sig( h - r ) > 0 )
return 0;
point tmp = point ( -(b - a).y, (b - a).x );
if ( sig( cross( a, c, b )) > 0 )
tmp.x = - tmp.x , tmp.y = -tmp.y;
tmp = tmp*( h / tmp.len() );
point mid = c + tmp;
if ( sig( h - r ) == 0 ){
v1 = mid;
return 1;
}
double l2 = ( r * r - h * h );
tmp = ( b - a );
tmp = tmp * sqrt( l2 / tmp.len2() );
v1 = tmp + mid; tmp.x = - tmp.x;
tmp.y = - tmp.y; v2 = tmp + mid;
return 2;
}
// 判断点在线段上，条件：pt,pl1,pl2共线
inline bool pt_in_seg( point pt, point pl1, point pl2 ){
return ( sig( dot(pl1,pt,pl2)) >= 0 && sig( dot(pl2,pt,pl1) ) >= 0 );
}
// 判断点在线段上（不含端点），条件：pt,pl1,pl2共线
inline bool pt_in_segs( point pt, point pl1, point pl2 ){
return ( sig( dot(pl1,pt,pl2)) > 0 && sig( dot(pl2,pt,pl1) ) > 0 );
}
inline int pt_line( point p1, point p2 , point &ret ){
ret.y = p2.x - p1.x; ret.x = p1.y - p2.y;
ret.dis = - dot( p1, ret );
return 0;
}
//求直线和园相交。0为相离，1为相切，2为相交，交点v1，v2 ,且 v2到v1是逆时针指向
int lcins(point po,double r,point pl1,point pl2,point &ret1,point &ret2)
{
point vl;
double tmp1, tmp2;
pt_line( pl1, pl2, vl );
tmp1 = dot( vl, po ) + vl.dis;
double ab = vl.len2(), mdis2 = ab * r * r - tmp1 * tmp1 ;
if( sig( mdis2 ) < 0 )
return 0;
else
{
tmp1 = ( vl.y * cross( vl, po ) - vl.dis * vl.x );
tmp2 = ( vl.x * cross( po, vl ) - vl.dis * vl.y );
if ( sig( mdis2 ) <= 0 )
{
ret1.x = tmp1 / ab ; ret1.y = tmp2 / ab ;
return 1;
}
else
{
double sq = sqrt( mdis2 );
ret1.x = ( tmp1 + vl.y * sq ) / ab; ret1.y = ( tmp2 - vl.x * sq ) / ab;
ret2.x = ( tmp1 - vl.y * sq ) / ab; ret2.y = ( tmp2 + vl.x * sq ) / ab;
return 2;
}
}
}
// 半径为r，圆心在原点的圆，与三角形 ps，pe，o的相交面积。
double area_c2p( double r, point ps, point pe ){
double d1 = ps.len2(), d2 = pe.len2(), r2 = r * r;
double rad1, rad2, drad;
if ( d1 <= r2 && d2 <= r2 ) return cross( pe, ps ) * .5;
point po(0,0), ins1, ins2;
rad1 = angle( ps ); rad2 = angle( pe );
drad = rad2 - rad1; fix_angle( drad );
int nins = cir_insert_line( po, r, ps, pe, ins1, ins2 );
if( nins <= 1 ) return drad*.5*r*r;
if(( sig( dot( ps, ins1, pe ) ) < 0 && sig( dot( ps, ins2, pe) ) < 0 )
||( sig( dot( pe, ins1, ps) ) < 0 && sig( dot( pe, ins2, ps ) ) < 0 ))
return drad*.5*r*r;
else if( pt_in_seg( ins1, ps, pe ) && pt_in_seg( ins2, ps, pe ))
{
double radi1, radi2, drad2, ret;
if ( pt_in_segs( ins1, ps, ins2 ))
{
radi1 = angle( ins1 ); radi2 = angle( ins2 );
ret = cross( ins2, ins1 ) * .5;
}
else
{
radi1 = angle( ins2 ); radi2 = angle( ins1 );
ret = cross( ins1, ins2 ) * .5;
}
drad = radi1 - rad1;
fix_angle( drad );
drad2 = rad2 - radi2;
fix_angle( drad2);
return ret + ( drad + drad2 ) * .5 * r * r;
}
else
{
point pins, pout;
double radi;
if ( sig( dot( ps, ins1, pe )) < 0
|| sig( dot( pe, ins1, ps )) < 0 )
{
pins = ins2, pout = ins1;
}
else
{
pins = ins1, pout = ins2;
}
radi = angle( pins );
if ( sig( dot( pe, pout, ps ) ) < 0 )
{
drad = radi - rad1; fix_angle( drad );
return cross ( pe, pins ) * .5 + drad * .5 * r * r;
}
else
{
drad = rad2 - radi; fix_angle( drad );
return cross( pins, ps ) * .5 + drad * .5 * r * r;
}
}
}
//半径为r圆，圆心在原点上的圆与多边形的相交面积
double area_cpoly( double r, point p[], int n ){
int i; double ret = 0;
p[n] = p[0];
for ( i = 1 ; i <= n ; i ++ )
ret += area_c2p( r, p[i-1], p[i] );
return ret;
}
//求圆和圆相交，返回交点个数, 两圆重合返回 -1
inline int cir_insert_cir( point c1, double r1,point c2, double r2, point &p1, point &p2 ){
double dist = c1.dist( c2 ) , tmp;
if ( sig ( dist - r1 - r2 ) > 0 )// 外离
return 0;
double r_min = min( r1, r2 ), r_max = max( r1, r2 );
if ( dist < r_max && sig ( dist + r_min - r_max ) < 0 )
return 0;// 内离
if ( sig( dist - r1 - r2 ) == 0 ){// 外切
p1 = c1 + (c2 - c1)*(r1)/(dist);
return 1;
}
if ( dist < r_max && sig ( dist + r_min - r_max ) == 0 ){// 内切
if( sig( dist ) == 0 )
return -1; // 两圆重合
p1 = c1 + (c2 - c1)*(r_max)/( dist );
return 1;
}
point u, v;
double t;
t = (1.+(r1*r1-r2*r2)/dist/dist)/2;
u.x = c1.x+(c2.x-c1.x)*t; u.y = c1.y+(c2.y-c1.y)*t;
v.x = u.x+c1.y-c2.y; v.y = u.y-c1.x+c2.x;
cir_insert_line(c1,r1,u,v,p1,p2);
return 2;
}
//求圆上两点之间的最小距离
inline double cir_dist( point a, point b, point c, double r){
double sita = point_angle( ( a - c ), ( b - c ) );
return r * sita;
}
//多边形
struct polygon{
static const int MAXN = 10510;//多边形大小
point p[MAXN] ;
int n ;
// 构造函数
polygon(){ n = 0 ; }
// 在多边形末尾加点
void add( point & t ){
this->p[ n ++ ] = t;
}
// 求面积，逆时针为正，顺时针为负
inline double area(){
double sum = 0 ;
p[n] = p[0];
for ( int i = 0 ; i < n ; i++ )
sum += cross( p[ i + 1 ] , p[ i ] );
return sum * .5 ;
}
// 将多边形转为反向
inline int turn_back(){
double tmp = area();
if( sig( tmp ) < 0 )
for ( int i = 0; i < n / 2 ; i ++ )
swap( p[i], p[n - i - 1] );
return sig( tmp );
}
//求点集的凸包，存于 多边形con中 ,原多边形将乱序
inline int convex( polygon & con ){
int t = 0, i;
point tmp;
for( i = 1 ; i < n ; i ++ )
{
if( p[i] < p[t] )//先 y 值最小，再x值最小处理
t = i;
}
swap( p[t], p[0] );
for( i = 0 ; i < n ; i ++ )
{
tmp = p[i] - p[0];
p[i].dis = tmp.len2();
p[i].angle = atan2( tmp.y, tmp.x );
}
sort( p, p + n, _cmp );
int k = 0;
con.p[ k ++ ] = p[n-1];
con.p[ k ++ ] = p[0];
if( sig( p[1].angle - p[n-1].angle ) == 0 )//凸包为一线段
con.p[ k ++ ] = p[ n - 1];
else
{
for ( i = 1 ; i < n ; i ++ )
if( sig( cross( con.p[k-1], con.p[k-2], p[i] )) > 0 )
con.p[ k ++ ] = p[i];
else
i -- , k -- ;
}
k -- ;
con.n = k;
return k;
}
// 求多边形重心
inline point tri_centroid( point a, point b, point c ){
point t( 0, 0 );
t.x = ( a.x + b.x + c.x ) / 3 ;
t.y = ( a.y + b.y + c.y ) / 3 ;
return t;
}
inline point poly_centroid()
{
int i;
double area = 0, tmp ;
point rel( 0, 0 );
for ( i = 1 ; i < n - 1 ; i ++ )
{
tmp = cross( p[0], p[i], p[ i + 1 ] );
area += tmp;
rel = ( tri_centroid( p[0], p[i], p[i+1] ) * tmp ) + rel;
}
if ( sig ( area ) != 0 )
rel = rel / area;
return rel;
}
inline int read(){
int t = scanf("%d",&n);
if ( t == EOF )
return t;
for( int i = 0 ; i < n ; i ++ )
p[ i ].read();
return n;
}
inline void out(){
printf("%d\n",n);
for ( int i = 0 ; i < n ; i ++ )
{
p[i].out();
cout<<endl;
}
}
// 直线st->ed，切割多边形，留下st->ed左边的部分，q[]为辅助数组。
void cut( point st, point ed, point q[] )
{
int i, j, k;
point tmp;
k = 0;
p[n] = p[0];
for ( j = 0 ; j < n ; j ++ )
{
if ( sig( cross( st, ed, p[j] ) ) <= 0 )
q[ k ++ ] = p[j];
if ( seg_insert_line( p[j], p[j+1], st, ed, tmp ) )
q[ k ++ ] = tmp;
}
n = k;
for ( j = 0 ; j < n ; j ++ ) p[j] = q[j];
}
// O(n^2) 的半平面求交，多边形上的点逆时针输入
polygon half_plane2(){
turn_back();
polygon tmp;
p[ n ] = p[ 0 ];
tmp = * this ;
point q[ MAXN ];
for ( int i = 0 ; i < n ; i ++ )
{
tmp.cut( p[i], p[i+1], q );
}
return tmp;
// tmp.n 为0的时候 ，无交集，即多边形的核
}
};
posted on 2009-08-24 22:58
Huicpc217 阅读(360)
评论(0) 编辑 收藏 引用 所属分类:
模板 、
计算几何