用C＃实现Hoey Shamos算法

好的，我现在从我当前的算法中获取正确的信息！ 但是，要检查700,000个多边形，这太慢了！ 上一期是固定的（My Line2D intersectsWith方法不正确）

现在这是确定我的瓶颈的问题！ 该算法假设为O（nlog-n），因此它应该更快。 我的intersectsWith方法看起来不能更快，但我会发布它的代码，万一我错了

编辑：添加了IComparable接口

我读取线段交叉点的方法。 为了便于阅读，省略了一些代码。

public class Line2D : IComparable { public Line2D(XYPoints p1, XYPoints p2) { } public bool intersectsLine(Line2D comparedLine) { if ((X2 == comparedLine.X1) && (Y2 == comparedLine.Y1)) return false; if ((X1 == comparedLine.X2) && (Y1 == comparedLine.Y2)) return false; if (X2 == comparedLine.X1 && Y2 == comparedLine.Y1) { return false; } if (X1 == comparedLine.X2 && Y1 == comparedLine.Y2) { return false; } double firstLineSlopeX, firstLineSlopeY, secondLineSlopeX, secondLineSlopeY; firstLineSlopeX = X2 - X1; firstLineSlopeY = Y2 - Y1; secondLineSlopeX = comparedLine.getX2() - comparedLine.getX1(); secondLineSlopeY = comparedLine.getY2() - comparedLine.getY1(); double s, t; s = (-firstLineSlopeY * (X1 - comparedLine.getX1()) + firstLineSlopeX * (getY1() - comparedLine.getY1())) / (-secondLineSlopeX * firstLineSlopeY + firstLineSlopeX * secondLineSlopeY); t = (secondLineSlopeX * (getY1() - comparedLine.getY1()) - secondLineSlopeY * (getX1() - comparedLine.getX1())) / (-secondLineSlopeX * firstLineSlopeY + firstLineSlopeX * secondLineSlopeY); if (s >= 0 && s = 0 && t <= 1) { return true; } return false; // No collision } int IComparable.CompareTo(object obj) { //return Y1.GetHashCode(); Line2D o1 = this; Line2D o2 = (Line2D)obj; if (o1.getY1()  o2.getY2()) { return 1; } else { if (o1.getY2()  o2.getY2()) { return 1; } else { return 0; } } } } 

我的算法实现的大部分，我意识到List不是算法的最快，但我需要索引！：

 //Create a new list, sort by Y values. List SortedList = events.OrderBy(o => o.getY()).ToList(); List sweepline = new List(); for (var g = 0; g  o.getY1()).ToList(); } else { Line2D nl = SortedList[g].line; Line2D above; Line2D below; /* Start generating above */ try { //grab index in sweepline int index = sweepline.IndexOf(nl); Console.Out.WriteLine("index:" + index); //add 1 to get above line above = sweepline[index + 1]; } catch (ArgumentOutOfRangeException) { above = null; } /* End generating above */ /* Start generating below */ try { //grab index in sweepline int index = sweepline.IndexOf(nl); //add 1 to get above line below = sweepline[index - 1]; } catch (ArgumentOutOfRangeException) { below = null; } /* End generating below */ sweepline = sweepline.OrderBy(o => o.getY1()).ToList(); sweepline.Remove(nl); if (above != null && below != null) { if (above.intersectsLine(below)) { return true; } } } Console.WriteLine(""); } } // end numofparts for-loop return false; 

============================================

更新：9月12日：从C5实现了TreeSet，对我的类实现了IComparable，并进一步降低了它的速度？ 如果重要的话，我还在编制索引吗？

http://www.itu.dk/research/c5/

使用TreeSet的代码：

 TreeSet sweepline = new TreeSet(); for (var g = 0; g  o.getY1()).ToList(); } else { Line2D nl = SortedList[g].line; Line2D above; Line2D below; /* Start generating above */ try { //grab index in sweepline int index = sweepline.IndexOf(nl); //Console.Out.WriteLine("index:" + index); //add 1 to get above line above = sweepline[index + 1]; } catch (IndexOutOfRangeException) { above = null; } /* End generating above */ /* Start generating below */ try { //grab index in sweepline int index = sweepline.IndexOf(nl); //add 1 to get above line below = sweepline[index - 1]; } catch (IndexOutOfRangeException) { below = null; } /* End generating below */ //sweepline = sweepline.OrderBy(o => o.getY1()).ToList(); sweepline.Remove(nl); if (above != null && below != null) { if (above.intersectsLine(below)) { return false; } } } //Console.WriteLine(""); 

}

Answer 1

首先，关于线交叉：您不需要实际的交点，只知道它们是否相交。 请参阅http://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/获取可实现此目的的算法。

---

关于List实现：

在使用List s的实现中，您在sweepline上调用indexOf来查找nl 。 这将从开始到结束搜索扫描线。 见List(T).IndexOf 。 如果您使用BinarySearch方法，那么应该大大加快搜索速度。

List的文档有一段称为性能注意事项。 他们敦促您使用实现IEquatable和IComparable的值类型。 所以，你的Line2D应该是一个struct并实现这些接口。

如果您遵循该建议，从扫描线检索端点应该是O（log n），这足以满足您的需要，并且应该更有效地使用内存。

插入和删除对于列表是O（n），因为底层数组需要在内存中移动。 SortedSet有更快的插入和删除，但我不太清楚如何在那里找到O（log n）中项目的邻居。 任何人？ （另请参阅为什么SortedSet .GetViewBetween不是O（log N）？ ）

---

无论如何，C5 TreeSet应该解决这个问题。

我在用户指南中查找了IndexOf和[i]的性能，它们都被列为O（log n）。 所以这不应该是问题。 它仍然可能稍微快一些，但不超过一个固定因素，调用特定的方法来寻找扫描线上的邻居，即Successor和Predecessor ，它们也是O（log n）。

所以

 [...] try { Line2D above = sweepline.Successor(nl); if (above.intersectsLine(nl)) { return false; } } catch (NoSuchItemException ignore) { } [...] 

我不喜欢他们没有一个不抛出exception的方法，因为抛出exception是非常昂贵的。 你的扫描线通常会非常充实，所以我最好的猜测是找不到它会很少见，并且调用Successor是最有效的方法。 或者，您可以像现在一样继续调用IndexOf ，但在检索[index + 1]之前检查它是否等于Count减[index + 1] ，并防止抛出exception：

 [...] int index = sweepline.IndexOf(nl); if( index < sweepline.Count-1 ) { Line2D above = sweepline[index + 1]; if (above.intersectsLine(nl)) { return false; } } [...] 

本手册的第二章介绍了C5集合的相等性和比较。 在这里，至少你必须实现IEquatable和IComparable ！

还有一个想法：你报告喂算算法700000行。 你可以从1000,2500,5000,10000行的时间开始，看看算法如何在它们不相交的情况下进行缩放？

---

关于如何比较扫描线上的线条：

您需要在Sweepline TreeSet上找到Line2D的某种自然排序，因为CompareTo方法会要求您将一个Line2D与另一个Line2D进行比较。 其中一个Line2D已经位于Sweepline TreeSet中，另一个刚刚遇到并正在添加。

我认为你的清扫线从下到上运行：

 List SortedList = events.OrderBy(o => o.getY()).ToList(); 

所以我们假设在事件1中将段S1添加到TreeSet中，我们希望将它与S2进行比较，S2现在正在事件2中添加。

线段可能在某个点交叉，这会改变排序，但算法会在插入之后检查这个，在上面和下面的检查中。 或许更好的称为左右检查，来考虑它。

无论如何..所以最简单的方法是比较两个线段的底端点。 左边较小，右边较大。 但是，我们需要查看扫描线位置的顺序，从那时起它们可能会改变位置，如图所示。

因此，我们需要将S2的底部端点与S1上的红点进行比较，看它是否位于该点的左侧或右侧。 它位于左侧，因此S2被认为小于S1。

通常它比这更简单：如果S1全部位于S2底部端点的左侧，则S1小于S2。 如果S1全部位于S2的底部端点的右侧，则S1大于S2。

我想你正在寻找类型更好的界面版本：

 public class Line2D : IComparable 

假设两个属性BottomY是两个Y值中最低的，而BottomX是最低端点的X值，这是一个有点测试的尝试：

 int IComparable.CompareTo(Line2D other) { if( BottomY < other.BottomY ) { return -other.CompareTo(this); } // we're the segment being added to the sweepline if( BottomX >= other.X1 && BottomX >= other.X2 ) { return 1; } if( BottomX <= other.X1 && BottomX <= other.X2 ) { return -1; } if( other.Y2 == other.Y1 ) { // Scary edge case: horizontal line that we intersect with. Return 0? return 0; } // calculate the X coordinate of the intersection of the other segment // with the sweepline // double redX = other.X1 + // (BottomY - other.Y1) * (other.X2 - other.X1) / (other.Y2 - other.Y1); // // return BottomX.CompareTo(redX); // But more efficient, and more along the lines of the orientation comparison: return Comparer.Default.Compare( (BottomX - other.X1) * (other.Y2 - other.Y1), (BottomY - other.Y1) * (other.X2 - other.X1) ); } 

Answer 2

[原始答案]

我不是C＃用户，但这应该可以加快速度。

• 少堆垃圾
• 不要两次计算相同的东西
• 如果可以，请避免所有子调用（删除函数）

码：

 public bool intersectsLine(const Line2D &comparedLine) { if ((X2==comparedLine.X1)&&(Y2==comparedLine.Y1)) return false; if ((X1==comparedLine.X2)&&(Y1==comparedLine.Y2)) return false; double dx1,dy1,dx2,dy2; dx1 = X2 - X1; dy1 = Y2 - Y1; dx2 = comparedLine.X2 - comparedLine.X1; dy2 = comparedLine.Y2 - comparedLine.Y1; double s,t,ax,ay,b; ax=X1-comparedLine.X1; ay=Y1-comparedLine.Y1; b=1.0/(-(dx2*dy1)+(dx1*dy2)); s = (-(dy1*ax)+(dx1*ay))*b; t = ( (dx2*ay)-(dy2*ax))*b; if ((s>=0)&&(s<=1)&&(t>=0)&&(t<=1)) return true; return false; // No collision } 

对于其余的代码，添加时间测量以找出确切缓慢的事情。 我的猜测是列表管理......不必要的重新分配会大大减慢速度。

[EDIT1]

经过对随机线数据的一些研究后，我得出结论：

1. 如果整个区域的线数太多，则优化没有效果
2. 如果有更多的小行，那么任何优化都会有更多的加速
3. 蛮力是T((N*NN)/2) ，对于700K线待处理（不可用），仍然是O(N*N)估计约35小时

4. 区域细分的优化powershell是T((((N/M)^2)-N)/2) - 优化~O((N/M)^2)其中

   • N是区域行数的最大值

   • M是每个轴的区域划分数，其思想是仅检查穿过某个区域的线（将数据集区域划分为M*M正方形/矩形）。 对于700K线路，最好细分为16x16区域。 测量时间：

     每32K线0.540s每64K线1.950s每128K线7.000s每256K线27.514s

   每700K线的估计运行时间为3.7分钟（对于整个区域的最大~10％长度的线路） 。 我认为这比你的19分钟要好。

5. 使用多CPU /核心可以实现另一种加速

   算法是完全并行的，4个CPU /核心3.7min/4 -> 56s或者你可以将它移植到GPU ...

我的优化powershell算法，区域细分O（（（（（N / M）^ 2）-N）/ 2） - 优化

1. 得到使用面积大小(xmin,xmax,ymin,ymax) O(N)
2. 选择细分M我尝试为我的随机数据集32K-256K行最好的是M=16

3. 循环遍历所有细分区域（均匀划分的数据集区域）

   创建穿过实际细分区域的线条列表，并检查该列表中所有线条的交叉点。 如果不想要重复的交叉点，则丢弃当前区域外的所有交叉点

我的代码（我正在使用BDS2006 C ++和我自己的列表，所以你需要移植它以与你的代码兼容）

 void Twin_GLView2D::main_intersect(int M=16) { int ia,ib,i,j,N; double zero=1e-6; glview2D::_lin *l; glview2D::_pnt p; struct _line { double bx0,by0,bx1,by1; // bounding rectangle double x0,y0,dx,dy; // precomputed params } *lin,*a,*b; struct _siz { double bx0,bx1,by0,by1; // zone bounding rectangle } sz,bz; List<_line*> zone; // load and precompute lines N=view.lin.num; lin=new _line[N]; if (lin==NULL) return; for (a=lin,l=view.lin.dat,ia=0;iap0.p[0]<=l->p1.p[0]) { a->bx0=l->p0.p[0]; a->bx1=l->p1.p[0]; } else { a->bx0=l->p1.p[0]; a->bx1=l->p0.p[0]; } if (l->p0.p[1]<=l->p1.p[1]) { a->by0=l->p0.p[1]; a->by1=l->p1.p[1]; } else { a->by0=l->p1.p[1]; a->by1=l->p0.p[1]; } a->x0=l->p0.p[0]; a->dx=l->p1.p[0]-l->p0.p[0]; a->y0=l->p0.p[1]; a->dy=l->p1.p[1]-l->p0.p[1]; // global image size for zone subdivision if (!ia) { sz.bx0=l->p0.p[0]; sz.by0=l->p0.p[1]; sz.bx1=sz.bx0; sz.by1=sz.by0; } if (sz.bx0>l->p0.p[0]) sz.bx0=l->p0.p[0]; if (sz.bx1p0.p[0]) sz.bx1=l->p0.p[0]; if (sz.by0>l->p0.p[1]) sz.by0=l->p0.p[1]; if (sz.by1p0.p[1]) sz.by1=l->p0.p[1]; if (sz.bx0>l->p1.p[0]) sz.bx0=l->p1.p[0]; if (sz.bx1p1.p[0]) sz.bx1=l->p1.p[0]; if (sz.by0>l->p1.p[1]) sz.by0=l->p1.p[1]; if (sz.by1p1.p[1]) sz.by1=l->p1.p[1]; } // process lines by zonal subdivision zone.allocate(N); view.pnt.num=0; view.pnt.allocate(view.lin.num); sz.bx1-=sz.bx0; sz.bx1/=double(M); sz.by1-=sz.by0; sz.by1/=double(M); for (bz.by0=sz.by0,bz.by1=sz.by0+sz.by1,i=0;ibx0<=bz.bx1)&&(a->bx1>=bz.bx0)) if ((a->by0<=bz.by1)&&(a->by1>=bz.by0)) zone.add(a); // add line to zone list // check for intersection within zone only // O((((N/M)^2)-N)/2) - optimizations for (ia= 0,a=zone.dat[ia];iabx1bx0) continue; if (a->bx0>b->bx1) continue; if (a->by1by0) continue; if (a->by0>b->by1) continue; // 2D lines a,b intersect ? double x0,y0,x1,y1,t0,t1; // compute intersection t1=divide(a->dx*(a->y0-b->y0)+a->dy*(b->x0-a->x0),(a->dx*b->dy)-(b->dx*a->dy)); x1=b->x0+(b->dx*t1); y1=b->y0+(b->dy*t1); if (fabs(a->dx)>=fabs(a->dy)) t0=divide(b->x0-a->x0+(b->dx*t1),a->dx); else t0=divide(b->y0-a->y0+(b->dy*t1),a->dy); x0=a->x0+(a->dx*t0); y0=a->y0+(a->dy*t0); // check if intersection exists if (fabs(x1-x0)>zero) continue; if (fabs(y1-y0)>zero) continue; if ((t0<0.0)||(t0>1.0)) continue; if ((t1<0.0)||(t1>1.0)) continue; // if yes add point pp[0]=x0; pp[1]=y0; pp[2]=0.0; // do not add points out of zone (allmost all duplicit points removal) if (x0bz.bx1) continue; if (y0bz.by1) continue; view.pnt.add(p); } } view.redraw=true; delete lin; } 

笔记：

1. List x; 与“无限”大小的T x[]相同
   • x.num; 是Ts x[]字节x[]实际大小x[] ！ index = <0,x.num-1>
   • x.add(q); 最后在列表中添加q
   • x.num=0; 清除列表
   • x.allocate(N); 为列表中的N项目分配空间以避免重定位
2. 输入List<>是view.lin ...包含点p0,p1每个都有double p[2] ... x,y
3. 输出List<>是view.pnt ...包含double p[2] ... x,y

[EDIT2]

另外我发现上述算法的最佳性能是当M=12+(N>>15)

• M是每个轴的细分区域数
• N是要检查的行数