GIS演算法基礎(二)計算幾何基礎(中)
地理資料在計算機中表示大致分為兩種,向量資料和柵格資料。
要計算地理資料的空間關係,一般是向量資料之間的比較。例如:點,線,面之間的比較。
如何判斷線段之間是否相交,線段與面的包含關係。點與面的包含關係等等這些空間關係,都用到計算幾何的演算法。
空間關係的判定演算法的內容有:
1、線段的拐向判斷
2、判斷兩線段是否相交
3、判斷矩形是否包含點
4、判斷線段,折線,多邊形是否在矩形中
5、判斷矩形是否在矩形中
6、判斷圓是否在矩形中
7、判斷點是否在多邊形內
8、判斷線段是否在多邊形內
9、判斷折線是否在多邊形內
。。。。。。等等
話不多說。一個個的來看。
一、線段的拐向的判斷
是現有p,q兩個線段,如何判斷p和q的方位問題呢?
可以利用向量的叉積判斷:
二維平面向量很的叉積很好計算 :
例如: p=(x1,y1),q=(x2,y2)
則 p×q=x1*y2-x2*y1;
有這麼三個關係:
若PxQ>0,則說明P在Q的順時針方向
若PxQ<0,則說明P在Q的逆時針方向
若PxQ=0,則說明PQ共線(共線有可能反向也可能正向)
//2個向量的向量積(叉積) public double crossProduct(Vector2D v) { return x * v.y - y * v.x; }
Vector2D是我構造的工具類,用於向量的表示
package math; import scau.gz.zhw.Line; import scau.gz.zhw.Point; //平面向量(x,y)的基本運算規則,角度弧度的轉換等實現 public class Vector2D { private double x; private double y; public Vector2D() { x = 0; y = 0; } public Vector2D(double _x, double _y) { x = _x; y = _y; } /** * * @param a 起點 * @param b 終點 */ public Vector2D(Point a,Point b) { this.x=b.getX()-a.getX(); this.y=b.getY()-a.getY(); } public Vector2D(Line line) { this(line.getStart(), line.getEnd()); } //獲取弧度 public double getRadian() { return Math.atan2(y, x); } //獲取角度 public double getAngle() { return getRadian() / Math.PI * 180; } public Vector2D clone() { return new Vector2D(x,y); } public double getLength() { return Math.sqrt(getLengthSQ()); } public double getLengthSQ() { return x * x + y * y; } //向量置零 public Vector2D Zero() { x = 0; y = 0; return this; } public boolean isZero() { return x == 0 && y == 0; } //向量的長度設定為我們期待的value public void setLength(double value) { double _angle = getAngle(); x = Math.cos(_angle) * value; y = Math.sin(_angle) * value; } //向量的標準化(方向不變,長度為1) public Vector2D normalize() { double length = getLength(); x = x / length; y = y / length; return this; } //是否已經標準化 public boolean isNormalized() { return getLength() == 1.0; } //向量的方向翻轉 public Vector2D reverse() { x = -x; y = -y; return this; } //2個向量的數量積(點積) public double dotProduct(Vector2D v) { return x * v.x + y * v.y; } //2個向量的向量積(叉積) public double crossProduct(Vector2D v) { return x * v.y - y * v.x; } //計算2個向量的夾角弧度 //參考點積公式:v1 * v2 = cos<v1,v2> * |v1| *|v2| public static double radianBetween(Vector2D v1, Vector2D v2) { if(!v1.isNormalized()) v1 = v1.clone().normalize(); // |v1| = 1 if(!v2.isNormalized()) v2 = v2.clone().normalize(); // |v2| = 1 return Math.acos(v1.dotProduct(v2)); } //弧度 = 角度乘以PI後再除以180、 推理可得弧度換算角度的公式 //弧度轉角度 public static double radian2Angle(double radian) { return radian / Math.PI * 180; } //向量加 public Vector2D add(Vector2D v) { return new Vector2D(x + v.x, y + v.y); } //向量減 public Vector2D subtract(Vector2D v) { return new Vector2D(x - v.x, y - v.y); } //向量乘 public Vector2D multiply(double value) { return new Vector2D(x * value, y * value); } //向量除 public Vector2D divide(double value) { return new Vector2D(x / value, y / value); } public double getX() { return x; } public double getY() { return y; } public Line toLine() { return new Line(new Point(0, 0),new Point(x, y)); } public void showGUI() { toLine().showGUI(); } @Override public String toString() { // TODO Auto-generated method stub return String.format("(%.2f,%.2f)", x,y); } }
判斷叉積函式:
public static void getDirection(Vector2D p,Vector2D q) {
if(p.crossProduct(q)>0) {
System.out.println("順時針");
}else if (p.crossProduct(q)<0) {
System.out.println("逆時針");
}else {
System.out.println("共線");
}
}
測試結果:
說明:為了方便測試,我用java Swing簡單寫了一個GUI介面 可以繪製向量的點,線,面,用來驗證拐向函式的結果是否正確
①測試資料:p=(100,100) q=(100,30)
②測試資料二:p=(100,100) q=(-100,-100)
控制檯結果:
二、判斷點是否線上段上
設點為Q,線段為P1P2,判斷點Q在該線段上的依據是(Q-P1)X(P2-1) = 0,這樣就保證Q在P1P2這條直線上,但是還是不能保證在P1P2的線段上,所以我們得多加個條件:且Q在P1,P2為對角頂點的矩形內
演算法:
/**
* 判斷點是否線上段上
* @param p1 線段端點
* @param p2 線段端點
* @param q 需要判斷的點
* @return
*/
public static boolean isPointAtSegment(Point p1,Point p2,Point q) {
//判斷點是否線上段圍成的矩形區域內
if(q.getX()<=Math.max(p1.getX(), p2.getX()) && q.getX()>=Math.min(p1.getX(), p2.getX())
&& q.getY()<= Math.max(p1.getY(), p2.getY()) && q.getY()>=Math.min(p1.getY(), p2.getY()))
{
Vector2D qp1 = getVector(q, p1);
Vector2D p2p1 = getVector(p2, p1);
return qp1.crossProduct(p2p1)==0?true:false;
}else {
return false;
}
}
三、判斷兩線段是否相交
判斷兩線段是否相交,我的第一反應就是解方程,看兩線段是否有交點,但是在GIS演算法中,我們面對的是海量的資料,這樣的演算法因為計算量大並不高效。
因此,我們最好得有個篩選的過程,把那些明顯不會相交的線段剔除掉,這樣就減小了一部分的計算量。
其次,我們判斷線段相交最好不要解方程,這樣涉及的計算量大,可以用向量的方法
所以,分為兩步確定兩天線段是否相交:
①快速排斥試驗
設以線段a,b為對角線的矩形為R,設以線段c,d為對角線的矩形為T,如果R和T不相交,顯然兩線段不會相交
②跨立試驗
如果ab,cd相交,那麼ab必跨過cd,那麼(ac x ab )x(bd xab)<=0
因為ab兩邊一定分別有個線段(在ab順時針方向的向量與ab的叉積小於0,在ab逆時針方向的向量與ab叉積大於0),所以乘積是小於0的。至於為什麼等於0,是因為ac,和bd有可能在cd上,即ac與ab共線或bd與ab共線或ab和cd共線,那麼這種情況就是等於0了,也算相交的一種情況
快速排斥試驗:
怎麼判斷兩個矩形相不相交,可以這樣判斷:
①ab的較低點低於cd的較高點 (y值比較)
②cd的左端小於ab的右端(x值比較)
③cd較低點低於ab的較高點(y值比較)
④ab的左端小於cd的右端(x值比較)
演算法實現:
//快速排斥試驗,判斷ab,cd圍成的兩矩形是否相交
if(Math.min(a.getY(), b.getY())<=Math.max(c.getY(), d.getY()) &&
Math.min(c.getX(), d.getX())<=Math.max(a.getX(), b.getX())&&
Math.min(c.getY(), d.getY())<= Math.max(a.getY(), b.getY()) &&
Math.min(a.getX(), b.getX())<= Math.max(c.getX(), d.getX())) {
flag1 = true;
}
跨立試驗
//跨立試驗
if(ac.crossProduct(ab) * bd.crossProduct(ab) <=0) {
flag2 = true;
}
完整演算法在這:
/**
* 判斷兩線段是否相交
* @param a
* @param b
* @param c
* @param d
* @return
*/
public static boolean isTwoSegmentIntersect(Point a,Point b,Point c,Point d) {
boolean flag1 = false;
boolean flag2 = false;
//快速排斥試驗
if(Math.min(a.getY(), b.getY())<=Math.max(c.getY(), d.getY()) && Math.min(c.getX(), d.getX())<=Math.max(a.getX(), b.getX())
&& Math.min(c.getY(), d.getY())<= Math.max(a.getY(), b.getY()) && Math.min(a.getX(), b.getX())<= Math.max(c.getX(), d.getX())) {
flag1 = true;
}
Vector2D ab = getVector(a, b);
Vector2D ac = getVector(a, c);
Vector2D bd = getVector(b, d);
//跨立試驗
if(ac.crossProduct(ab) * bd.crossProduct(ab) <=0) {
flag2 = true;
}
return flag1&&flag2;
}
三、判斷點是否在多邊形內
判斷點是否在多邊形內有兩種方法,一種是射線法,另一種是轉角法。
①射線法:引一條從P開始,穿過多邊形邊界的次數為交點數目。當交點數目為偶數時,點P在多邊形外部
為方便計算選取一條水平的、從被判斷點的右邊延伸的,平行於x軸的射線。
為了使在某些特殊情況下判斷穿越是否有效,有效穿越要符合以下幾個規則:
- 方向向上的邊包括開始點,不包括終止點
- 方向向下的邊包括終止點,不包括開始點
- 水平邊不參與測試
- 射線和多邊形的邊的交點必須嚴格在點P的右邊
- 如果點在多邊形邊上,則直接判斷為在多邊形內部(這一條可以根據不同需要設定為不同的結果)
射線法演算法步驟:
說明:p射線是由點p延伸出的水平x軸正方向射線
- 開始遍歷多邊形的每條邊
- 判斷邊是否水平,如果是就跳出本次迴圈
- 如果點在多邊形上,直接返回true
- 判斷p射線是否在邊的左邊,如果不在,就直接跳過該邊的計算
- 判斷p射線是否穿過邊的端點,是則利用上面所述規則1,2進行穿越測試
- 遍歷結束,根據穿越數得出結果
話不多說,看演算法實現吧:
一、射線法的實現
public static boolean isPointAtPolygon(Point[] points,Point p) {
int crossNum = 0;
boolean flag = false;
if(Double.doubleToLongBits(points[points.length-1].getX())==Double.doubleToLongBits(points[0].getX()) &&
Double.doubleToLongBits(points[points.length-1].getY()) == Double.doubleToLongBits((points[0].getY()))){
flag = true;
}
//如果最後一個點不等於第一個點
//自動閉合
ArrayList<Line> lines = new ArrayList<Line>();
if(!flag) {
for(int i=0;i<points.length;i++) {
Line line = new Line(points[i%points.length], points[(i+1)%points.length]);
lines.add(line);
}
}else {
for(int i=0;i<points.length-1;i++) {
Line line = new Line(points[i], points[i+1]);
lines.add(line);
}
}
//遍歷每條邊
if(type==0) {
for(Line line : lines) {
//rule#1:方向向上的邊包括其開始點,不包括其終止點
//rule#2:方向向下的邊包括其終止點,不包括其開始點
//rule#3:水平邊不參與穿越測試
//rule#4:射線和多邊形的邊的交點必須嚴格在點p的右邊
//如果點線上段上則直接判定為在多邊形內部
if(isPointAtSegment(line.getStart(), line.getEnd(), p)) {
return true;
}
//rule#3
if(!line.isHorizontal()) {
//rule#4,保證p在邊的右邊
if(Double.doubleToLongBits(p.getX()) < Double.doubleToLongBits(line.getXByY(p.getY()))
&& Double.doubleToLongBits(line.getXByY(p.getY()))>=Double.doubleToLongBits(Math.min(line.getStart().getX(),line.getEnd().getX()))
&& Double.doubleToLongBits(line.getXByY(p.getY()))<=Double.doubleToLongBits(Math.max(line.getStart().getX(), line.getEnd().getX()))) {
//rule#1,2
//當點的射線穿過每條邊的端點時,規則1,2起作用
if(Double.doubleToLongBits(p.getY()) == Double.doubleToLongBits(line.getStart().getY())
||Double.doubleToLongBits(p.getY()) == Double.doubleToLongBits(line.getEnd().getY()) ) {
if(line.isUp()) {
//判斷是穿過開始點還是終止點
if(Double.doubleToLongBits(p.getY()) == Double.doubleToLongBits(line.getStart().getY())) {
//方向為上,穿過開始點,則有效穿越
crossNum++;
}else {
//方向為上,穿過終止點,則無效穿越
}
}else {
//判斷是穿過開始點還是終止點
if(Double.doubleToLongBits(p.getY()) == Double.doubleToLongBits(line.getStart().getY())) {
//方向為下,穿過開始點,則無效穿越
}else {
//方向為下,穿過終止點,則有效穿越
crossNum++;
}
}
}else {
//直接計算
++crossNum;
}
}
}
}
//奇數在內,偶數在外
//System.out.println("crossNum : " +crossNum);
return crossNum%2==0?false:true;
}
註釋說的很詳細了,就不加贅述了。
轉角法
轉角法可以很精確地判斷一個點是否在非簡單多邊形內部(就是多邊形內部還有一些複雜的構造,後面有對比說明)。它需要計算多邊形繞點有多少次。如果環繞數為零,那麼點在多邊形外部,非零,則點在多邊形內部。
轉角法也和射線法類似,遵守這麼些規則:
- 方向向上的邊包括開始點,不包括終止點
- 方向向下的邊包括終止點,不包括開始點
- 水平邊不參與測試
- 射線和多邊形的邊的交點必須嚴格在點P的右邊
- 如果點在多邊形邊上,則直接判斷為在多邊形內部(這一條可以根據不同需要設定為不同的結果)
看上面這個多邊形,如果從多邊形的a邊開始遍歷 (水平邊不參與穿越測試),因此a跳過
規定點在邊的左邊(以邊的前進方向為準)時環繞數+1,點在邊的右邊時(以邊的前進方向為準)環繞數-1;
那麼上邊這個p點環繞數(從a邊開始)的計算過程就為(-1,-1,-1,-1) 結果為-4,不等於零,說明在多邊形內部
q點的環繞數計算過程(從a邊開始)為(+1,-1,-1,+1),結果為0,說明在多邊形外部。
關於如何判斷點在邊的右邊還是左邊,可以從點向右做一條水平射線(為了保證邊與射線的交點在點的右邊) 判斷該射線與邊的叉積即可
二、轉角法的實現
public static boolean isPointAtPolygon(Point[] points,Point p){
//轉角法需要判斷從p向右出發的水平射線與線段方向的關係,即p是否在邊的左邊
//int count = 0;
for(Line line:lines) {
//如果點線上段上則直接判定為在多邊形內部
if(isPointAtSegment(line.getStart(), line.getEnd(), p)) {
return true;
}
if(!line.isHorizontal()) {
//構造從p出發的水平向右射線,rule#3
Vector2D pVector2d = new Vector2D(p.getX(),0);
//保證p在邊的左邊,rule#4
if(Double.doubleToLongBits(p.getX()) < Double.doubleToLongBits(line.getXByY(p.getY()))
&& Double.doubleToLongBits(line.getXByY(p.getY()))>=Double.doubleToLongBits(Math.min(line.getStart().getX(),line.getEnd().getX()))
&& Double.doubleToLongBits(line.getXByY(p.getY()))<=Double.doubleToLongBits(Math.max(line.getStart().getX(), line.getEnd().getX()))) {
//System.out.println("id : "+count++ +" direction: "+line.getDirection());
//邊向上,p在邊左邊->p的向右出發的水平射線在邊的順時針方向
if(line.isUp()) {
//rule#1,2
//這種情況只在點的射線穿過每條邊的端點才有用
if(Double.doubleToLongBits(p.getY()) == Double.doubleToLongBits(line.getStart().getY())
||Double.doubleToLongBits(p.getY()) == Double.doubleToLongBits(line.getEnd().getY()) ) {
//判斷是穿過開始點還是終止點
if(Double.doubleToLongBits(p.getY()) == Double.doubleToLongBits(line.getStart().getY())) {
//方向為上,穿過開始點,則有效穿越
if(pVector2d.crossProduct(line.getVector2D())>0) {
crossNum++;
}else {
crossNum--;
}
}else {
//方向為上,穿過終止點,則無效穿越
}
}else {
//不穿過端點的情況
if(pVector2d.crossProduct(line.getVector2D())>0) {
crossNum++;
}else {
crossNum--;
}
}
}
//邊向下,p在邊左邊->p的向右出發的水平射線在邊的順時針方向
if (line.isDown()) {
//rule#1,2
//這種情況只在點的射線穿過每條邊的端點才有用
if(Double.doubleToLongBits(p.getY()) == Double.doubleToLongBits(line.getStart().getY())
||Double.doubleToLongBits(p.getY()) == Double.doubleToLongBits(line.getEnd().getY()) ) {
//判斷是穿過開始點還是終止點
if(Double.doubleToLongBits(p.getY()) == Double.doubleToLongBits(line.getStart().getY())) {
//方向為下,穿過開始點,則無效穿越
}else {
//方向為下,穿過終止點,則有效穿越
if(pVector2d.crossProduct(line.getVector2D())>0) {
crossNum++;
}else {
crossNum--;
}
}
}else {
//不穿過端點的情況
if(pVector2d.crossProduct(line.getVector2D())>0) {
crossNum++;
}else {
crossNum--;
}
}
}
}
}
}
上面也說到,轉角法更為精確一些,為什麼呢?
因為在一些多邊形中可能會有複雜的構造,這種情況下,射線法的計算結果會不準確
例如:
這種情況下
可以看到同樣的多邊形,同樣的點,用不同的方法會得到不同的結果。這是因為轉角法更精確,而射線法的缺點在於沒有考慮到多邊形內部的複雜構造可能使結果出現偏差。
轉角法是在射線法的基礎上進行的優化,他具有和射線法相同的效率,並且,它更為精確,因此,判斷點是否在任意多邊形時,轉角法是首選
射線法的思路比較簡單,只需要判斷穿越次數的奇偶性即可。但需要注意的是,判斷穿越為有效穿越需要嚴格遵守4個規則。但是射線法的缺陷也是很明顯的,它的判斷不夠精確。在左圖的情況下,射線法的結果是false,而轉角法的結果是true。
射線法是沒有考慮到多邊形內部的構造的。而轉角法是射線法的優化(因為他們使用了相同的穿越規則),它通過計算環繞次數來得出結果,會考慮多邊形內部的構造問題。
在同樣的效率下,轉角演算法比射線演算法準確度更高。
同時,我認為這也是個仁者見仁智者見智的過程,因為不同的場景下,可能把多邊形形成的內部構造的認為是屬於該多邊形的,也可能認為是不屬於的。因此具體的使用要根據具體場景而選擇不同的方法。