折れ線を間引く

読み込んだGPSログのデータを間引きたい、と思って調べたところ、 (Ramer-)Douglas-Peuckerのアルゴリズムというものがあることが分かった。

基本的な考え方は、

  1. 折れ線の始点と終点を結ぶ線分と各点の距離を求める。
  2. すべての点との距離が許容誤差ε以内に入っていれば始点と終点だけを返して終了。
  3. そうでなければ距離が最大の点Pを1つ選択。
  4. 始点から点Pまでの折れ線と、点Pから終点までの折れ線のそれぞれについてまた1から処理する。

という再帰的なもの。

再帰的なものはMathematicaの得意分野なので、MathematicaでRamer-Douglas-Peuckerのアルゴリズムを実装してみた。

線分と点の距離

さて、まずは線分と点の距離を求める関数を作成する。

ネットで検索すると直線と点の距離を求めるものが多いが、線分と点の距離のほうが形状をよく保存できるのではないだろうか、と個人的に考えている。

線分ABに点Pから下ろした垂線の足が点Hだとすると、

  • 線分ABと直線PHは直交⇔内積が0
  • 点Hは直線AB上にある⇔実数tに対してB = t H(B、Hは点Aから見た点B、点Hの位置ベクトル)

の条件から、Hの座標と実数tの値が求まる。 Mathematicaを使うと

A = {ax, ay}; B = {bx, by}; P = {px, py}; H = {x, y};

Solve[{
     (H - A).(P - H) == 0,
     H == A + (B - A) t
     }, {x, y, t}][[2]] // FullSimplify
{x -> (ay^2 bx + ax^2 px + ax ay (-by + py) - 
     ay bx (by + py) + ax (by^2 - 2 bx px - by py) + 
     bx (bx px + by py))/((ax - bx)^2 + (ay - by)^2), 
 y -> ((ax - bx) (by (ax - px) + ay (-bx + px)) + (ay - 
      by)^2 py)/((ax - bx)^2 + (ay - by)^2), 
 t -> (ax^2 + ay^2 + bx px - ax (bx + px) + by py - 
   ay (by + py))/((ax - bx)^2 + (ay - by)^2)}

となる。 これは2次元の場合だが、一般に

  • t == (A - B).(A - P) / (A - B).(A - B)
  • H == A + (B - A)t

となるので、3次元以上の場合はこれから計算すればよい。

実数tが0から1までの間にあるとき、垂線の足、点Hは線分AB上にある。

よって線分ABと点Pの距離は

  • 0<=t<=1ならばHP
  • t<0ならばAP
  • 1<tならばBP

となる。以上を実装すると、

SegmentPointDistance[{{ax_, ay_}, {bx_, by_}}, {px_, py_}] := Module[
  {x, y, t},
  t = (ax^2 + ay^2 + bx px - ax (bx + px) + by py - ay (by + py))/(
   ax^2 + ay^2 - 2 ax bx + bx^2 - 2 ay by + by^2);
  x = ax - (ax - bx) t;
  y = ay - (ay - by) t;
  Which[
   TrueQ[0 <= t <= 1], Norm[{x, y} - {px, py}],
   TrueQ[1 < t], Norm[{bx, by} - {px, py}],
   True, Norm[{ax, ay} - {px, py}]
   ]
  ]

となる。本当は点Aと点Bが一致してしまう場合を回避しなければならない。

折れ線を間引く

これで本題に入れる。とはいっても最初のアルゴリズムを素直に実装すればよく、

LineSimplify[{p1 : {_, _}, midpoints : {_, _} .., p2 : {_, _}}, \[Epsilon]_] := Module[
  {distance, maxdistance, maxdistanceposition},
  distance = SegmentPointDistance[{p1, p2}, #] & /@ {midpoints};
  maxdistanceposition = Ordering[distance, -1][[1]];
  maxdistance = distance[[maxdistanceposition]];
  If[TrueQ[maxdistance < \[Epsilon]],
   {p1, p2},
   Join[
    LineSimplify[
     Join[{p1}, {midpoints}[[1 ;; maxdistanceposition]]], \[Epsilon]],
    Rest[LineSimplify[
      Join[{midpoints}[[
        maxdistanceposition ;; -1]], {p2}], \[Epsilon]]]
    ]
   ]
  ]
LineSimplify[{p1 : {_, _}, p2 : {_, _}}, \[Epsilon]_] := {p1, p2}

となる。距離の最大値と最大値を与える点の番号を別々に求めているところがちょっと冗長だが仕方ない。

MathematicaのCountryDataから取り出せる日本の形状データを使ってこのアルゴリズムをテストしてみた例。

間引き前(2005点) 間引き後(212点) ε=0.1

間引いた後の点の数を指定できるようにするには

さて、Douglas-Peuckerのアルゴリズムは許容誤差εを与えて間引くので、間引いた後の点の数がいくらになるのかは指定できない。 しかし、GPSログデータの編集・加工ツールとして有名な轍(Wadachi)では、間引いた後の点の数を指定できる。

これを実現するには、以下のようなアルゴリズムになる。

  1. 始点と終点を結ぶ線分と各点の距離を求める。
  2. 距離が最大の点Pとその距離の情報をリスト(優先度付きキュー)に加える
  3. リストからもっとも距離の大きい要素を取り出し、点Pを残す点として登録し、残す点の個数をカウントする
  4. 点Pを新たな終点として、最大距離を与える点とその距離をリストに加える
  5. 点Pを新たな始点として、最大距離を与える点とその距離をリストに加える
  6. リストが空になるか、点の個数が指定した数になるまで3以降を繰り返す。
  7. 残す点として登録したものだけを取り出す

というアルゴリズムになる。

GPXデータを簡略化ではこれを実装した。


Jun. 26 2014 追記:

Mathematica 8で実装してみた。簡単のため優先度付きキューではなく配列を毎回Sortしている。

SegmentPointDistance2 = 
  Compile[{{a, _Real, 1}, {b, _Real, 1}, {p, _Real, 1}}, 
   Module[{ax, ay, bx, by, px, py, x, y, t},
    ax = a[[1]]; ay = a[[2]];
    bx = b[[1]]; by = b[[2]];
    px = p[[1]]; py = p[[2]];
    t = (ax^2 + ay^2 + bx px - ax (bx + px) + by py - 
        ay (by + py))/(ax^2 + ay^2 - 2 ax bx + bx^2 - 2 ay by + 
        by^2);
    x = ax - (ax - bx) t;
    y = ay - (ay - by) t;
    Which[
     0 <= t <= 1, Sqrt[((x - px)^2 + (y - py)^2)],
     1 < t, Sqrt[((bx - px)^2 + (by - py)^2)],
     True, Sqrt[((ax - px)^2 + (ay - py)^2)]
     ]
    ]
   , CompilationTarget -> "C"
   , Parallelization -> True
   , RuntimeAttributes -> {Listable}
   ];
findFarthest[data_, pos1_, pos2_] /; pos2 - pos1 >= 2 := 
  Module[{dist, maxdist, maxpos},
   dist = 
    SegmentPointDistance2[data[[pos1]], data[[pos2]], 
     data[[pos1 + 1 ;; pos2 - 1]]];
   maxpos = Ordering[dist, -1][[1]];
   maxdist = dist[[maxpos]];
   {{maxdist, pos1, pos1 + maxpos, pos2}}
   ];
findFarthest[data_, pos1_, pos2_] := {};
lineSimplifyNStep[data_, {first_, rest___}] := 
 Module[{dist, pos1, maxpos, pos2},
  {dist, pos1, maxpos, pos2} = first;
  Sow[maxpos];
  [email protected][Join[
     {rest},
     findFarthest[data, pos1, maxpos],
     findFarthest[data, maxpos, pos2]
     ]]
  ];
LineSimplifyN[line : {{_, _} ..}, n_ /; n <= 0] := {};
LineSimplifyN[line : {{_, _} ..}, 1] := {First[line]};
LineSimplifyN[line : {{_, _} ..}, 2] := {First[line], Last[line]};
LineSimplifyN[line : {{_, _} ..}, n_] /; Length[line] <= n := line;
LineSimplifyN[line : {{_, _} ..}, n_Integer] := 
  line[[{1, Sequence @@ Sort[
       Reap[Nest[lineSimplifyNStep[line, #] &,
          findFarthest[line, 1, Length[line]], n - 2]][[2, 1]]
       ], Length[line]}]];

例(アメリカ合衆国)

line = CountryData["USA", "Polygon"][[1, 1]];
Graphics[{Polygon[line], Red, [email protected][line, 50]}]

元は4672点、赤線が50点に削減したもの。