4

必ずしも凸であるとは限らない、交点のない多角形と、この多角形の外側の点があります。2 次元空間で最も効率的にユークリッド距離を計算する方法を考えています。に標準的な方法はありRますか?

私の最初のアイデアは、多角形のすべての線の最小距離を計算し (無限に拡張されるため、線分ではなく線分になります)、次に、線分とピタゴラスの開始点を使用して、点から個々の線までの距離を計算することでした。

効率的なアルゴリズムを実装するパッケージについて知っていますか?

4

4 に答える 4

8

rgeos パッケージとメソッドを使用できますgDistance。これには、ジオメトリを準備し、持っているデータからオブジェクトを作成する必要がありますspgeom(data.frame などと仮定します)。rgeos のドキュメントは非常に詳細です (CRAN ページのパッケージの PDF マニュアルを参照してください)。これは、gDistanceドキュメントの関連する例の 1 つです。

pt1 = readWKT("POINT(0.5 0.5)")
pt2 = readWKT("POINT(2 2)")
p1 = readWKT("POLYGON((0 0,1 0,1 1,0 1,0 0))")
p2 = readWKT("POLYGON((2 0,3 1,4 0,2 0))")
gDistance(pt1,pt2)
gDistance(p1,pt1)
gDistance(p1,pt2)
gDistance(p1,p2)

readWKTrgeos にも含まれています。

Rgeos は、ジオメトリック コンピューティングのデファクト スタンダードの 1 つである GEOS ライブラリに基づいています。車輪を再発明したくない場合は、これが良い方法です。

于 2013-03-08T13:19:27.367 に答える
6

後世のために、戻って理論的な解決策を書き留めることにしました。これは最も簡潔な例ではありませんが、このような問題を手動で解決する方法を知りたい人にとっては完全に透過的です。

理論上のアルゴリズム

まず、私たちの仮定です。

  1. 多角形の頂点は、時計回りまたは反時計回りの回転順序で多角形の点を指定し、多角形の線は交差できないと仮定します。これは、奇妙に定義されたベクトル グラフィック形状ではなく、通常の幾何学的多角形があることを意味します。
  2. これは、2 次元平面上の位置を表す「x」値と「y」値を使用したデカルト座標のセットであると想定しています。
  3. ポイントは、ポリゴンの内部領域の外側にある必要があると想定しています。
  4. 最後に、必要な距離は、ポイントとポリゴンの周囲にある無数のポイントすべてとの間の最小距離であると仮定します。

コーディングする前に、やりたいことを基本的な言葉で書き出す必要があります。多角形と多角形の外側の点との間の最短距離は、常に、多角形の頂点または 2 つの頂点間の線上の点のいずれかであると想定できます。これを念頭に置いて、次の手順を実行します。

  1. すべての頂点とターゲット ポイントの間の距離を計算します。
  2. ターゲット ポイントに最も近い 2 つの頂点を見つけます。
  3. (a) 2 つの最も近い頂点が隣接していないか、(b) 2 つの頂点のいずれかの内角が 90 度以上である場合、最も近い頂点が最も近い点になります。最も近いポイントとターゲット ポイントの間の距離を計算します。
  4. それ以外の場合は、2 点間で形成される三角形の高さを計算します。

基本的には、頂点がポイントに最も近いかどうか、または線上のポイントがポイントに最も近いかどうかを確認するだけです。これを機能させるには、いくつかの三角関数を使用する必要があります。

コード

これを適切に機能させるには、「for」ループを回避し、ポリゴン頂点のリスト全体を表示するときにベクトル化された関数のみを使用する必要があります。幸いなことに、これは R では非常に簡単です。ポリゴンの頂点として「x」列と「y」列を持つデータ フレームを受け入れ、点の位置として 1 つの「x」値と「y」値を持つベクトルを受け入れます。

get_Point_Dist_from_Polygon <- function(.polygon, .point){

    # Calculate all vertex distances from the target point.
    vertex_Distance <- sqrt((.point[1] - .polygon$x)^2 + (.point[2] - .polygon$y)^2)

    # Select two closest vertices.
    min_1_Index <- which.min(vertex_Distance)
    min_2_Index <- which.min(vertex_Distance[-min_1_Index])

    # Calculate lengths of triangle sides made of
    # the target point and two closest points.
    a <- vertex_Distance[min_1_Index]
    b <- vertex_Distance[min_2_Index]
    c <- sqrt(diff(.polygon$x[c(min_1_Index, min_2_Index)])^2 + diff(.polygon$y[c(min_1_Index, min_2_Index)])^2)

    if(abs(min_1_Index - min_2_Index) != 1 |
        acos((b^2 + c^2 - a^2)/(2*b*c)) >= pi/2 | 
        acos((a^2 + c^2 - b^2)/(2*a*c)) >= pi/2
        ){
        # Step 3 of algorithm.
        return(vertex_Distance[min_1_Index])
    } else {
        # Step 4 of algorithm.
        # Here we are using the law of cosines.
        return(sqrt((a+b-c) * (a-b+c) * (-a+b+c) * (a+b+c)) / (2 * c))
    }
}

デモ

polygon <- read.table(text="
x,  y
0,  1
1,  0.8
2,  1.3
3,  1.4
2.5,0.3
1.5,0.5
0.5,0.1", header=TRUE, sep=",")

point <- c(3.2, 4.1)

get_Point_Dist_from_Polygon(polygon, point)
# 2.707397
于 2013-03-08T16:08:07.560 に答える
1

さもないと：

p2poly <- function(pt, poly){
    # Closing the polygon
    if(!identical(poly[1,],poly[nrow(poly),])){poly<-rbind(poly,poly[1,])}
    # A simple distance function
    dis <- function(x0,x1,y0,y1){sqrt((x0-x1)^2 +(y0-y1)^2)}
    d <- c()   # Your distance vector
    for(i in 1:(nrow(poly)-1)){
        ba <- c((pt[1]-poly[i,1]),(pt[2]-poly[i,2])) #Vector BA
        bc <- c((poly[i+1,1]-poly[i,1]),(poly[i+1,2]-poly[i,2])) #Vector BC
        dbc <- dis(poly[i+1,1],poly[i,1],poly[i+1,2],poly[i,2]) #Distance BC
        dp <- (ba[1]*bc[1]+ba[2]*bc[2])/dbc          #Projection of A on BC
        if(dp<=0){ #If projection is outside of BC on B side
            d[i] <- dis(pt[1],poly[i,1],pt[2],poly[i,2])
            }else if(dp>=dbc){ #If projection is outside of BC on C side
                d[i] <- dis(poly[i+1,1],pt[1],poly[i+1,2],pt[2])
                }else{ #If projection is inside of BC
                    d[i] <- sqrt(abs((ba[1]^2 +ba[2]^2)-dp^2))
                    }
        }
    min(d)
    }

例：

pt <- c(3,2)
triangle <- matrix(c(1,3,2,3,4,2),byrow=T, nrow=3)
p2poly(pt,triangle)
[1] 0.3162278
于 2013-03-08T13:46:11.773 に答える
0

distm()パッケージ内の関数を使用geosphereして、点と頂点が座標系で表示されている場合の距離を計算しました。dis <- function(x0,x1,y0,y1){sqrt((x0-x1)^2 +(y0-y1)^2)} また、 を代入することで簡単に変更することができますdistm()。

algo.p2poly <- function(pt, poly){
  if(!identical(poly[1,],poly[nrow(poly),])){poly<-rbind(poly,poly[1,])}
  library(geosphere)
  n <- nrow(poly) - 1
  pa <- distm(pt, poly[1:n, ])
  pb <- distm(pt, poly[2:(n+1), ])
  ab <- diag(distm(poly[1:n, ], poly[2:(n+1), ]))
  p <- (pa + pb + ab) / 2
  d <- 2 * sqrt(p * (p - pa) * (p - pb) * (p - ab)) / ab
  cosa <- (pa^2 + ab^2 - pb^2) / (2 * pa * ab)
  cosb <- (pb^2 + ab^2 - pa^2) / (2 * pb * ab)
  d[which(cosa <= 0)] <- pa[which(cosa <= 0)]
  d[which(cosb <= 0)] <- pb[which(cosb <= 0)]
  return(min(d))
}

例：

poly <- matrix(c(114.33508, 114.33616,
                 114.33551, 114.33824,
                 114.34629, 114.35053,
                 114.35592, 114.35951, 
                 114.36275, 114.35340, 
                 114.35391, 114.34715,
                 114.34385, 114.34349,
                 114.33896, 114.33917,
                 30.48271, 30.47791,
                 30.47567, 30.47356, 
                 30.46876, 30.46851,
                 30.46882, 30.46770, 
                 30.47219, 30.47356,
                 30.47499, 30.47673,
                 30.47405, 30.47723, 
                 30.47872, 30.48320),
               byrow = F, nrow = 16)
pt1 <- c(114.33508, 30.48271)
pt2 <- c(114.6351, 30.98271)
algo.p2poly(pt1, poly)
algo.p2poly(pt2, poly)

結果：

> algo.p2poly(pt1, poly)
[1] 0
> algo.p2poly(pt2, poly)
[1] 62399.81
于 2017-04-12T01:52:01.920 に答える