Library of Geometry

Geometry.cpp

#include <iostream>
#include <algorithm>
#include <vector>
#include <cmath>
#include <limits.h>
#include <float.h>

#define eps 1e-10
#define M INT_MAX
using namespace std;

double add(double a, double b){
    if(abs(a + b) < eps * (abs(a)) + abs(b)) return 0;
    else return a + b;
}

struct P{
    double x, y;
    P() {}
    P(double x, double y) : x(x), y(y){
    }
    P operator + ( P t){
	return P(add(x, t.x), add(y, t.y)); 
    }
    P operator - (P t){
	return P(add(x, -t.x) , add(y, -t.y));
    }
    P operator * (double d){
	return P(x * d , y * d); 
    }
};

//内積を求める a・b
double dot(P a, P b){
    return add(a.x * b.x, a.y * b.y);
}

//外積を求める a×b
double det(P a, P b){
    return add(a.x * b.y, -a.y * b.x);
}

//距離を求める
double dis(P t){
    return sqrt(t.x*t.x+t.y*t.y);
}

//単位ベクトルを求める
P unit_vector(P t){
	double u=fabs(dis(t));
	return P(t.x/u , t.y/u);
}

//cosθを求める
double getcos(P a , P b){
	return (dot(a,b) / (dis(a)*dis(b)));
}

//sinθを求める
double getsin(P a , P b){
	double t=getcos(a,b);
	return sqrt(1.0-t*t);
}

//点pが直線ab上にあるかどうか
bool pol(P p , P a , P b){
	return ( det( b-a , p-a ) < eps );
}

//点pが線分ab上にあるかどうか
bool pos(P p , P a , P b){
	return ( fabs(dis(a-p)) + fabs(dis(p-b)) < fabs(dis(a-b)) + eps);
}

//点pと直線abとの距離を求める
double pld(P p , P a , P b){
	return fabs(det( b-a , p-a ))/dis( b-a );
}

//点pと線分abとの距離を求める
double psd(P p , P a , P b){
	if( dot( b-a , p-a ) < eps) return fabs(dis(p-a));
	if( dot( a-b , p-b ) < eps) return fabs(dis(p-b));
	return fabs(det( b-a , p-a )) / fabs(dis(b-a));
}

//線分交差判定
bool intersect_s(P a1 , P a2 , P b1 , P b2){
	return ( det(a2-a1 , b1-a1)*det(a2-a1 , b2-a1) < eps) &&
		   ( det(b2-b1 , a1-b1)*det(b2-b1 , a2-b1) < eps);
}

//線分交点計算
P interpoint_s(P a1 , P a2 , P b1 , P b2){
	P b=b2-b1;
	double d1=fabs(det(b , a1-b1));
	double d2=fabs(det(b , a2-b1));
	double t=d1 / (d1 + d2);
	P a=a2-a1;
	P v=a*t;
	return a1+v;
}

//直線交差判定
bool intersect_l(P a1 , P a2 , P b1 , P b2){
	return (det(a1-a2 , b1-b2) < eps);
}

//直線交点計算
P interpoint_l(P a1 , P a2 , P b1 , P b2){
	P a=a2-a1;
	P b=b2-b1;
	double t=det(b , b1-a1) / det(b , a);
	P v=a*t;
	return a1+v;
}

//直線abと円及び球 |x-c|=rの交点を求める 　球の場合は型P→Sに変更するだけ
void interpoint_lc(P a , P b , P c , double r , P ans[]){
	if(pld(c , a , b) > r + eps) return;
	P v=unit_vector(b-a);
	double delta=dot(v,a-c)*dot(v,a-c)-dis(a-c)*dis(a-c)+r*r;
	double t=-dot(v,a-c);
	double s=sqrt(delta);
	ans[0]=a+v*(t+s);
	ans[1]=a+v*(t-s);
}

//１次変換集

//ｘ（y=k）に関する対象変換 　k=0でｘ軸による変換
P x_translate(P t,double k){
	return P(t.x , 2*k-t.y);
}

//ｙ(x=k)に関する対象変換	k=0でｙ軸による変換
P y_translate(P t,double k){
	return P(2*k-t.x , t.y);
}

//点P kに関する対象変換		P(0,0)で原点による変換
P o_translate(P t,P k){
	return k+(k-t);
}

//点pを中心としてr(radian)回転		p(0,0)で原点を中心として回転
P rotate(P t , P p , double r){
	//double r=radians(angle);
	double ta=cos(r)*(t.x-p.x)-sin(r)*(t.y-p.y)+p.x;
	double tb=sin(r)*(t.x-p.x)+cos(r)*(t.y-p.y)+p.y;
	return P(ta , tb);
}

//応用
//２円 |x-a|=raと|x-b|=rbの交点計算
void interpoint_cc(P a , double ra , P b , double rb , P ans[]){
	double di=fabs(dis(a-b));
	if(di > ra+rb || di < fabs(ra-rb)) return;
	double t=(ra*ra-rb*rb+di*di)/(di+di);
	double rd=acos(t/ra);

	P dv=unit_vector(b-a);
	P g1=rotate(dv , P(0,0) , rd);
	P g2=rotate(dv , P(0,0) , -rd);
	ans[0]=a+g1*ra;
	ans[1]=a+g2*ra;
}

//法線ベクトルを求める
P normal_vector(P p,P a,P b){
	P v=unit_vector(b-a);
	v = det(v , p-a) > 0 ?  P(v.y,(-1)*v.x) : P((-1)*v.y , v.x);
	return v*pld(p,a,b);
}

//直線abに関する対象変換
P f_translate(P t , P a , P b){
	return t+normal_vector(t,a,b)*2;
}

//３角形の面積を求める
double area(P a,P b,P c){
	return fabs(det(c-a , b-a)*0.5);
}

//多角形の面積を求める //凸包のソート済みが前提
double polygon_area(vector<P> t){
	double ans=0.0;
	for(unsigned int i=0;i<t.size();i++)
		ans+=det(t[i] , t[(i+1)%t.size()]);
	return ans/2;
}

//凸包
vector<P> conhel(vector<P> t){
	vector<P> left,right;
	bool jd[10000];
	fill(jd ,&jd[10000],false);
	
	P mio,mao;
	double midy=DBL_MAX,midx=DBL_MAX;
	double mady=DBL_MIN,madx=DBL_MIN;
	int mik=0,mak=0;

	for(unsigned int i=0;i<t.size();i++){
		if(t[i].y < midy)
			mio=t[i], midy=t[i].y, midx=t[i].x , mik=i;
		else if(fabs(t[i].y-midy) < eps && t[i].x < midx)
			mio=t[i], midy=t[i].y, midx=t[i].x , mik=i;
		if(t[i].y > mady)
			mao=t[i], mady=t[i].y, madx=t[i].x , mak=i;
		else if(fabs(t[i].y-mady) < eps && t[i].x > madx)
			mao=t[i], mady=t[i].y, madx=t[i].x , mak=i;
	}
	right.push_back(mio);
	left.push_back(mao);
	bool ld=false , rd=false;

	while(1){
		int mig=-1,mag=-1;
		double mih=DBL_MAX,mah=DBL_MAX;
		for(unsigned int i=0;i<t.size();i++){
			if(jd[i] || fabs(dis(t[i]-t[mik])) < eps || fabs(dis(t[i]-t[mak])) < eps) 
				continue;
			double mir=acos(getcos(t[i]-t[mik] , P(1,0)));
			double mar=acos(getcos(t[i]-t[mak] , P(-1,0)));

			if(mih > mir) mih=mir, mig=i;
			if(mah > mar) mah=mar, mag=i;
		}
		jd[mig]=true , jd[mag]=true;
		if(mig==-1) rd=true;
		if(mag==-1) ld=true;
		if(!rd)
			if(fabs(dis(t[mig]-mao)) >= eps)  right.push_back(t[mig]);
			else rd=!rd;
		if(!ld)
			if(fabs(dis(t[mag]-mio)) >= eps)  left.push_back(t[mag]);
			else ld=!ld;
		if(rd && ld) break;
		mik=mig , mak=mag;
		
	}
	for(unsigned int i=0;i<left.size();i++)
		right.insert(right.end() , left[i]);

	return right;
}