3维计算几何模板
#include <iostream>
#include <cmath>
#include <stdio.h>
#include <vector>
#include <string.h>
#include <stdlib.h>
#include <algorithm>
using namespace std;#define MAX_N 110/*------------------常量區-------------------*/const double INF = 1e10; // 無窮大
const double EPS = 1e-5; // 計算精度
const double PI = acos(-1.0);// PI
const int PINXING = 0; // 平行
const int XIANGJIAO = 1; // 相交
const int XIANGLI = 0; // 相離
const int GONGXIAN = 2; // 共線
const int CHONGDIE = -1; // 重疊
const int INSIDE = 1; // 點在圖形內部
const int OUTSIDE = 0; // 點在圖形外部
const int BORDER = 2; // 點在圖形邊界/*-----------------類型定義區----------------*/struct Point { // 二維點或矢量double x, y;//double angle, dis;
Point() {}Point(double x0, double y0): x(x0), y(y0) {}void read(){scanf("%lf%lf",&x,&y);}
};struct Line { // 二維的直線或線段
Point p1, p2;Line() {}Line(Point p10, Point p20): p1(p10), p2(p20) {}void read(){scanf("%lf%lf%lf%lf",&p1.x,&p1.y,&p2.x,&p2.y);}
};struct Rect { // 用長寬表示矩形的方法 w, h分別表示寬度和高度double w, h;Rect() {}Rect(double _w,double _h) : w(_w),h(_h) {}
};
struct Rect_2 { // 表示矩形，左下角坐標是(xl, yl)，右上角坐標是(xh, yh)double xl, yl, xh, yh;Rect_2() {}Rect_2(double _xl,double _yl,double _xh,double _yh) : xl(_xl),yl(_yl),xh(_xh),yh(_yh) {}
};
struct Circle { //圓
Point c;double r;Circle() {}Circle(Point _c,double _r) :c(_c),r(_r) {}
};typedef vector<Point> Polygon; // 二維多邊形
typedef vector<Point> Points; // 二維點集/*-------------------基本函數區---------------------*/inline double max(double x,double y)
{return x > y ? x : y;
}
inline double min(double x, double y)
{return x > y ? y : x;
}
inline bool ZERO(double x) // x == 0
{return (fabs(x) < EPS);
}
inline bool ZERO(Point p) // p == 0
{return (ZERO(p.x) && ZERO(p.y));
}inline bool EQ(double x, double y) // eqaul, x == y
{return (fabs(x - y) < EPS);
}
inline bool NEQ(double x, double y) // not equal, x != y
{return (fabs(x - y) >= EPS);
}
inline bool LT(double x, double y) // less than, x < y
{return ( NEQ(x, y) && (x < y) );
}
inline bool GT(double x, double y) // greater than, x > y
{return ( NEQ(x, y) && (x > y) );
}
inline bool LEQ(double x, double y) // less equal, x <= y
{return ( EQ(x, y) || (x < y) );
}
inline bool GEQ(double x, double y) // greater equal, x >= y
{return ( EQ(x, y) || (x > y) );
}// 輸出浮點數前，防止輸出-0.00調用該函數進行修正！
inline double FIX(double x)
{return (fabs(x) < EPS) ? 0 : x;
}/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
//-------------------3D 區域----------------------------//
struct Point3D { //三維點或矢量double x, y, z;Point3D() {}Point3D(double x0, double y0, double z0): x(x0), y(y0), z(z0) {}void read(){scanf("%lf%lf%lf",&x,&y,&z);}
};struct Line3D { // 三維的直線或線段
Point3D p1, p2;Line3D() {}Line3D(Point3D p10, Point3D p20): p1(p10), p2(p20) {}void read(){scanf("%lf%lf%lf%lf%lf%lf",&p1.x,&p1.y,&p1.z,&p2.x,&p2.y,&p2.z);}
};struct Area3D{Point3D p1,p2,p3;Area3D(){}Area3D(Point3D p10, Point3D p20,Point3D p30): p1(p10), p2(p20), p3(p30){}void read(){p1.read(); p2.read(); p3.read();//scanf("%lf%lf%lf%lf%lf%lf%lf%lf%lf",&p1.x,&p1.y,&p1.z,&p2.x,&p2.y,&p2.z,&p3.x,&p3.y,&p3.z);
}
};inline bool ZERO(Point3D p) // p == 0
{return (ZERO(p.x) && ZERO(p.y) && ZERO(p.z));
}//
//三維矢量運算
bool operator==(Point3D p1, Point3D p2)
{return ( EQ(p1.x, p2.x) && EQ(p1.y, p2.y) && EQ(p1.z, p2.z) );
}
bool operator<(Point3D p1, Point3D p2)
{if (NEQ(p1.x, p2.x)) {return (p1.x < p2.x);} else if (NEQ(p1.y, p2.y)) {return (p1.y < p2.y);} else {return (p1.z < p2.z);}
}
Point3D operator+(Point3D p1, Point3D p2)
{return Point3D(p1.x + p2.x, p1.y + p2.y, p1.z + p2.z);
}
Point3D operator-(Point3D p1, Point3D p2)
{return Point3D(p1.x - p2.x, p1.y - p2.y, p1.z - p2.z);
}
Point3D operator * (const Point3D& A, double p) { return Point3D(A.x*p, A.y*p, A.z*p); }
Point3D operator / (const Point3D& A, double p) { return Point3D(A.x/p, A.y/p, A.z/p); }Point3D operator*(Point3D p1, Point3D p2) // 計算叉乘 p1 x p2
{return Point3D(p1.y * p2.z - p1.z * p2.y,p1.z * p2.x - p1.x * p2.z,p1.x * p2.y - p1.y * p2.x );
}
double operator&(Point3D p1, Point3D p2) { // 計算點積 p1·p2return (p1.x * p2.x + p1.y * p2.y + p1.z * p2.z);
}
double Norm(Point3D p) // 計算矢量p的模
{return sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
}//取平面法向量
Point3D GetV(Area3D s){return (s.p1-s.p2)*(s.p2-s.p3);
}//判三點共線
int PointOnLine(Point3D p1,Point3D p2,Point3D p3){return ZERO( (p1-p2)*(p2-p3) );
}//判四點共面
int PointOnArea(Point3D a,Point3D b,Point3D c,Point3D d){return ZERO( GetV(Area3D(a, b, c))&(d-a) );
}//求三維空間中兩點間的距離
double Dis(Point3D p1, Point3D p2)
{return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z));
}
// 求三維空間中點到直線的距離
double Dis(Point3D p, Line3D L)
{if(L.p1==L.p2) return Dis(p, L.p1);return Norm((p - L.p1) * (L.p2 - L.p1)) / Norm(L.p2 - L.p1);
}bool OnLine(Point3D p, Line3D L) // 判斷三維空間中點p是否在直線L上
{if(L.p1==L.p2 && p==L.p1) return true;//共點時，返回truereturn ZERO( (p - L.p1) * (L.p2 - L.p1) );
}
bool OnLineSeg(Point3D p, Line3D L) // 判斷三維空間中點p是否在線段l上
{return ( ZERO((L.p1 - p) * (L.p2 - p)) &&EQ( Norm(p - L.p1) + Norm(p - L.p2), Norm(L.p2 - L.p1)) );
}// 點p到平面Ap-Al的距離。
double Dis(Point3D p, Point3D Ap, Point3D Al) {return fabs((p-Ap)&Al)/Norm(Al); // 如果不取絕對值，得到的是有向距離
}// 點p在平面Ap-Al上的投影。
Point3D PointToArea(Point3D p,Point3D Ap, Point3D Al) {Al=Al/(Norm(Al));//把Al變成法向量。return p-Al*((p-Ap)&Al);
}
//得到點p到直線L的距離，并返回p到直直線L的最近點rep
double PointToLine(Point3D p,Line3D L,Point3D& rep)
{if(L.p1==L.p2){rep=L.p1;return Dis(p,L.p1);}Point3D a,b;a = L.p2-L.p1;b = p-L.p1;double dis12 = Dis(L.p1,L.p2);double dis = ( Norm(a*b) )/dis12;double k = (a&b)/(Norm(a)*dis12) ;rep = L.p1+(L.p2-L.p1)*k;return dis;
}
//求兩條直線之間的關系（三維）
//輸入：兩條不為點的直線
//輸出：相交返回XIANGJIAO和交點p,平行返回PINGXING,共線返回GONGXIAN
int LineAndLine(Line3D L1,Line3D L2,Point3D &p)
{Point3D px,py;px = L1.p1 - L1.p2;py = L2.p1 - L2.p2;if( ZERO(px*py) )//平行或者共線
{if( ZERO( (L2.p1-L1.p1)*py ) ) //共線
{return GONGXIAN;}return PINXING;}//判斷是否共面Point3D tp=(L1.p1-L2.p1)*py;if( !ZERO(tp&px) ) return XIANGLI;//XIANGLI與平行相同
p = L1.p1;Point3D tp1=(L2.p1-L1.p1)*(L2.p1-L2.p2);Point3D tp2=(L1.p2-L1.p1)*(L2.p1-L2.p2);double _t = Norm(tp1)/Norm(tp2);//tp1和tp2肯定是共線的，如果反向則_t 為負if( LT( (tp1&tp2),0 ) ) _t*=-1;p.x += (L1.p2.x-L1.p1.x)*_t;p.y += (L1.p2.y-L1.p1.y)*_t;p.z += (L1.p2.z-L1.p1.z)*_t;return XIANGJIAO;
}//空間兩直線最近點對。直線不能平行，直線不能為點.
//ans1為直線a1,b1上的最近點
Point3D ans1,ans2;
double SegSegDistance(Point3D a1, Point3D b1, Point3D a2, Point3D b2)
{Point3D v1 = (a1-b1), v2 = (a2-b2);Point3D N = v1*v2;Point3D ab = (a1-a2);double ans = (N&ab) / Norm(N);Point3D p1 = a1, p2 = a2;Point3D d1 = b1-a1, d2 = b2-a2;double t1, t2;t1 = ((p2-p1)*d2 )&(d1*d2);t2 = ((p2-p1)*d1 )&(d1*d2);double dd = Norm( (d1*d2) );t1 /= dd*dd;t2 /= dd*dd;ans1=a1+(b1-a1)*t1;ans2=a2+(b2-a2)*t2;return fabs(ans);
}//直線與平面交
int LineAndArea(Line3D l1,Point3D Ap,Point3D Al,Point3D &p)
{//輸入直線，和平面的點法式//第一步，判斷直線與平面是否平行。if( ZERO((l1.p2-l1.p1)&Al) ) return 0;//直線與平面平行。double _t =( (Ap-l1.p1)&Al ) / ((l1.p1-l1.p2)&Al);p = l1.p1+(l1.p1-l1.p2)*_t;return 1;
}void dfs(int x,double &len)
{len++;dfs(x-1,len);dfs(x-2,len);
}//空間兩直線最近點對
//注意：直線不能平行
double LineAndLine(Line3D l1,Line3D l2,Point3D &p1,Point3D &p2)
{//先求出法向量
Point3D v1,v2;v1 = l1.p2-l1.p1;v2 = l2.p2-l2.p1;Point3D vt=v1*v2;//然后先把l2投影到 l1所在的平面上double len = ((l2.p1-l1.p1)&vt)/Norm(vt);double normvt = -len/Norm(vt);vt.x = vt.x*normvt;vt.y = vt.y*normvt;vt.z = vt.z*normvt;Line3D tl2;tl2.p1 = l2.p1+vt;tl2.p2 = l2.p2+vt;int sign=LineAndLine(l1, tl2, p1);/*//測試用if(sign!=XIANGJIAO){int x=0;printf("%lf\n",len/x);dfs(100000000,len);}*/return fabs(len);
}//已知空間四面體6條邊，求體積
double P( double a,double b,double c,double d,double e ){ return a*(b*c-d*e); }
double Get4V(int OA,int OB,int OC,int AB,int CA,int BC)
{OA*=OA;OB*=OB;OC*=OC;AB*=AB;CA*=CA;BC*=BC;double ans=0;ans+=P( OA,OB,OC,(OB+OC-BC)/2.0,(OB+OC-BC)/2.0 );ans-=P( (OA+OB-AB)/2.0,(OA+OB-AB)/2.0,OC,(OA+OC-CA)/2.0,(OB+OC-BC)/2.0 );ans+=P( (OA+OC-CA)/2.0,(OA+OB-AB)/2.0,(OB+OC-BC)/2.0,OB,(OA+OC-CA)/2.0);return sqrt(ans/36);
}//求兩面相交,平行或共面返回PINGXING,否則返回XIANGJIAO和直線rel
int AreaAndArea(Area3D a1,Area3D a2,Line3D &rel)
{Point3D va1 = GetV(a1),va2 = GetV(a2);Point3D lv = va1*va2;//相交直線的方向向量if( ZERO(lv) )//平行
{return PINXING;}//然后得到某一個交點
Point3D p1;if(LineAndArea(Line3D(a1.p1,a1.p2), a2.p1, va2, p1) == 0)if(LineAndArea(Line3D(a1.p1,a1.p3), a2.p1, va2, p1) == 0)LineAndArea(Line3D(a1.p2,a1.p3), a2.p1, va2, p1);rel.p1 = p1; rel.p2 = p1 + (lv*10);return XIANGJIAO;
}//已知3點坐標，求平面ax+by+cz+d=0;
void GetAreaABCD(Point3D p1,Point3D p2,Point3D p3,double &a,double &b,double &c,double &d)
{a = ( (p2.y-p1.y)*(p3.z-p1.z)-(p2.z-p1.z)*(p3.y-p1.y) );b = ( (p2.z-p1.z)*(p3.x-p1.x)-(p2.x-p1.x)*(p3.z-p1.z) );c = ( (p2.x-p1.x)*(p3.y-p1.y)-(p2.y-p1.y)*(p3.x-p1.x) );d = ( 0-(a*p1.x+b*p1.y+c*p1.z) );
}///*---------------------代碼區---------------------------*/
?
與50位技術專家面對面20年技術見證，附贈技術全景圖總結
- 上一篇： Java发邮件-QQ服务器认证
- 下一篇： MB与MQ简介