golang 二维平面求多重遮挡三角形总面积
解决问题描述：二维平面有很多三角形错落，可能会相互叠加落在一起，也可能互相远离。目标求出这些三角形的总占地面积。
我最开始想的解决方案是用总面积-总重叠面积 = 总占地面积。后来实现起来发现当面临多次重叠时，这样算是错误的。
后来参考了一些文献，得出结论：这个问题抽象出来就是求n个集合的并集问题。公式如下：
A1∪A2∪......∪An = A1 + A2 + ......+ An ﹣(A1∩A2 + A1∩A3 + ......+ A[n-1]∩An) + (A1∩A2∩A3 + A1∩A2∩A4 + ......+ A[n-2]∩A[n-1]∩An) + ......+ [(-1)^(n-1)]A1∩A2∩......∩An
第二个需要解决的问题是如何判断两个三角形关系并求两个三角形的重叠面积
思路如下：
第一步：两个三角形各边分别求交点，交点放入集合a中
第二步：两个三角形分别求一个点是否在另一个三角形内部，内点放入集合a中
第三步：集合a去重
第四部：如果a中元素数量>=3，则两三角形相交，重叠部分面积为a中点集所构成的凸多边形的面积。如果a中元素数量<3，则两三角形不相交。
这两个问题解决，问题就迎刃而解了，中间还涉及了很多小函数，整体代码如下：
type Vector2 struct { X float64 `json:"x"` Y float64 `json:"y"` } type Triangle struct { A Vector2 B Vector2 C Vector2 } type Line struct { X Vector2 Y Vector2 } //全局变量 var triangleArea = float64(0) //所有三角片面积 var n_getAcreage = 0 func GetAcreage(triangles []Triangle) float64 { triangleArea = 0 n_getAcreage = len(triangles) for i:=0;i<n_getAcreage-1 ;i++ { for j:=i+1;j< n_getAcreage;j++ { triangleA:=triangles[i] triangleB:=triangles[j] if QuickCross(TriangleToList(triangleA),TriangleToList(triangleB)) { //快速判断，不准确，后面还需再次判断 GetShadow(TriangleToList(triangleA),TriangleToList(triangleB),triangles,1,j+1) } } } } func GetShadow(list1 []Vector2,list2 []Vector2,trangles []Triangle,sign float64,index int){ //两个多边形求遮挡多边形和遮挡面积 var list []Vector2 //相交点集 list1 = ClockwiseSortPoints(list1) list2 = ClockwiseSortPoints(list2) //求相交点 for i:=0;i< len(list1);i++ { for j:=0;j<len(list2);j++ { var line1,line2 Line if i==len(list1)-1 { line1.X = list1[0] line1.Y = list1[i] }else { line1.X = list1[i] line1.Y = list1[i+1] } if j==len(list2)-1 { line2.X = list2[0] line2.Y = list2[j] }else { line2.X = list2[j] line2.Y = list2[j+1] } point := CrossPoint(line1.X,line1.Y,line2.X,line2.Y) if point[0].X!=-1000 && point[0].Y!=-1000{ list = append(list, point...) } } } //求内部点 for _,v1 := range list1{ if IsInNedge(list2,v1) { list = append(list, v1) } } for _,v2 := range list2{ if IsInNedge(list1,v2) { list = append(list, v2) } } //去重 list = DeduplicationVec2(list) if len(list)<3 { return } shadow := GetNedgeArea(ClockwiseSortPoints(list)) triangleArea = triangleArea+math.Pow(-1,sign)*shadow for k:=index;k<n_getAcreage ;k++ { if QuickCross(list,TriangleToList(trangles[k])){ GetShadow(list,TriangleToList(trangles[k]),trangles,sign+1,k+1) } } return } func (this *Vector2)Set(x, y float64) { this.X = x this.Y = y } func (this *Vector2)Equal(x Vector2) bool { if math.Abs(this.X - x.X)<0.1 && math.Abs(this.Y - x.Y)<0.1{ return true } else { return false } } func TriangleToList(triangle Triangle) []Vector2 { var res []Vector2 var a Vector2 a.Set(triangle.A.X,triangle.A.Y) res = append(res, a) a.Set(triangle.B.X,triangle.B.Y) res = append(res, a) a.Set(triangle.C.X,triangle.C.Y) res = append(res, a) return res } func QuickCross(list1,list2 []Vector2) bool { //判断两多边形是否相交 var maxx1,maxy1,minx1,miny1,maxx2,maxy2,minx2,miny2 float64 for _,v1 :=range list1{ if v1.X>=maxx1 { maxx1 = v1.X } if v1.X<=minx1 { minx1 = v1.X } if v1.Y>=maxy1 { maxy1 = v1.Y } if v1.Y<=miny1 { miny1 = v1.Y } } for _,v2 :=range list2{ if v2.X>=maxx2 { maxx2 = v2.X } if v2.X<=minx2 { minx2 = v2.X } if v2.Y>=maxy2 { maxy2 = v2.Y } if v2.Y<=miny2 { miny2 = v2.Y } } Lx := math.Abs((minx1+maxx1)/2-(minx2+maxx2)/2) Ly := math.Abs((miny1+maxy1)/2-(miny2+maxy2)/2) if Lx<math.Abs(maxx1-minx1)+math.Abs(maxx2-minx2) && Ly<math.Abs(maxy1-miny1)+math.Abs(maxy2-miny2) { return true }else { return false } } func GetNedgeArea(list []Vector2) float64 { //得到任意N边型面积 前提是点的排列方向为顺时针或逆时针 var x,y float64 for _,v := range list{ x+=v.X y+=v.Y } var point Vector2 point.Set(x/ float64(len(list)),y/ float64(len(list))) res := float64(0) for q:=0;q< len(list);q++ { if q==len(list)-1 { res+= GetTriangleAreaByVector(point,list[q],list[0]) }else { res+= GetTriangleAreaByVector(point,list[q],list[q+1]) } } return res } func ClockwiseSortPoints(list []Vector2) []Vector2 { //多边形点集排序--针对凸多边形，按逆时针方向进行排序 var res []Vector2 var center Vector2 var a Vector2 var first []Vector2 var second []Vector2 var third []Vector2 var four []Vector2 a.Set(1,0) x := float64(0) y := float64(0) for i:=0;i < len(list);i++{ x += list[i].X y += list[i].Y } center.X = x/ float64(len(list)) center.Y = y/ float64(len(list)) for _,v := range list{ var b Vector2 b.Set(v.X-center.X,v.Y-center.Y) if b.X>=0 && b.Y>0 { first = append(first, b) } if b.X<0 && b.Y>=0 { second = append(second, b) } if b.X<=0 && b.Y<0 { third = append(third, b) } if b.X>0 && b.Y<=0 { four = append(four, b) } } if first!=nil { sort1(first) for _,v := range first{ v.X = v.X+center.X v.Y = v.Y+center.Y res = append(res, v) } } if second!=nil { sort1(second) for _,v := range second{ v.X = v.X+center.X v.Y = v.Y+center.Y res = append(res, v) } } if third!=nil { sort2(third) for _,v := range third{ v.X = v.X+center.X v.Y = v.Y+center.Y res = append(res, v) } } if four!=nil { sort2(four) for _,v := range four{ v.X = v.X+center.X v.Y = v.Y+center.Y res = append(res, v) } } return res } func sort1(buf[]Vector2) { //从小到大排序 var a Vector2 a.Set(1,0) for i := 0; i < len(buf)-1; i++ { flag := false for j := 1; j < len(buf)-i; j++ { if GetAngelByVector2(buf[j-1],a) > GetAngelByVector2(buf[j],a) { tmp := buf[j-1] buf[j-1] = buf[j] buf[j] = tmp flag = true } } if !flag { break } } } func sort2(buf[]Vector2) { //从大到小排序 var a Vector2 a.Set(1,0) for i := 0; i < len(buf)-1; i++ { flag := false for j := 1; j < len(buf)-i; j++ { if GetAngelByVector2(buf[j-1],a) < GetAngelByVector2(buf[j],a) { tmp := buf[j-1] buf[j-1] = buf[j] buf[j] = tmp flag = true } } if !flag { break } } } func GetAngelByVector2(a,b Vector2) float64 { //二维得到两向量夹角 a_sqrt := math.Sqrt(a.X*a.X + a.Y*a.Y) b_sqrt := math.Sqrt(b.X*b.X + b.Y*b.Y) cos := (a.X*b.X + a.Y*b.Y)/(a_sqrt*b_sqrt) res := math.Acos(cos) return res } func DeduplicationVec2(crossPoint []Vector2) []Vector2 { //去重 if crossPoint==nil { return crossPoint } var res []Vector2 var temp = true res = append(res, crossPoint[0]) for _,v1 := range crossPoint{ for _,v2:= range res{ if v2.Equal(v1) { temp = false break } } if temp { res = append(res, v1) } temp = true } return res } func IsInNedge(list []Vector2,point Vector2) bool { nArea := GetNedgeArea(list) var a float64 for i:=0;i<len(list) ;i++ { if i==len(list)-1 { a+=GetTriangleAreaByVector(point,list[i],list[0]) }else { a+=GetTriangleAreaByVector(point,list[i],list[i+1]) } } if math.Abs(a-nArea)<0.1 { return true }else { return false } } func CrossPoint(line1,line2,line3,line4 Vector2) []Vector2 { var res []Vector2 if line1.Equal(line3) && line2.Equal(line4) { res = append(res, line1) res = append(res, line2) return res } if line1.Equal(line4) && line2.Equal(line3) { res = append(res, line1) res = append(res, line2) return res } if line1.Equal(line3) && !line2.Equal(line4) { res = append(res,line1) return res } if line1.Equal(line4) && !line2.Equal(line3) { res = append(res,line1) return res } if line2.Equal(line3) && !line1.Equal(line4) { res = append(res,line2) return res } if line2.Equal(line4) && !line1.Equal(line3) { res = append(res,line2) return res } var cross_point Vector2 var a = float64(0) var b = float64(0) var state = 0 if math.Abs(line1.X - line2.X)>1e-6 { a = (line2.Y - line1.Y) / (line2.X - line1.X) state |= 1 } if math.Abs(line3.X - line4.X)>1e-6 { b = (line4.Y - line3.Y) / (line4.X - line3.X) state |= 2 } switch state { case 0: //L1与L2都平行Y轴 if math.Abs(line1.X - line3.X)<1e-6{ //throw new Exception("两条直线互相重合，且平行于Y轴，无法计算交点。"); cross_point.Set(-1000,-1000) }else{ //throw new Exception("两条直线互相平行，且平行于Y轴，无法计算交点。"); cross_point.Set(-1000,-1000) } case 1: //L1存在斜率, L2平行Y轴 x := line3.X y := (line1.X - x) * (-a) + line1.Y cross_point.Set(x,y) case 2://L1 平行Y轴，L2存在斜率 x := line1.X y := (line3.X - x) * (-b) + line3.Y cross_point.Set(x,y) case 3://L1，L2都存在斜率 if math.Abs(a - b)<1e-6 { // throw new Exception("两条直线平行或重合，无法计算交点。"); cross_point.Set(-1000,-1000) } x := (a * line1.X - b * line3.X - line1.Y + line3.Y) / (a - b) y := a * x - a * line1.X + line1.Y cross_point.Set(x,y) } //beego.Debug(cross_point) if !((OnSegment(line1,line2,cross_point)) && (OnSegment(line3,line4,cross_point))){ cross_point.Set(-1000,-1000) } res = append(res, cross_point) return res //没有交点返回(-1000,-1000) } func GetTriangleAreaByVector(x Vector2,y Vector2,z Vector2) float64 { //根据三点坐标获得三角形面积 area := (x.X*y.Y+y.X*z.Y+z.X*x.Y-x.X*z.Y-y.X*x.Y-z.X*y.Y)/2 if area<0 { area = -area } return area } func OnSegment( p1,p2,Q Vector2) bool { //判断一个点是不是在一个线段内 maxx := math.Max(p1.X,p2.X) minx := math.Min(p1.X,p2.X) maxy := math.Max(p1.Y,p2.Y) miny := math.Min(p1.Y,p2.Y) index := (Q.X -p1.X )*(p2.Y -p1.Y) - (p2.X -p1.X) *(Q.Y -p1.Y) if index<=1e-5 && ( (Q.X > minx || math.Abs(Q.X-minx)<1e-5) && (Q.X < maxx || math.Abs(Q.X-maxx)<1e-5)&& (Q.Y > miny || math.Abs(Q.Y-miny)<1e-5) && (Q.Y < maxy || math.Abs(Q.Y-maxy)<1e-5)) { return true }else { return false } }
参考文献：徐元铭,龙伟,王永庆.军机易损性分析中多重遮挡投影面积计算[J].北京航空航天大学学报,2002(02):245-248.
