| 1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702 |
- /////////////////////////////////////////////
- //Copyright (c) 2011, 北京探索者软件公司
- //All rights reserved.
- //文件名称:
- //文件描述:
- //创 建 者: mjy
- //创建日期: 2011-11-02
- //版 本 号:4.0.0.0
- /////////////////////////////////////////////
- using System;
- using System.Collections.Generic;
- using System.Drawing;
- namespace SAGA.DotNetUtils.Geometry
- {
- /*
- 计算几何
- 目录
- ㈠ 点的基本运算
- 1. 平面上两点之间距离 1
- 2. 判断两点是否重合 1
- 3. 矢量叉乘 1
- 4. 矢量点乘 2
- 5. 判断点是否在线段上 2
- 6. 求一点饶某点旋转后的坐标 2
- 7. 求矢量夹角 2
- ㈡ 线段及直线的基本运算
- 1. 点与线段的关系 3
- 2. 求点到线段所在直线垂线的垂足 4
- 3. 点到线段的最近点 4
- 4. 点到线段所在直线的距离 4
- 5. 点到折线集的最近距离 4
- 6. 判断圆是否在多边形内 5
- 7. 求矢量夹角余弦 5
- 8. 求线段之间的夹角 5
- 9. 判断线段是否相交 6
- 10.判断线段是否相交但不交在端点处 6
- 11.求线段所在直线的方程 6
- 12.求直线的斜率 7
- 13.求直线的倾斜角 7
- 14.求点关于某直线的对称点 7
- 15.判断两条直线是否相交及求直线交点 7
- 16.判断线段是否相交,如果相交返回交点 7
- ㈢ 多边形常用算法模块
- 1. 判断多边形是否简单多边形 8
- 2. 检查多边形顶点的凸凹性 9
- 3. 判断多边形是否凸多边形 9
- 4. 求多边形面积 9
- 5. 判断多边形顶点的排列方向,方法一 10
- 6. 判断多边形顶点的排列方向,方法二 10
- 7. 射线法判断点是否在多边形内 10
- 8. 判断点是否在凸多边形内 11
- 9. 寻找点集的graham算法 12
- 10.寻找点集凸包的卷包裹法 13
- 11.判断线段是否在多边形内 14
- 12.求简单多边形的重心 15
- 13.求凸多边形的重心 17
- 14.求肯定在给定多边形内的一个点 17
- 15.求从多边形外一点出发到该多边形的切线 18
- 16.判断多边形的核是否存在 19
- ㈣ 圆的基本运算
- 1 .点是否在圆内 20
- 2 .求不共线的三点所确定的圆 21
- ㈤ 矩形的基本运算
- 1.已知矩形三点坐标,求第4点坐标 22
- ㈥ 常用算法的描述 22
- ㈦ 补充
- 1.两圆关系: 24
- 2.判断圆是否在矩形内: 24
- 3.点到平面的距离: 25
- 4.点是否在直线同侧: 25
- 5.镜面反射线: 25
- 6.矩形包含: 26
- 7.两圆交点: 27
- 8.两圆公共面积: 28
- 9. 圆和直线关系: 29
- 10. 内切圆: 30
- 11. 求切点: 31
- 12. 线段的左右旋: 31
- 13.公式: 32
- */
- public class Geometry2D
- {
- /* 常用的常量定义 */
- public const double INFCONST = 1E200;
- public const double EQUALSCONST = 1E-10;
- public const int MAXV = 300;
- public const double PI = Math.PI;
- /// <summary>
- /// 弧度转角度
- /// </summary>
- /// <param name="v"></param>
- /// <returns></returns>
- public static double ToAngle(double v)
- {
- return v * (180d / Math.PI);
- }
- /// <summary>
- /// 角度转弧度
- /// </summary>
- /// <param name="v"></param>
- /// <returns></returns>
- public static double ToRadian(double v)
- {
- return v * (Math.PI / 180d);
- }
- public static PointF ToPointF(PointD p)
- {
- return new PointF((float)p.X, (float)p.Y);
- }
- public static PointD ToPointF(PointF p)
- {
- return new PointD(p.X, p.Y);
- }
- /// <summary>
- /// 平面点D
- /// </summary>
- public class PointD
- {
- public double X;
- public double Y;
- public PointD()
- {
- }
- public PointD(double a = 0, double b = 0)
- {
- X = a;
- Y = b;
- }
- public bool IsEquals(PointD other)
- {
- if (other == null)
- return false;
- return Math.Abs(other.X - this.X) < 0.0000001 &&
- Math.Abs(other.Y - this.Y) < 0.0000001;
- }
- public override bool Equals(object obj)
- {
- return IsEquals(obj as PointD);
- }
- public override int GetHashCode()
- {
- return base.GetHashCode();
- }
- public override string ToString()
- {
- return string.Format("{{{0},{1}}}", X, Y);
- }
- }
- /// <summary>
- /// 线段
- /// </summary>
- public class LineSeg
- {
- public PointD Start;
- public PointD End;
- public LineSeg()
- {
- }
- public LineSeg(PointD a, PointD b)
- {
- Start = a;
- End = b;
- }
- public LineSeg(double x1, double y1, double x2, double y2)
- {
- Start = new PointD(x1, y1);
- End = new PointD(x2, y2);
- }
- }
- /// <summary>
- /// 直线的解析方程 a*x+b*y+c=0 为统一表示,约定 a >= 0
- /// </summary>
- public class LineAbc
- {
- public double A;
- public double B;
- public double C;
- public LineAbc()
- {
- }
- public LineAbc(double d1 = 1, double d2 = -1, double d3 = 0)
- {
- A = d1;
- B = d2;
- C = d3;
- }
- public double X(double y)
- {
- //a*x+b*y+c=0
- double x = ((-C) - (B * y)) / A;
- return x;
- }
- public double Y(double x)
- {
- //a*x+b*y+c=0
- double y = ((-C) - (A * x)) / B;
- return y;
- }
- }
- /**********************
- * *
- * 点的基本运算 *
- * *
- **********************/
- /// <summary>
- /// 返回两点之间欧氏距离
- /// </summary>
- /// <param name="p1"></param>
- /// <param name="p2"></param>
- /// <returns></returns>
- public static double Dist(PointD p1, PointD p2)
- {
- //return (Math.Sqrt((p1.X - p2.X) * (p1.X - p2.X) + (p1.Y - p2.Y) * (p1.Y - p2.Y)));
- return Math.Sqrt(Math.Pow(p1.X - p2.X, 2) + Math.Pow(p1.Y - p2.Y, 2));
- }
- public static double GetAngleBySinV(double sinValue)
- {
- List<double> listResult = new List<double>();
- double max = Math.PI / 2;
- for (double x = 0.001d; x < max; x = x + 0.001d)
- {
- double v = Math.Sin(x);
- if (Math.Abs(v - sinValue) < 0.001d)
- {
- listResult.Add(x);
- }
- }
- if (listResult.Count > 0)
- {
- return listResult[listResult.Count - 1];
- }
- return 0d;
- }
- public static double Dist(double x1, double y1, double x2, double y2)
- {
- return Math.Sqrt(Math.Pow(x1 - x2, 2) + Math.Pow(y1 - y2, 2));
- }
- /// <summary>
- /// 判断两个点是否重合
- /// </summary>
- /// <param name="p1"></param>
- /// <param name="p2"></param>
- /// <returns></returns>
- public static bool EqualPoint(PointD p1, PointD p2)
- {
- return ((Math.Abs(p1.X - p2.X) < EQUALSCONST) && (Math.Abs(p1.Y - p2.Y) < EQUALSCONST));
- }
- /******************************************************************************
- r=multiply(sp,ep,op),得到(sp-op)和(ep-op)的叉积
- r>0:ep在矢量opsp的逆时针方向;
- r=0:opspep三点共线;
- r<0:ep在矢量opsp的顺时针方向
- *******************************************************************************/
- public static double Multiply(PointD sp, PointD ep, PointD op)
- {
- double v = ((sp.X - op.X) * (ep.Y - op.Y) - (ep.X - op.X) * (sp.Y - op.Y));
- if (IsZero(v))
- return 0d;
- return v;
- }
- /*
- r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积,如果两个矢量都非零矢量
- r<0:两矢量夹角为锐角;
- r=0:两矢量夹角为直角;
- r>0:两矢量夹角为钝角
- *******************************************************************************/
- public static double DotMultiply(PointD p1, PointD p2, PointD p0)
- {
- double v = ((p1.X - p0.X) * (p2.X - p0.X) + (p1.Y - p0.Y) * (p2.Y - p0.Y));
- if (IsZero(v))
- return 0d;
- return v;
- }
- public static bool IsZero(double v)
- {
- return Math.Abs(v) < EQUALSCONST;
- }
- public static bool IsEqual(double v1, double v2)
- {
- return (Math.Abs(v1 - v2) < EQUALSCONST);
- }
- /******************************************************************************
- 判断点p是否在线段l上
- 条件:(p在线段l所在的直线上) && (点p在以线段l为对角线的矩形内)
- *******************************************************************************/
- public static bool Online(LineSeg l, PointD p)
- {
- return ((IsZero(Multiply(l.End, p, l.Start))) &&
- (((p.X - l.Start.X) * (p.X - l.End.X) <= 0) && ((p.Y - l.Start.Y) * (p.Y - l.End.Y) <= 0)));
- }
- /// <summary>
- /// 返回点p以点center为圆心逆时针旋转alpha(单位:弧度)后所在的位置
- /// </summary>
- /// <param name="center"></param>
- /// <param name="alpha"></param>
- /// <param name="p"></param>
- /// <returns></returns>
- public static PointD Rotate(PointD center, double alpha, PointD p)
- {
- PointD tp = new PointD();
- p.X -= center.X;
- p.Y -= center.Y;
- tp.X = p.X * Math.Cos(alpha) - p.Y * Math.Sin(alpha) + center.X;
- tp.Y = p.Y * Math.Cos(alpha) + p.X * Math.Sin(alpha) + center.Y;
- return tp;
- }
- /* 返回顶角在o点,起始边为os,终止边为oe的夹角(单位:弧度)
- 角度小于pi,返回正值
- 角度大于pi,返回负值
- 可以用于求线段之间的夹角
- 原理:
- r = dotmultiply(s,e,o) / (dist(o,s)*dist(o,e))
- r'= multiply(s,e,o)
- r >= 1 angle = 0;
- r <= -1 angle = -PI
- -1<r<1 && r'>0 angle = arccos(r)
- -1<r<1 && r'<=0 angle = -arccos(r)
- */
- public static double Angle(PointD o, PointD s, PointD e)
- {
- double cosfi, fi, norm;
- double dsx = s.X - o.X;
- double dsy = s.Y - o.Y;
- double dex = e.X - o.X;
- double dey = e.Y - o.Y;
- cosfi = dsx * dex + dsy * dey;
- norm = (dsx * dsx + dsy * dsy) * (dex * dex + dey * dey);
- cosfi /= Math.Sqrt(norm);
- if (cosfi >= 1.0) return 0;
- if (cosfi <= -1.0) return -PI;
- fi = Math.Acos(cosfi);
- // 说明矢量os 在矢量 oe的顺时针方向
- if (dsx * dey - dsy * dex > 0)
- return fi;
- return -fi;
- }
- /*****************************\
- * *
- * 线段及直线的基本运算 *
- * *
- \*****************************/
- /* 判断点与线段的关系,用途很广泛
- 本函数是根据下面的公式写的,P是点C到线段AB所在直线的垂足
- AC dot AB
- r = ---------
- ||AB||^2
- (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
- = -------------------------------
- L^2
- r has the following meaning:
- r=0 P = A
- r=1 P = B
- r<0 P is on the backward extension of AB
- r>1 P is on the forward extension of AB
- 0<r<1 P is interior to AB
- */
- public static double Relation(PointD p, LineSeg l)
- {
- LineSeg tl = new LineSeg();
- tl.Start = l.Start;
- tl.End = p;
- return DotMultiply(tl.End, l.End, l.Start) / (Dist(l.Start, l.End) * Dist(l.Start, l.End));
- }
- // 求点C到线段AB所在直线的垂足 P
- public static PointD PerpenDicular(PointD p, LineSeg l)
- {
- double r = Relation(p, l);
- PointD tp = new PointD();
- tp.X = l.Start.X + r * (l.End.X - l.Start.X);
- tp.Y = l.Start.Y + r * (l.End.Y - l.Start.Y);
- return tp;
- }
- /* 求点p到线段l的最短距离,并返回线段上距该点最近的点np
- 注意:np是线段l上到点p最近的点,不一定是垂足 */
- public static double PtoLineSegDist(PointD p, LineSeg l, ref PointD np)
- {
- double r = Relation(p, l);
- if (r < 0)
- {
- np = l.Start;
- return Dist(p, l.Start);
- }
- if (r > 1)
- {
- np = l.End;
- return Dist(p, l.End);
- }
- np = PerpenDicular(p, l);
- return Dist(p, np);
- }
- // 求点p到线段l所在直线的距离,请注意本函数与上个函数的区别
- public static double PtoLDist(PointD p, LineSeg l)
- {
- return Math.Abs(Multiply(p, l.End, l.Start)) / Dist(l.Start, l.End);
- }
- /* 计算点到折线集的最近距离,并返回最近点.
- 注意:调用的是ptolineseg()函数 */
- public static double PtoPointSet(int vcount, PointD[] pointset, PointD p, ref PointD q)
- {
- int i;
- double cd = INFCONST, td;
- LineSeg l = new LineSeg();
- PointD tq = new PointD();
- PointD cq = new PointD();
- for (i = 0; i < vcount - 1; i++)
- {
- l.Start = pointset[i];
- l.End = pointset[i + 1];
- td = PtoLineSegDist(p, l, ref tq);
- if (td < cd)
- {
- cd = td;
- cq = tq;
- }
- }
- q = cq;
- return cd;
- }
- /* 判断圆是否在多边形内.ptolineseg()函数的应用2 */
- public static bool CircleInsidePolygon(int vcount, PointD center, double radius, PointD[] polygon)
- {
- PointD q = new PointD();
- double d;
- q.X = 0;
- q.Y = 0;
- d = PtoPointSet(vcount, polygon, center, ref q);
- if (d < radius || (Math.Abs(d - radius) < EQUALSCONST))
- return true;
- else
- return false;
- }
- /* 返回两个矢量l1和l2的夹角的余弦(-1 --- 1)注意:如果想从余弦求夹角的话,注意反余弦函数的定义域是从 0到pi */
- public static double Cosine(LineSeg line1, LineSeg line2)
- {
- return (((line1.End.X - line1.Start.X) * (line2.End.X - line2.Start.X) +
- (line1.End.Y - line1.Start.Y) * (line2.End.Y - line2.Start.Y)) /
- (Dist(line1.End, line1.Start) * Dist(line2.End, line2.Start)));
- }
- /// <summary>
- /// 返回线段line1与line2之间的夹角 单位:弧度 范围(-pi,pi)
- /// </summary>
- /// <param name="line1"></param>
- /// <param name="line2"></param>
- /// <returns></returns>
- public static double LinesAngle(LineSeg line1, LineSeg line2)
- {
- PointD o = new PointD(), s = new PointD(), e = new PointD();
- o.X = o.Y = 0;
- s.X = line1.End.X - line1.Start.X;
- s.Y = line1.End.Y - line1.Start.Y;
- e.X = line2.End.X - line2.Start.X;
- e.Y = line2.End.Y - line2.Start.Y;
- return Angle(o, s, e);
- }
- // 如果线段u和v相交(包括相交在端点处)时,返回true
- //
- //判断P1P2跨立Q1Q2的依据是:( P1 - Q1 ) × ( Q2 - Q1 ) * ( Q2 - Q1 ) × ( P2 - Q1 ) >= 0。
- //判断Q1Q2跨立P1P2的依据是:( Q1 - P1 ) × ( P2 - P1 ) * ( P2 - P1 ) × ( Q2 - P1 ) >= 0。
- public static bool Intersect(LineSeg u, LineSeg v)
- {
- return ((Math.Max(u.Start.X, u.End.X) >= Math.Min(v.Start.X, v.End.X)) && //排斥实验
- (Math.Max(v.Start.X, v.End.X) >= Math.Min(u.Start.X, u.End.X)) &&
- (Math.Max(u.Start.Y, u.End.Y) >= Math.Min(v.Start.Y, v.End.Y)) &&
- (Math.Max(v.Start.Y, v.End.Y) >= Math.Min(u.Start.Y, u.End.Y)) &&
- (Multiply(v.Start, u.End, u.Start) * Multiply(u.End, v.End, u.Start) >= 0) && //跨立实验
- (Multiply(u.Start, v.End, v.Start) * Multiply(v.End, u.End, v.Start) >= 0));
- }
- // (线段u和v相交)&&(交点不是双方的端点) 时返回true
- public static bool IntersectA(LineSeg u, LineSeg v)
- {
- return ((Intersect(u, v)) &&
- (!Online(u, v.Start)) &&
- (!Online(u, v.End)) &&
- (!Online(v, u.End)) &&
- (!Online(v, u.Start)));
- }
- // 线段v所在直线与线段u相交时返回true;方法:判断线段u是否跨立线段v
- public static bool Intersectl(LineSeg u, LineSeg v)
- {
- return Multiply(u.Start, v.End, v.Start) * Multiply(v.End, u.End, v.Start) >= 0;
- }
- // 根据已知两点坐标,求过这两点的直线解析方程: a*x+b*y+c = 0 (a >= 0)
- public static LineAbc MakeLine(PointD p1, PointD p2)
- {
- LineAbc tl = new LineAbc();
- int sign = 1;
- tl.A = p2.Y - p1.Y;
- if (tl.A < 0)
- {
- sign = -1;
- tl.A = sign * tl.A;
- }
- tl.B = sign * (p1.X - p2.X);
- tl.C = sign * (p1.Y * p2.X - p1.X * p2.Y);
- return tl;
- }
- // 根据直线解析方程返回直线的斜率k,水平线返回 0,竖直线返回 1e200
- public static double Slope(LineAbc line)
- {
- if (Math.Abs(line.A) < 1e-20)
- return 0;
- if (Math.Abs(line.B) < 1e-20)
- return INFCONST;
- return -(line.A / line.B);
- }
- // 返回直线的倾斜角alpha ( 0 - pi)
- public static double Alpha(LineAbc line)
- {
- if (Math.Abs(line.A) < EQUALSCONST)
- return 0;
- if (Math.Abs(line.B) < EQUALSCONST)
- return PI / 2;
- double k = Slope(line);
- if (k > 0)
- return Math.Atan(k);
- else
- return PI + Math.Atan(k);
- }
- // 求点p关于直线l的对称点
- public static PointD Symmetry(LineAbc line, PointD p)
- {
- PointD tp = new PointD();
- tp.X = ((line.B * line.B - line.A * line.A) * p.X - 2 * line.A * line.B * p.Y - 2 * line.A * line.C) /
- (line.A * line.A + line.B * line.B);
- tp.Y = ((line.A * line.A - line.B * line.B) * p.Y - 2 * line.A * line.B * p.X - 2 * line.B * line.C) /
- (line.A * line.A + line.B * line.B);
- return tp;
- }
- // 如果两条直线 l1(a1*x+b1*y+c1 = 0), l2(a2*x+b2*y+c2 = 0)相交,返回true,且返回交点p
- public static bool LineIntersect(LineAbc line1, LineAbc line2, ref PointD p) // 是 L1,L2
- {
- double d = line1.A * line2.B - line2.A * line1.B;
- if (Math.Abs(d) < EQUALSCONST) // 不相交
- return false;
- p.X = (line2.C * line1.B - line1.C * line2.B) / d;
- p.Y = (line2.A * line1.C - line1.A * line2.C) / d;
- return true;
- }
- // 如果线段l1和l2相交,返回true且交点由(inter)返回,否则返回false
- public static bool Intersection(LineSeg line1, LineSeg line2, ref PointD p)
- {
- LineAbc ll1 = new LineAbc(), ll2 = new LineAbc();
- ll1 = MakeLine(line1.Start, line1.End);
- ll2 = MakeLine(line2.Start, line2.End);
- if (LineIntersect(ll1, ll2, ref p))
- return Online(line1, p);
- else
- return false;
- }
- /******************************\
- * *
- * 多边形常用算法模块 *
- * *
- \******************************/
- // 如果无特别说明,输入多边形顶点要求按逆时针排列
- /*
- 返回值:输入的多边形是简单多边形,返回true
- 要 求:输入顶点序列按逆时针排序
- 说 明:简单多边形定义:
- 1:循环排序中相邻线段对的交是他们之间共有的单个点
- 2:不相邻的线段不相交
- 本程序默认第一个条件已经满足
- */
- public static bool IsSimple(int vcount, PointD[] polygon)
- {
- int i, cn;
- LineSeg line1 = new LineSeg(), line2 = new LineSeg();
- for (i = 0; i < vcount; i++)
- {
- line1.Start = polygon[i];
- line1.End = polygon[(i + 1) % vcount];
- cn = vcount - 3;
- while (cn != 0)
- {
- line2.Start = polygon[(i + 2) % vcount];
- line2.End = polygon[(i + 3) % vcount];
- if (Intersect(line1, line2))
- break;
- cn--;
- }
- if (cn != 0)
- return false;
- }
- return true;
- }
- /// <summary>
- /// 返回值:按输入顺序返回多边形顶点的凸凹性判断,bc[i]=1,iff:第i个顶点是凸顶点
- /// </summary>
- public static void CheckConvex(int vcount, PointD[] polygon, bool[] bc)
- {
- int i, index = 0;
- PointD tp = polygon[0];
- for (i = 1; i < vcount; i++) // 寻找第一个凸顶点
- {
- if (polygon[i].Y < tp.Y || (polygon[i].Y == tp.Y && polygon[i].X < tp.X))
- {
- tp = polygon[i];
- index = i;
- }
- }
- int count = vcount - 1;
- bc[index] = true;
- while (count != 0) // 判断凸凹性
- {
- if (Multiply(polygon[(index + 1) % vcount], polygon[(index + 2) % vcount], polygon[index]) >= 0)
- bc[(index + 1) % vcount] = true;
- else
- bc[(index + 1) % vcount] = false;
- index++;
- count--;
- }
- }
- /// <summary>
- /// 返回值:多边形polygon是凸多边形时,返回true
- /// </summary>
- public static bool IsConvex(int vcount, PointD[] polygon)
- {
- bool[] bc = new bool[MAXV];
- CheckConvex(vcount, polygon, bc);
- for (int i = 0; i < vcount; i++) // 逐一检查顶点,是否全部是凸顶点
- if (!bc[i])
- return false;
- return true;
- }
- /// <summary>
- /// 返回多边形面积(signed);输入顶点按逆时针排列时,返回正值;否则返回负值
- /// </summary>
- public static double AreaOfPolygon(int vcount, PointD[] polygon)
- {
- int i;
- double s;
- if (vcount < 3)
- return 0;
- s = polygon[0].Y * (polygon[vcount - 1].X - polygon[1].X);
- for (i = 1; i < vcount; i++)
- s += polygon[i].Y * (polygon[(i - 1)].X - polygon[(i + 1) % vcount].X);
- return s / 2;
- }
- /// <summary>
- /// 如果输入顶点按逆时针排列,返回true
- /// </summary>
- public static bool IsConterClock(int vcount, PointD[] polygon)
- {
- return AreaOfPolygon(vcount, polygon) > 0;
- }
- /// <summary>
- /// 另一种判断多边形顶点排列方向的方法
- /// </summary>
- public static bool IsCcwize(int vcount, PointD[] polygon)
- {
- int i, index;
- PointD a, b, v;
- v = polygon[0];
- index = 0;
- for (i = 1; i < vcount; i++) // 找到最低且最左顶点,肯定是凸顶点
- {
- if (polygon[i].Y < v.Y || polygon[i].Y == v.Y && polygon[i].X < v.X)
- {
- index = i;
- }
- }
- a = polygon[(index - 1 + vcount) % vcount]; // 顶点v的前一顶点
- b = polygon[(index + 1) % vcount]; // 顶点v的后一顶点
- return Multiply(v, b, a) > 0;
- }
- /********************************************************************************************
- 射线法判断点q与多边形polygon的位置关系,要求polygon为简单多边形,顶点逆时针排列
- 如果点在多边形内: 返回0
- 如果点在多边形边上: 返回1
- 如果点在多边形外: 返回2
- *********************************************************************************************/
- public static int InSidePolygon(PointD[] polygon, PointD q)
- {
- int c = 0, i, n;
- LineSeg l1 = new LineSeg(), l2 = new LineSeg();
- l1.Start = q;
- l1.End = new PointD(INFCONST, q.Y);
- n = polygon.Length;
- PointD originStart = null;
- for (i = 0; i < n; i++)
- {
- var start = polygon[i];
- var end = polygon[(i + 1) % n];
- var end2 = polygon[(i + 2) % n];
- var end3 = polygon[(i + 3) % n];
- l2.Start = start;
- l2.End = end;
- if (Online(l2, q))
- return 1; // 如果点在边上,返回1
- if (IntersectA(l1, l2)) // 相交且不在端点
- {
- c++;
- }
- else if(Online(l1, end))// end点在射线上
- {
- //判断下一个点end2是否在射线上
- if (Online(l1, end2))
- {
- /* 下一条边是水平线,前一个端点和后一个端点在射线两侧 */
- originStart = originStart ?? start;
- var multiplyresult = Multiply(originStart, end2, l1.Start) *
- Multiply(end2, end3, l1.Start);
- if (multiplyresult > 0)
- {
- c++;
- originStart = null;
- }
- }
- else
- {
- /* 前一个端点和后一个端点在射线两侧 */
- if (Multiply(start, end, l1.Start) *
- Multiply(end, end2, l1.Start) > 0)
- c++;
- }
- }
- }
- if (c % 2 == 1)
- return 0;
- else
- return 2;
- }
- //点q是凸多边形polygon内时,返回true;注意:多边形polygon一定要是凸多边形
- public static bool InsideConvexPolygon(int vcount, PointD[] polygon, PointD q) // 可用于三角形!
- {
- PointD p = new PointD();
- LineSeg l = new LineSeg();
- int i;
- p.X = 0;
- p.Y = 0;
- for (i = 0; i < vcount; i++) // 寻找一个肯定在多边形polygon内的点p:多边形顶点平均值
- {
- p.X += polygon[i].X;
- p.Y += polygon[i].Y;
- }
- p.X /= vcount;
- p.Y /= vcount;
- for (i = 0; i < vcount; i++)
- {
- l.Start = polygon[i];
- l.End = polygon[(i + 1) % vcount];
- if (Multiply(p, l.End, l.Start) * Multiply(q, l.End, l.Start) < 0) /* 点p和点q在边l的两侧,说明点q肯定在多边形外 */
- break;
- }
- return (i == vcount);
- }
- /**********************************************
- 寻找凸包的graham 扫描法
- PointSet为输入的点集;
- ch为输出的凸包上的点集,按照逆时针方向排列;
- n为PointSet中的点的数目
- len为输出的凸包上的点的个数
- **********************************************/
- public static void GrahamScan(PointD[] PointSet, PointD[] ch, int n, ref int len)
- {
- int i, j, k = 0, top = 2;
- PointD tmp = new PointD();
- // 选取PointSet中y坐标最小的点PointSet[k],如果这样的点有多个,则取最左边的一个
- for (i = 1; i < n; i++)
- if (PointSet[i].Y < PointSet[k].Y || (PointSet[i].Y == PointSet[k].Y) && (PointSet[i].X < PointSet[k].X))
- k = i;
- tmp = PointSet[0];
- PointSet[0] = PointSet[k];
- PointSet[k] = tmp; // 现在PointSet中y坐标最小的点在PointSet[0]
- for (i = 1; i < n - 1; i++) /* 对顶点按照相对PointSet[0]的极角从小到大进行排序,极角相同的按照距离PointSet[0]从近到远进行排序 */
- {
- k = i;
- for (j = i + 1; j < n; j++)
- if (Multiply(PointSet[j], PointSet[k], PointSet[0]) > 0 || // 极角更小
- (Multiply(PointSet[j], PointSet[k], PointSet[0]) == 0) && /* 极角相等,距离更短 */
- Dist(PointSet[0], PointSet[j]) < Dist(PointSet[0], PointSet[k])
- )
- k = j;
- tmp = PointSet[i];
- PointSet[i] = PointSet[k];
- PointSet[k] = tmp;
- }
- ch[0] = PointSet[0];
- ch[1] = PointSet[1];
- ch[2] = PointSet[2];
- for (i = 3; i < n; i++)
- {
- while (Multiply(PointSet[i], ch[top], ch[top - 1]) >= 0)
- top--;
- ch[++top] = PointSet[i];
- }
- len = top + 1;
- }
- // 卷包裹法求点集凸壳,参数说明同graham算法
- public static void ConvexClosure(PointD[] PointSet, PointD[] ch, int n, ref int len)
- {
- int top = 0, i, index, first;
- double curmax, curcos, curdis;
- PointD tmp = new PointD();
- LineSeg l1 = new LineSeg(), l2 = new LineSeg();
- bool[] use = new bool[MAXV];
- tmp = PointSet[0];
- index = 0;
- // 选取y最小点,如果多于一个,则选取最左点
- for (i = 1; i < n; i++)
- {
- if (PointSet[i].Y < tmp.Y || PointSet[i].Y == tmp.Y && PointSet[i].X < tmp.X)
- {
- index = i;
- }
- use[i] = false;
- }
- tmp = PointSet[index];
- first = index;
- use[index] = true;
- index = -1;
- ch[top++] = tmp;
- tmp.X -= 100;
- l1.Start = tmp;
- l1.End = ch[0];
- l2.Start = ch[0];
- while (index != first)
- {
- curmax = -100;
- curdis = 0;
- // 选取与最后一条确定边夹角最小的点,即余弦值最大者
- for (i = 0; i < n; i++)
- {
- if (use[i]) continue;
- l2.End = PointSet[i];
- curcos = Cosine(l1, l2); // 根据cos值求夹角余弦,范围在 (-1 -- 1 )
- if (curcos > curmax || Math.Abs(curcos - curmax) < 1e-6 && Dist(l2.Start, l2.End) > curdis)
- {
- curmax = curcos;
- index = i;
- curdis = Dist(l2.Start, l2.End);
- }
- }
- use[first] = false; //清空第first个顶点标志,使最后能形成封闭的hull
- use[index] = true;
- ch[top++] = PointSet[index];
- l1.Start = ch[top - 2];
- l1.End = ch[top - 1];
- l2.Start = ch[top - 1];
- }
- len = top - 1;
- }
- /*********************************************************************************************
- 判断线段是否在简单多边形内(注意:如果多边形是凸多边形,下面的算法可以化简)
- 必要条件一:线段的两个端点都在多边形内;
- 必要条件二:线段和多边形的所有边都不内交;
- 用途: 1. 判断折线是否在简单多边形内
- 2. 判断简单多边形是否在另一个简单多边形内
- **********************************************************************************************/
- public static bool LinesegInsidePolygon(int vcount, PointD[] polygon, LineSeg l)
- {
- // 判断线端l的端点是否不都在多边形内
- //if(!insidepolygon(vcount,polygon,l.s)||!insidepolygon(vcount,polygon,l.e))
- if (!(InSidePolygon( polygon, l.Start) != 0) || !(InSidePolygon(polygon, l.End) != 0))
- return false;
- int top = 0, i, j;
- PointD[] PointSet = new PointD[MAXV];
- PointD tmp = new PointD();
- LineSeg s = new LineSeg();
- for (i = 0; i < vcount; i++)
- {
- s.Start = polygon[i];
- s.End = polygon[(i + 1) % vcount];
- if (Online(s, l.Start)) //线段l的起始端点在线段s上
- PointSet[top++] = l.Start;
- else if (Online(s, l.End)) //线段l的终止端点在线段s上
- PointSet[top++] = l.End;
- else
- {
- if (Online(l, s.Start)) //线段s的起始端点在线段l上
- PointSet[top++] = s.Start;
- else if (Online(l, s.End)) // 线段s的终止端点在线段l上
- PointSet[top++] = s.End;
- else
- {
- if (Intersect(l, s)) // 这个时候如果相交,肯定是内交,返回false
- return false;
- }
- }
- }
- for (i = 0; i < top - 1; i++) /* 冒泡排序,x坐标小的排在前面;x坐标相同者,y坐标小的排在前面 */
- {
- for (j = i + 1; j < top; j++)
- {
- if (PointSet[i].X > PointSet[j].X ||
- Math.Abs(PointSet[i].X - PointSet[j].X) < EQUALSCONST && PointSet[i].Y > PointSet[j].Y)
- {
- tmp = PointSet[i];
- PointSet[i] = PointSet[j];
- PointSet[j] = tmp;
- }
- }
- }
- for (i = 0; i < top - 1; i++)
- {
- tmp.X = (PointSet[i].X + PointSet[i + 1].X) / 2; //得到两个相邻交点的中点
- tmp.Y = (PointSet[i].Y + PointSet[i + 1].Y) / 2;
- if (!(InSidePolygon( polygon, tmp) != 0))
- return false;
- }
- return true;
- }
- /*********************************************************************************************
- 求任意简单多边形polygon的重心
- 需要调用下面几个函数:
- void AddPosPart(); 增加右边区域的面积
- void AddNegPart(); 增加左边区域的面积
- void AddRegion(); 增加区域面积
- 在使用该程序时,如果把xtr,ytr,wtr,xtl,ytl,wtl设成全局变量就可以使这些函数的形式得到化简,
- 但要注意函数的声明和调用要做相应变化
- **********************************************************************************************/
- public static void AddPosPart(double x, double y, double w, ref double xtr, ref double ytr, ref double wtr)
- {
- if (Math.Abs(wtr + w) < 1e-10) return; // detect zero regions
- xtr = (wtr * xtr + w * x) / (wtr + w);
- ytr = (wtr * ytr + w * y) / (wtr + w);
- wtr = w + wtr;
- return;
- }
- public static void AddNegPart(double x, double y, double w, ref double xtl, ref double ytl, ref double wtl)
- {
- if (Math.Abs(wtl + w) < 1e-10)
- return; // detect zero regions
- xtl = (wtl * xtl + w * x) / (wtl + w);
- ytl = (wtl * ytl + w * y) / (wtl + w);
- wtl = w + wtl;
- return;
- }
- public static void AddRegion(double x1, double y1, double x2, double y2, ref double xtr, ref double ytr,
- ref double wtr, ref double xtl, ref double ytl, ref double wtl)
- {
- if (Math.Abs(x1 - x2) < 1e-10)
- return;
- if (x2 > x1)
- {
- AddPosPart((x2 + x1) / 2, y1 / 2, (x2 - x1) * y1, ref xtr, ref ytr, ref wtr); /* rectangle 全局变量变化处 */
- AddPosPart((x1 + x2 + x2) / 3, (y1 + y1 + y2) / 3, (x2 - x1) * (y2 - y1) / 2, ref xtr, ref ytr, ref wtr);
- // triangle 全局变量变化处
- }
- else
- {
- AddNegPart((x2 + x1) / 2, y1 / 2, (x2 - x1) * y1, ref xtl, ref ytl, ref wtl);
- // rectangle 全局变量变化处
- AddNegPart((x1 + x2 + x2) / 3, (y1 + y1 + y2) / 3, (x2 - x1) * (y2 - y1) / 2, ref xtl, ref ytl, ref wtl);
- // triangle 全局变量变化处
- }
- }
- public static PointD CgSimple(int vcount, PointD[] polygon)
- {
- double xtr, ytr, wtr, xtl, ytl, wtl;
- //注意: 如果把xtr,ytr,wtr,xtl,ytl,wtl改成全局变量后这里要删去
- PointD p1 = new PointD(), p2 = new PointD(), tp = new PointD();
- xtr = ytr = wtr = 0.0;
- xtl = ytl = wtl = 0.0;
- for (int i = 0; i < vcount; i++)
- {
- p1 = polygon[i];
- p2 = polygon[(i + 1) % vcount];
- AddRegion(p1.X, p1.Y, p2.X, p2.Y, ref xtr, ref ytr, ref wtr, ref xtl, ref ytl, ref wtl); //全局变量变化处
- }
- tp.X = (wtr * xtr + wtl * xtl) / (wtr + wtl);
- tp.Y = (wtr * ytr + wtl * ytl) / (wtr + wtl);
- return tp;
- }
- // 求凸多边形的重心,要求输入多边形按逆时针排序
- public static PointD GravityCenter(int vcount, PointD[] polygon)
- {
- PointD tp = new PointD();
- double x, y, s, x0, y0, cs, k;
- x = 0;
- y = 0;
- s = 0;
- for (int i = 1; i < vcount - 1; i++)
- {
- x0 = (polygon[0].X + polygon[i].X + polygon[i + 1].X) / 3;
- y0 = (polygon[0].Y + polygon[i].Y + polygon[i + 1].Y) / 3; //求当前三角形的重心
- cs = Multiply(polygon[i], polygon[i + 1], polygon[0]) / 2;
- //三角形面积可以直接利用该公式求解
- if (Math.Abs(s) < 1e-20)
- {
- x = x0;
- y = y0;
- s += cs;
- continue;
- }
- k = cs / s; //求面积比例
- x = (x + k * x0) / (1 + k);
- y = (y + k * y0) / (1 + k);
- s += cs;
- }
- tp.X = x;
- tp.Y = y;
- return tp;
- }
- /************************************************
- 给定一简单多边形,找出一个肯定在该多边形内的点
- 定理1 :每个多边形至少有一个凸顶点
- 定理2 :顶点数>=4的简单多边形至少有一条对角线
- 结论 : x坐标最大,最小的点肯定是凸顶点
- y坐标最大,最小的点肯定是凸顶点
- ************************************************/
- public static PointD APointInsidePoly(int vcount, PointD[] polygon)
- {
- PointD v = new PointD(), a = new PointD(), b = new PointD(), r = new PointD();
- int i, index;
- v = polygon[0];
- index = 0;
- for (i = 1; i < vcount; i++) //寻找一个凸顶点
- {
- if (polygon[i].Y < v.Y)
- {
- v = polygon[i];
- index = i;
- }
- }
- a = polygon[(index - 1 + vcount) % vcount]; //得到v的前一个顶点
- b = polygon[(index + 1) % vcount]; //得到v的后一个顶点
- PointD[] tri = new PointD[3];
- PointD q = new PointD();
- tri[0] = a;
- tri[1] = v;
- tri[2] = b;
- double md = INFCONST;
- int in1 = index;
- bool bin = false;
- for (i = 0; i < vcount; i++) //寻找在三角形avb内且离顶点v最近的顶点q
- {
- if (i == index) continue;
- if (i == (index - 1 + vcount) % vcount) continue;
- if (i == (index + 1) % vcount) continue;
- if (!InsideConvexPolygon(3, tri, polygon[i])) continue;
- bin = true;
- if (Dist(v, polygon[i]) < md)
- {
- q = polygon[i];
- md = Dist(v, q);
- }
- }
- if (!bin) //没有顶点在三角形avb内,返回线段ab中点
- {
- r.X = (a.X + b.X) / 2;
- r.Y = (a.Y + b.Y) / 2;
- return r;
- }
- r.X = (v.X + q.X) / 2; //返回线段vq的中点
- r.Y = (v.Y + q.Y) / 2;
- return r;
- }
- /***********************************************************************************************
- 求从多边形外一点p出发到一个简单多边形的切线,如果存在返回切点,其中rp点是右切点,lp是左切点
- 注意:p点一定要在多边形外 ,输入顶点序列是逆时针排列
- 原 理: 如果点在多边形内肯定无切线;凸多边形有唯一的两个切点,凹多边形就可能有多于两个的切点;
- 如果polygon是凸多边形,切点只有两个只要找到就可以,可以化简此算法
- 如果是凹多边形还有一种算法可以求解:先求凹多边形的凸包,然后求凸包的切线
- /***********************************************************************************************/
- public static void PointTangentPoly(int vcount, PointD[] polygon, PointD p, ref PointD rp, ref PointD lp)
- {
- LineSeg ep = new LineSeg(), en = new LineSeg();
- bool blp, bln;
- rp = polygon[0];
- lp = polygon[0];
- for (int i = 1; i < vcount; i++)
- {
- ep.Start = polygon[(i + vcount - 1) % vcount];
- ep.End = polygon[i];
- en.Start = polygon[i];
- en.End = polygon[(i + 1) % vcount];
- blp = Multiply(ep.End, p, ep.Start) >= 0; // p is to the left of pre edge
- bln = Multiply(en.End, p, en.Start) >= 0; // p is to the left of next edge
- if (!blp && bln)
- {
- if (Multiply(polygon[i], rp, p) > 0) // polygon[i] is above rp
- rp = polygon[i];
- }
- if (blp && !bln)
- {
- if (Multiply(lp, polygon[i], p) > 0) // polygon[i] is below lp
- lp = polygon[i];
- }
- }
- return;
- }
- // 如果多边形polygon的核存在,返回true,返回核上的一点p.顶点按逆时针方向输入
- public static bool CoreExist(int vcount, PointD[] polygon, ref PointD p)
- {
- int i, j, k;
- LineSeg l = new LineSeg();
- LineAbc[] lineset = new LineAbc[MAXV];
- for (i = 0; i < vcount; i++)
- {
- lineset[i] = MakeLine(polygon[i], polygon[(i + 1) % vcount]);
- }
- for (i = 0; i < vcount; i++)
- {
- for (j = 0; j < vcount; j++)
- {
- if (i == j) continue;
- if (LineIntersect(lineset[i], lineset[j], ref p))
- {
- for (k = 0; k < vcount; k++)
- {
- l.Start = polygon[k];
- l.End = polygon[(k + 1) % vcount];
- if (Multiply(p, l.End, l.Start) > 0)
- //多边形顶点按逆时针方向排列,核肯定在每条边的左侧或边上
- break;
- }
- if (k == vcount) //找到了一个核上的点
- break;
- }
- }
- if (j < vcount) break;
- }
- if (i < vcount)
- return true;
- else
- return false;
- }
- /*************************\
- * *
- * 圆的基本运算 *
- * *
- \*************************/
- /******************************************************************************
- 返回值 : 点p在圆内(包括边界)时,返回true
- 用途 : 因为圆为凸集,所以判断点集,折线,多边形是否在圆内时,
- 只需要逐一判断点是否在圆内即可。
- *******************************************************************************/
- public static bool PointInCircle(PointD o, double r, PointD p)
- {
- double d2 = (p.X - o.X) * (p.X - o.X) + (p.Y - o.Y) * (p.Y - o.Y);
- double r2 = r * r;
- return d2 < r2 || Math.Abs(d2 - r2) < EQUALSCONST;
- }
- /******************************************************************************
- 用 途 :求不共线的三点确定一个圆
- 输 入 :三个点p1,p2,p3
- 返回值 :如果三点共线,返回false;反之,返回true。圆心由c返回,半径由r返回
- *******************************************************************************/
- public static bool GetCircle(PointD p1, PointD p2, PointD p3, ref PointD c, ref double r)
- {
- double x12 = p2.X - p1.X;
- double y12 = p2.Y - p1.Y;
- double x13 = p3.X - p1.X;
- double y13 = p3.Y - p1.Y;
- double z2 = x12 * (p1.X + p2.X) + y12 * (p1.Y + p2.Y);
- double z3 = x13 * (p1.X + p3.X) + y13 * (p1.Y + p3.Y);
- double d = 2.0 * (x12 * (p3.Y - p2.Y) - y12 * (p3.X - p2.X));
- if (Math.Abs(d) < EQUALSCONST) //共线,圆不存在
- return false;
- c.X = (y13 * z2 - y12 * z3) / d;
- c.Y = (x12 * z3 - x13 * z2) / d;
- r = Dist(p1, c);
- return true;
- }
- public static int LineCircle(LineAbc l, PointD o, double r, ref PointD p1, ref PointD p2)
- {
- return 1;
- }
- /**************************\
- * *
- * 矩形的基本运算 *
- * *
- \**************************/
- /*
- 说明:因为矩形的特殊性,常用算法可以化简:
- 1.判断矩形是否包含点
- 只要判断该点的横坐标和纵坐标是否夹在矩形的左右边和上下边之间。
- 2.判断线段、折线、多边形是否在矩形中
- 因为矩形是个凸集,所以只要判断所有端点是否都在矩形中就可以了。
- 3.判断圆是否在矩形中
- 圆在矩形中的充要条件是:圆心在矩形中且圆的半径小于等于圆心到矩形四边的距离的最小值。
- */
- // 已知矩形的三个顶点(a,b,c),计算第四个顶点d的坐标. 注意:已知的三个顶点可以是无序的
- public static PointD Rect4th(PointD a, PointD b, PointD c)
- {
- PointD d = new PointD();
- if (Math.Abs(DotMultiply(a, b, c)) < EQUALSCONST) // 说明c点是直角拐角处
- {
- d.X = a.X + b.X - c.X;
- d.Y = a.Y + b.Y - c.Y;
- }
- if (Math.Abs(DotMultiply(a, c, b)) < EQUALSCONST) // 说明b点是直角拐角处
- {
- d.X = a.X + c.X - b.X;
- d.Y = a.Y + c.Y - b.X;
- }
- if (Math.Abs(DotMultiply(c, b, a)) < EQUALSCONST) // 说明a点是直角拐角处
- {
- d.X = c.X + b.X - a.X;
- d.Y = c.Y + b.Y - a.Y;
- }
- return d;
- }
- /*************************\
- * *
- * 常用算法的描述 *
- * *
- \*************************/
- /*
- 尚未实现的算法:
- 1. 求包含点集的最小圆
- 2. 求多边形的交
- 3. 简单多边形的三角剖分
- 4. 寻找包含点集的最小矩形
- 5. 折线的化简
- 6. 判断矩形是否在矩形中
- 7. 判断矩形能否放在矩形中
- 8. 矩形并的面积与周长
- 9. 矩形并的轮廓
- 10.矩形并的闭包
- 11.矩形的交
- 12.点集中的最近点对
- 13.多边形的并
- 14.圆的交与并
- 15.直线与圆的关系
- 16.线段与圆的关系
- 17.求多边形的核监视摄象机
- 18.求点集中不相交点对 railwai
- */
- /*
- 寻找包含点集的最小矩形
- 原理:该矩形至少一条边与点集的凸壳的某条边共线
- First take the convex hull of the points. Let the resulting convex
- polygon be P. It has been known for some time that the minimum
- area rectangle enclosing P must have one rectangle side flush with
- (i.e., collinear with and overlapping) one edge of P. This geometric
- fact was used by Godfried Toussaint to develop the "rotating calipers"
- algorithm in a hard-to-find 1983 paper, "Solving Geometric Problems
- with the Rotating Calipers" (Proc. IEEE MELECON). The algorithm
- rotates a surrounding rectangle from one flush edge to the next,
- keeping track of the minimum area for each edge. It achieves O(n)
- time (after hull computation). See the "Rotating Calipers Homepage"
- http://www.cs.mcgill.ca/~orm/rotcal.frame.html for a description
- and applet.
- */
- /*
- 折线的化简 伪码如下:
- Input: tol = the approximation tolerance
- L = {V0,V1,...,Vn-1} is any n-vertex polyline
- Set start = 0;
- Set k = 0;
- Set W0 = V0;
- for each vertex Vi (i=1,n-1)
- {
- if Vi is within tol from Vstart
- then ignore it, and continue with the next vertex
- Vi is further than tol away from Vstart
- so add it as a new vertex of the reduced polyline
- Increment k++;
- Set Wk = Vi;
- Set start = i; as the new initial vertex
- }
- Output: W = {W0,W1,...,Wk-1} = the k-vertex simplified polyline
- */
- /********************\
- * *
- * 补充 *
- * *
- \********************/
- //两圆关系:
- /* 两圆:
- 相离: return 1;
- 外切: return 2;
- 相交: return 3;
- 内切: return 4;
- 内含: return 5;
- */
- public static int CircleRelation(PointD p1, double r1, PointD p2, double r2)
- {
- double d = Math.Sqrt((p1.X - p2.X) * (p1.X - p2.X) + (p1.Y - p2.Y) * (p1.Y - p2.Y));
- if (Math.Abs(d - r1 - r2) < EQUALSCONST) // 必须保证前两个if先被判定!
- return 2;
- if (Math.Abs(d - Math.Abs(r1 - r2)) < EQUALSCONST)
- return 4;
- if (d > r1 + r2)
- return 1;
- if (d < Math.Abs(r1 - r2))
- return 5;
- if (Math.Abs(r1 - r2) < d && d < r1 + r2)
- return 3;
- return 0; // indicate an error!
- }
- //判断圆是否在矩形内:
- // 判定圆是否在矩形内,是就返回true(设矩形水平,且其四个顶点由左上开始按顺时针排列)
- // 调用ptoldist函数,在第4页
- public static bool CircleRecRelation(PointD pc, double r, PointD pr1, PointD pr2, PointD pr3, PointD pr4)
- {
- if (pr1.X < pc.X && pc.X < pr2.X && pr3.Y < pc.Y && pc.Y < pr2.Y)
- {
- LineSeg line1 = new LineSeg(pr1, pr2);
- LineSeg line2 = new LineSeg(pr2, pr3);
- LineSeg line3 = new LineSeg(pr3, pr4);
- LineSeg line4 = new LineSeg(pr4, pr1);
- if (r < PtoLDist(pc, line1) && r < PtoLDist(pc, line2) && r < PtoLDist(pc, line3) &&
- r < PtoLDist(pc, line4))
- return true;
- }
- return false;
- }
- //点到平面的距离:
- //点到平面的距离,平面用一般式表示ax+by+cz+d=0
- public static double P2planeDist(double x, double y, double z, double a, double b, double c, double d)
- {
- return Math.Abs(a * x + b * y + c * z + d) / Math.Sqrt(a * a + b * b + c * c);
- }
- //点是否在直线同侧:
- //两个点是否在直线同侧,是则返回true
- public static bool SameSide(PointD p1, PointD p2, LineAbc line)
- {
- return (line.A * p1.X + line.B * p1.Y + line.C) * (line.A * p2.X + line.B * p2.Y + line.C) > 0;
- }
- //镜面反射线:
- // 已知入射线、镜面,求反射线。
- // a1,b1,c1为镜面直线方程(a1 x + b1 y + c1 = 0 ,下同)系数;
- //a2,b2,c2为入射光直线方程系数;
- //a,b,c为反射光直线方程系数.
- // 光是有方向的,使用时注意:入射光向量:<-b2,a2>;反射光向量:<b,-a>.
- // 不要忘记结果中可能会有"negative zeros"
- public static void Reflect(double a1, double b1, double c1, double a2, double b2, double c2, ref double a,
- ref double b, ref double c)
- {
- double n, m;
- double tpb, tpa;
- tpb = b1 * b2 + a1 * a2;
- tpa = a2 * b1 - a1 * b2;
- m = (tpb * b1 + tpa * a1) / (b1 * b1 + a1 * a1);
- n = (tpa * b1 - tpb * a1) / (b1 * b1 + a1 * a1);
- if (Math.Abs(a1 * b2 - a2 * b1) < 1e-20)
- {
- a = a2;
- b = b2;
- c = c2;
- return;
- }
- double xx, yy; //(xx,yy)是入射线与镜面的交点。
- xx = (b1 * c2 - b2 * c1) / (a1 * b2 - a2 * b1);
- yy = (a2 * c1 - a1 * c2) / (a1 * b2 - a2 * b1);
- a = n;
- b = -m;
- c = m * yy - xx * n;
- }
- //矩形包含:
- // 矩形2(C,D)是否在1(A,B)内
- public static bool R2InrL(double A, double B, double C, double D)
- {
- double X, Y, L, K, DMax;
- if (A < B)
- {
- double tmp = A;
- A = B;
- B = tmp;
- }
- if (C < D)
- {
- double tmp = C;
- C = D;
- D = tmp;
- }
- if (A > C && B > D) // trivial case
- return true;
- else if (D >= B)
- return false;
- else
- {
- X = Math.Sqrt(A * A + B * B); // outer rectangle's diagonal
- Y = Math.Sqrt(C * C + D * D); // inner rectangle's diagonal
- if (Y < B) // check for marginal conditions
- return true; // the inner rectangle can freely rotate inside
- else if (Y > X)
- return false;
- else
- {
- L = (B - Math.Sqrt(Y * Y - A * A)) / 2;
- K = (A - Math.Sqrt(Y * Y - B * B)) / 2;
- DMax = Math.Sqrt(L * L + K * K);
- if (D >= DMax)
- return false;
- else
- return true;
- }
- }
- }
- //两圆交点:
- // 两圆已经相交(相切)
- public static void C2Point(PointD p1, double r1, PointD p2, double r2, ref PointD rp1, ref PointD rp2)
- {
- double a, b, r;
- a = p2.X - p1.X;
- b = p2.Y - p1.Y;
- r = (a * a + b * b + r1 * r1 - r2 * r2) / 2;
- if (a == 0 && b != 0)
- {
- rp1.Y = rp2.Y = r / b;
- rp1.X = Math.Sqrt(r1 * r1 - rp1.Y * rp1.Y);
- rp2.X = -rp1.X;
- }
- else if (a != 0 && b == 0)
- {
- rp1.X = rp2.X = r / a;
- rp1.Y = Math.Sqrt(r1 * r1 - rp1.X * rp2.X);
- rp2.Y = -rp1.Y;
- }
- else if (a != 0 && b != 0)
- {
- double delta;
- delta = b * b * r * r - (a * a + b * b) * (r * r - r1 * r1 * a * a);
- rp1.Y = (b * r + Math.Sqrt(delta)) / (a * a + b * b);
- rp2.Y = (b * r - Math.Sqrt(delta)) / (a * a + b * b);
- rp1.X = (r - b * rp1.Y) / a;
- rp2.X = (r - b * rp2.Y) / a;
- }
- rp1.X += p1.X;
- rp1.Y += p1.Y;
- rp2.X += p1.X;
- rp2.Y += p1.Y;
- }
- public static void Swap<T>(ref T a, ref T b)
- {
- T c = a;
- a = b;
- b = c;
- }
- //两圆公共面积:
- // 必须保证相交
- public static double C2Area(PointD p1, double r1, PointD p2, double r2)
- {
- PointD rp1 = new PointD(), rp2 = new PointD();
- C2Point(p1, r1, p2, r2, ref rp1, ref rp2);
- if (r1 > r2) //保证r2>r1
- {
- Swap(ref p1, ref p2);
- Swap(ref r1, ref r2);
- }
- double a, b, rr;
- a = p1.X - p2.X;
- b = p1.Y - p2.Y;
- rr = Math.Sqrt(a * a + b * b);
- double dx1, dy1, dx2, dy2;
- double sita1, sita2;
- dx1 = rp1.X - p1.X;
- dy1 = rp1.Y - p1.Y;
- dx2 = rp2.X - p1.X;
- dy2 = rp2.Y - p1.Y;
- sita1 = Math.Acos((dx1 * dx2 + dy1 * dy2) / r1 / r1);
- dx1 = rp1.X - p2.X;
- dy1 = rp1.Y - p2.Y;
- dx2 = rp2.X - p2.X;
- dy2 = rp2.Y - p2.Y;
- sita2 = Math.Acos((dx1 * dx2 + dy1 * dy2) / r2 / r2);
- double s = 0;
- if (rr < r2) //相交弧为优弧
- s = r1 * r1 * (PI - sita1 / 2 + Math.Sin(sita1) / 2) + r2 * r2 * (sita2 - Math.Sin(sita2)) / 2;
- else //相交弧为劣弧
- s = (r1 * r1 * (sita1 - Math.Sin(sita1)) + r2 * r2 * (sita2 - Math.Sin(sita2))) / 2;
- return s;
- }
- //圆和直线关系:
- //0----相离 1----相切 2----相交
- public static int ClPoint(PointD p, double r, double a, double b, double c, ref PointD rp1, ref PointD rp2)
- {
- int res = 0;
- c = c + a * p.X + b * p.Y;
- double tmp;
- if (a == 0 && b != 0)
- {
- tmp = -c / b;
- if (r * r < tmp * tmp)
- res = 0;
- else if (r * r == tmp * tmp)
- {
- res = 1;
- rp1.Y = tmp;
- rp1.X = 0;
- }
- else
- {
- res = 2;
- rp1.Y = rp2.Y = tmp;
- rp1.X = Math.Sqrt(r * r - tmp * tmp);
- rp2.X = -rp1.X;
- }
- }
- else if (a != 0 && b == 0)
- {
- tmp = -c / a;
- if (r * r < tmp * tmp)
- res = 0;
- else if (r * r == tmp * tmp)
- {
- res = 1;
- rp1.X = tmp;
- rp1.Y = 0;
- }
- else
- {
- res = 2;
- rp1.X = rp2.X = tmp;
- rp1.Y = Math.Sqrt(r * r - tmp * tmp);
- rp2.Y = -rp1.Y;
- }
- }
- else if (a != 0 && b != 0)
- {
- double delta;
- delta = b * b * c * c - (a * a + b * b) * (c * c - a * a * r * r);
- if (delta < 0)
- res = 0;
- else if (delta == 0)
- {
- res = 1;
- rp1.Y = -b * c / (a * a + b * b);
- rp1.X = (-c - b * rp1.Y) / a;
- }
- else
- {
- res = 2;
- rp1.Y = (-b * c + Math.Sqrt(delta)) / (a * a + b * b);
- rp2.Y = (-b * c - Math.Sqrt(delta)) / (a * a + b * b);
- rp1.X = (-c - b * rp1.Y) / a;
- rp2.X = (-c - b * rp2.Y) / a;
- }
- }
- rp1.X += p.X;
- rp1.Y += p.Y;
- rp2.X += p.X;
- rp2.Y += p.Y;
- return res;
- }
- //内切圆:
- public static void InCircle(PointD p1, PointD p2, PointD p3, ref PointD rp, ref double r)
- {
- double dx31, dy31, dx21, dy21, d31, d21, a1, b1, c1;
- dx31 = p3.X - p1.X;
- dy31 = p3.Y - p1.Y;
- dx21 = p2.X - p1.X;
- dy21 = p2.Y - p1.Y;
- d31 = Math.Sqrt(dx31 * dx31 + dy31 * dy31);
- d21 = Math.Sqrt(dx21 * dx21 + dy21 * dy21);
- a1 = dx31 * d21 - dx21 * d31;
- b1 = dy31 * d21 - dy21 * d31;
- c1 = a1 * p1.X + b1 * p1.Y;
- double dx32, dy32, dx12, dy12, d32, d12, a2, b2, c2;
- dx32 = p3.X - p2.X;
- dy32 = p3.Y - p2.Y;
- dx12 = -dx21;
- dy12 = -dy21;
- d32 = Math.Sqrt(dx32 * dx32 + dy32 * dy32);
- d12 = d21;
- a2 = dx12 * d32 - dx32 * d12;
- b2 = dy12 * d32 - dy32 * d12;
- c2 = a2 * p2.X + b2 * p2.Y;
- rp.X = (c1 * b2 - c2 * b1) / (a1 * b2 - a2 * b1);
- rp.Y = (c2 * a1 - c1 * a2) / (a1 * b2 - a2 * b1);
- r = Math.Abs(dy21 * rp.X - dx21 * rp.Y + dx21 * p1.Y - dy21 * p1.X) / d21;
- }
- //求切点:
- // p---圆心坐标, r---圆半径, sp---圆外一点, rp1,rp2---切点坐标
- public static void CutPoint(PointD p, double r, PointD sp, ref PointD rp1, ref PointD rp2)
- {
- PointD p2 = new PointD();
- p2.X = (p.X + sp.X) / 2;
- p2.Y = (p.Y + sp.Y) / 2;
- double dx2, dy2, r2;
- dx2 = p2.X - p.X;
- dy2 = p2.Y - p.Y;
- r2 = Math.Sqrt(dx2 * dx2 + dy2 * dy2);
- C2Point(p, r, p2, r2, ref rp1, ref rp2);
- }
- //线段的左右旋:
- /* l2在l1的左/右方向(l1为基准线)
- 返回 0 : 重合;
- 返回 1 : 右旋;
- 返回 –1 : 左旋;
- */
- public static int Rotat(LineSeg l1, LineSeg l2)
- {
- double dx1, dx2, dy1, dy2;
- dx1 = l1.Start.X - l1.End.X;
- dy1 = l1.Start.Y - l1.End.Y;
- dx2 = l2.Start.X - l2.End.X;
- dy2 = l2.Start.Y - l2.End.Y;
- double d;
- d = dx1 * dy2 - dx2 * dy1;
- if (d == 0)
- return 0;
- else if (d > 0)
- return -1;
- else
- return 1;
- }
- }
- }
|