earcutより高速: 370万の地図ポリゴンを489ミリ秒で三角形分割

earcutより高速: 370万の地図ポリゴンを489ミリ秒で三角形分割

当社の地図エンジンのポリゴン三角形分割器を、earcutと同じISCライセンスの下で、小さなライブラリearstripとして公開しました。これは1つのC++ヘッダーファイルと、earcutと同じ引数を受け取るnpm上のWebAssemblyビルドです。C++ではearcut.hppより14%高速、ブラウザではearcut.jsの1.4～2.1倍高速です。比較はほぼ最後にあります。この記事では、その方法について説明します。

GPUは三角形を塗りつぶすことしかできません。一方、地図は主にポリゴンで構成されています。建物、公園、駐車場、2つの路面電車の線路の間にある草地などです。以下は、当社のレンダラーが地図ファイルから読み込む、ワルシャワのŚródmieście地区の数ブロックです。

これらが塗りつぶされる前に、すべての輪郭は三角形に分割される必要があります。ポリゴンにカーソルを合わせるか、タップすると、それがどのように変化するかを確認できます。これらのブロックには331個のポリゴンと5,151個の頂点が含まれており、4,514個の三角形になります。

マゾフシェ県（ワルシャワ周辺地域）のオフライン地図のズームレベル13には、3,740,425個のポリゴンがあります。M4 Proの1コアでは、当社のエンジンの制約付きDelaunay三角形分割器がそれらすべてを処理するのに2.93秒かかります。今年の10月、これをear clippingに置き換えました。この記事を書くために、私はear clipperをアイデアごとに、1つのコミットにつき1つのアイデアで再構築し、それらの同じ370万ポリゴンで各コミットをタイミング測定しました。そのため、以下の各セクションの最後には数値があります。

Delaunay​

点の集合の三角形分割において、Delaunayとは、どの点も任意の三角形の3つの頂点を通る円の内側に入らない場合を指します。辺を共有する2つの三角形があれば、この規則がどのように機能するかを確認できます。点Dをドラッグしてください。

45°ABCD

最小角度45°; BDにすると37°になります

DがA、B、Cを通る破線の円の外側にある限り、共有辺はACです。そして辺はBDに反転し、元に戻ります。各反転は局所的な修正であり、どの辺も反転したがらない三角形分割はDelaunayです。同じ点を三角形分割するすべての方法の中で、最小角度が最大になります。オレンジ色の弧は2つの三角形の最小角度を示し、図の下の線はもう一方の対角線が与えるものを示しています。Dをどこに置いても、円の規則が選択する辺が勝ちます。

「制約付き」とは、ポリゴンの辺自体が、円の規則がそれらを排除しようとする場所でも、三角形分割の辺として維持されなければならないことを意味します。当社の三角形分割器は、頂点を1つずつ挿入します。各新しい点について、それまでに構築された三角形を横断して、点が落ちる三角形を見つけ、それを分割し、円の規則が再び成り立つまで新しい点の周りの辺を反転させます。また、ポリゴンの辺を強制的に挿入し、最後にポリゴンの内側で穴の外側にある三角形のみを保持します。

形状の良い三角形は、それらを介して何かを補間したり、それらでシミュレーションを実行したりする場合に重要です。平坦な塗りつぶしは気にしません。GPUは、細長い三角形も太い三角形も正しく塗りつぶします。そして、ほとんどの作業は、ほとんど必要としないポリゴンに費やされます。マゾフシェ県のポリゴンのほぼ半分、180万個は、頂点が4つ以下です。三角形はまったく切断する必要がなく、凸四角形は2つの三角形ですが、Delaunayコードはそれらに419ミリ秒を費やし、各三角形分割をゼロから構築します。

Ears​

ポリゴンの「耳」とは、1本の直線で切り取ることができる角のことです。連続する3つの頂点a、b、cを取ります。bの角が凸であり、ポリゴンの他の部分が三角形abcに食い込んでいない場合、その三角形は耳です。それを切り取ると、残りは頂点が1つ少ないポリゴンになります。1975年、Gary Meistersは、3つより多い頂点を持つ単純なポリゴンには、重なり合わない耳が少なくとも2つあることを証明しました。したがって、常に別の耳を切り取ることができ、n個の頂点を持つポリゴンはn-2個の三角形に分解されます。

上記のブロックにある建物の一つです。スライダーを動かして耳を1つずつ切り取るか、以下のようにしてください。

0 / 25

緑色の三角形が次に切り取られる耳です。アルゴリズムは1つのテスト、「この角は耳か？」を必要とします。その最も単純なバージョンは、角をチェックし、次に境界の他のすべての頂点をチェックします。crossは、すべてのジオメトリコードベースにある2x2の行列式です。aからbへの線の左側にcがある場合、正になります。したがって、反時計回りに実行される境界の場合、cross(a, b, c)が正であれば、それは凸角を意味します。

static double cross(double ax, double ay, double bx,                    double by, double cx, double cy){
    return (bx - ax) * (cy - ay) - (by - ay) * (cx - ax);
}

bool isEar(uint32_t v) const{
    const Node &b = _nodes[v];
    const Node &a = _nodes[b.prev], &c = _nodes[b.next];

    if (cross(a, b, c) <= 0)
        return false; // turns the wrong way, or not at all

    for (uint32_t w = c.next; w != b.prev;
         w = _nodes[w].next) {
        const Node &p = _nodes[w];
        if (cross(a.x, a.y, b.x, b.y, p.x, p.y) < 0 ||
            cross(b.x, b.y, c.x, c.y, p.x, p.y) < 0 ||
            cross(c.x, c.y, a.x, a.y, p.x, p.y) < 0)
            continue; // outside the triangle

        if (!samePos(p, a) && !samePos(p, b) &&
            !samePos(p, c))
            return false;
    }
    return true;
}

三角形の辺上の点は内側にあると見なされます。耳が境界の残りの部分に触れていた場合、それを切り取るとポリゴンが自己接触することになります。唯一の例外は、a、b、またはcの位置に正確にある点です。これは、輪郭が単一の頂点で自分自身と交差する場合に発生します。穴があると、それは一般的になります。

Holes​

ear clippingは1つの閉じた境界で機能しますが、マゾフシェ県のポリゴンの11,345個には、中庭や森の空き地のような穴があります。それらを処理する方法は、各穴を開いて、穴の頂点から境界の頂点への橋（両方向に1回ずつ歩く辺）で外側のリングに縫い付けることです。その後、境界は外側を回り、橋を渡り、反対方向に穴の周りを回って戻ります。2つの橋のコピーの間には領域がないため、ポリゴンは同じですが、現在は1つのリングになっています。

橋は何も横断してはなりません。David Eberlyのear clippingに関する解説には、それを見つけるための巧妙な方法があります。穴の最も右にある頂点mを取り、右にレイをキャストします。レイが最初に当たる辺、I、の端Pはさらに右にあります。境界の頂点が三角形m、I、Pの内側になければ、橋m-Pは安全です。穴をドラッグできます。

IPm

へこみの先端が三角形m、I、P内にあります。Pへの橋は境界を横断するため、へこみに向かいます。

三角形は空であり、橋はPに向かいます。へこみの先端が三角形内にあり、橋は代わりにへこみに向かいます。へこみのうち、Eberlyはレイに対して角度が最も近いものを選びます。穴は右端から順に橋渡しされるため、左側にある穴はすでに境界の一部になっている穴に橋渡しできます。

地図データは独自の2つのケースを追加します。当社の座標はタイルグリッド上の整数であるため、グリッド線に沿ったレイは常に実行されます。交差テストはこのような頂点を通過し、私の最初のバージョンはこのケースを正確に間違えました。また、マッパーはしばしば中庭を描画しますが、それはある時点で単一の頂点になります。そのような穴には橋は必要ありません。2つのリングは共有頂点で結合されます。

The naive clipper (−1,812 ms)​

穴が橋渡しされた後、最初の完全なバージョンは境界を歩き回り、上記のisEarで角をテストし、見つかったすべての耳を切り取り、各三角形を3つのインデックスとして書き出します。マゾフシェ県のポリゴンでは、Delaunayの2.6倍速い1.12秒かかります。チャートはポリゴンのサイズごとに分割されます。バーの一部にカーソルを合わせるとその数を確認できます。行をクリックしてください。

Delaunay2,931 ms

1. Naive ears1,119 ms
2. Only dents882 ms
3. Convex shortcut880 ms
4. Zig-zag strips762 ms
5. Buckets697 ms
6. No branches519 ms
7. Marked hits489 ms

頂点数4つまで5～8頂点9～32頂点33～256頂点257以上

最後のサイズクラス以外はすべて短くなっています。257以上の頂点を持つ1,123個の最大のポリゴンは、クリップに199ミリ秒、Delaunayでは90ミリ秒かかります。各耳のテストでは数百の頂点を歩き、ポリゴンは数百のテストを必要とするため、コストはサイズに比例して増加します。この記事の残りは、earテストをより安価にし、頻度を減らす方法についてです。

Test only the dents (−237 ms)​

どの頂点が実際に凸角が耳になるのを妨げる可能性があるかを見てみましょう。たとえば、いくつかの頂点がその三角形内にあり、角自体に最も近いものを取り上げます。その2つの辺は三角形の辺を横切ることはできないため、両方とも角から離れて戻り、ポリゴンはその角の側からその頂点の周りをラップします。これにより、それは凹角（へこみ）になります。したがって、三角形内に何かがあれば、へこみもその内側にあり、テストはへこみのみを調べる必要があります。

エンジンはそれらのリストを保持します。へこみが隣の角を切り取って凸にするとリストから外れ、凸角はへこみになることはないため、リストは縮小するだけです。テストは次のようになります。

bool isEar(uint32_t v) const{
    const Node &b = _nodes[v];
    if (b.reflex)
        return false;

    const Node &a = _nodes[b.prev], &c = _nodes[b.next];

    for (uint32_t i = 0; i < _live; ++i) {
        const Item &p = _items[i]; // a dent
        if (p.node == b.prev || p.node == b.next)
            continue;

        if (cross(a.x, a.y, b.x, b.y, p.x, p.y) < 0 ||
            cross(b.x, b.y, c.x, c.y, p.x, p.y) < 0 ||
            cross(c.x, c.y, a.x, a.y, p.x, p.y) < 0)
            continue;

        if ((p.x != a.x || p.y != a.y) &&
            (p.x != b.x || p.y != b.y) &&
            (p.x != c.x || p.y != c.y))
            return false;
    }
    return true;
}

2810万回のearテストは、3億5300万点ではなく1億8500万点を見るようになり、合計は減少します。これは期待されるよりも少ないです。なぜなら、輪郭の構築には多くのへこみがあるからです。上記のアーチのある建物には、27個の頂点のうち15個のへこみがあります。また、最小のポリゴンは遅くなり、34ミリ秒から54ミリ秒になります。三角形や四角形の場合、角をへこみとそれ以外にソートするコストは、節約されるウォークよりも高くなります。

Skip tests on convex polygons (−2 ms)​

へこみのないポリゴンは凸であり、凸ポリゴンのすべての角は耳であるため、単一のテストなしで切り取ることができます。それをチェックするには注意が必要です。「すべての角が左に曲がる」だけでは十分ではありません。

左に曲がる角: 5/5。x方向の変化: 4。凸ではない。

すべての5つの角で左に曲がります。なぜなら、中心を2回回るからです。エンジンはもう1つの条件を追加します。凸のリングを歩くと、辺は下部で右に、上部で左に進むため、x方向の変化は左端と右端の点で正確に2回変化します。スターは円で囲まれた角で4回方向を変え、通常のearテストに進みます。両方のカウントは、へこみをマークするのと同じループから取得されるため、チェックのコストはほぼゼロです。三角形と四角形は、そのループの前に独自のショートカットを取得します。

ショートカットは2810万回のearテストから440万回を削除しますが、それらは最も安価なテストでした。へこみのないリングにはスキャンするものがありません。それ自体では合計にほとんど影響しませんが、頂点が4つ以下のポリゴンは、へこみがコストをかけた一部を取り戻し、54ミリ秒から48ミリ秒になります。次のステップでは、凸リングが単一のストリップになるため、より重要になります。

Write strips (−118 ms)​

これまでのところ、すべての三角形は3つのインデックスを必要とします。GPUは三角形ストリップも描画できます。各新しいインデックスは、それと前の2つの頂点からなる三角形を追加します。凸リングは、それをジグザグに横切ると1つのストリップを作成します。

01234567

0 1 7 2 6 3 5 4

ear clipperもストリップを書き込むことができます。耳を切り取ると、その2つの隣接頂点の角がより鋭くなるため、それらは次の耳の候補として適しています。それらのいずれかを切り取ると、ストリップが1つのインデックスだけ拡張されます。エンジンはそれらを最初に試し、耳が続く限り両端を交互に使用します。コードでは、切り取られた頂点の隣接頂点はaとbです。

// The next ear alternates between the ends of the cut.
while (true) {
    if (sideA) {
        const uint32_t z = n[a].prev;
        _strips.push_back(n[z].index);
        unlink(a);
        a = z;
    } else {
        const uint32_t z = n[b].next;
        _strips.push_back(n[z].index);
        unlink(b);
        b = z;
    }
    count--;
    sideA = !sideA;
    if (count < 3 || !isEar(sideA ? a : b))
        break;
}

次の角が耳でない場合、ストリップは終了し、ウォークは停止した場所から続行されます。2150万個の三角形は、6450万個の代わりに3580万個のインデックスを必要とします。earテストは減りません。ナイーブウォークはすでにカットの後の次の角で続行されており、今では各ストリップも失敗するテストで終了します。2370万回ではなく2500万回です。それでも、合計は減少し、最大のポリゴンは146ミリ秒から103ミリ秒に低下します。

Sort the dents into buckets (−65 ms)​

森や湖には数千の頂点と数百のへこみがある場合がありますが、各earテストはそれらすべてを読み取ります。ポリゴンに32個以上のへこみがある場合、エンジンはそれらを境界ボックスの長い辺に沿ってバケットにソートします。1バケットあたり約2個のへこみです。earテストは、三角形がまたがるバケットのみを調べます。へこみが凸になり離れると、インデックスは半分がなくなった時点で再構築されます。

スキャンされたポイントは1億8100万から1億100万に減少し、そのほとんどは大きなポリゴンで発生します。257以上の頂点を持つポリゴンは、103ミリ秒から57ミリ秒になります。合計:。

Count without branches (−178 ms)​

へこみテストのループは、三角形内の最初のへこみで終了します。ほとんどの場合、へこみはありません。テストはすべてのへこみを読み取り、はいと答えます。早期終了は何も得られず、ポイントごとの3つの比較はCPUが予測する必要のある分岐です。エンジンは最初に、分岐なしで三角形内のへこみを数えるだけのパスを実行します。

const double abx = b.x - a.x, aby = b.y - a.y;
const double bcx = c.x - b.x, bcy = c.y - b.y;
const double cax = a.x - c.x, cay = a.y - c.y;

uint32_t hits = 0;
for (uint32_t i = begin; i < end; ++i) {
    const Item &p = items[i];
    const double side = std::min(
        std::min(abx * (p.y - a.y) - aby * (p.x - a.x),
                 bcx * (p.y - b.y) - bcy * (p.x - b.x)),
        cax * (p.y - c.y) - cay * (p.x - c.x));
    hits += (uint32_t)(side >= 0) & (uint32_t)(p.node != ia) 
          & (uint32_t)(p.node != ic);
}

if (hits == 0)
    return true;

3つの行列式の最小値が非負になるのは、点が閉じた三角形の内側にある場合のみです。カウントがゼロでない場合にのみ、前の正確なループが実行され、角にあるへこみが除外されます。バケットバージョンは、64個のポイントを一度にビットマスクにまとめて、同じことを行います。合計は減少し、中間サイズが最も多く得られます。33～256個の頂点を持つポリゴンは、327ミリ秒から216ミリ秒になります。

Test only the marked dents (−30 ms)​

カウントがゼロでない場合、正確なループはすべてのへこみを再度読み取りますが、最初のパスではどのへこみが内側にあるかを確認しています。バケットバージョンはビットマスクにそれを保持します。最大32個のへこみを持つ短いリストも同様に行うことができます。最初のパスはカウントの代わりにヒットごとに1ビットを設定し、正確なテストはそのへこみのみを調べます。

    hits |= ((uint32_t)(side >= 0) & (uint32_t)(p.node != ia) 
             & (uint32_t)(p.node != ic)) << (i - begin);
}

while (hits) {
    const Item &p = items[begin + lowestBit(hits)];
    hits &= hits - 1;
    // the exact test from before
}

読み取られるポイントは1億2800万から1億2000万に減少します。これは最後のステップであり、公開されているearstripのコードです。同じポリゴンでDelaunayよりも6.0倍少ない時間で完了します。

When there are no ears​

Meistersの定理は単純なポリゴンに関するものです。地図データはタイルグリッドに丸められ、丸めによって単純なポリゴンが壊れることがあります。これは、上記のブロックから1キロメートル未満の鉄道駅であるワルシャワ中央駅の輪郭で、さまざまなサイズのグリッドに丸められたものです。オレンジ色のリングは、隣接する頂点と同じ位置に着地する頂点をマークしています。

グリッドステップ73cm。154個の頂点のうち68個が隣接する頂点と同じ位置に着地します。84個の三角形。

タイルあたり18cmのグリッド（地図ファイルで使用されるグリッド）では、グリッドステップは18cmです。それでも、154個の頂点のうち31個が隣接する頂点とマージされます。なぜなら、マッパーはしばしば数センチメートル離れた場所に点を配置し、輪郭はそのまま残るからです。頂点のほぼ半分がマージされ、クリッパーはそれをまだカットします。丸められた輪郭は、耳がまったく残らないように自己接触し、10ビットでは再び機能します。

クリッパーが境界を一周しても耳が見つからない場合、面積のない角を削除してもう一度歩きます。それでも失敗した場合、ギブアップし、エンジンはDelaunayコードでそのポリゴンを三角形分割します。これはまだ存在します。マゾフシェデータでは、76個のポリゴンがそこに到達します。それらは最後の変更の前に構築された地図から来ています。当社の地図ジェネレーターは現在、GEOSのスナップ丸めを使用してすべてのポリゴンを最終グリッドに丸め、ear clippingでカットできることを確認しています。マゾフシェ地図の新しいビルドでは273個のポリゴンが修復され、ear clippingはファイル内のすべてのポリゴンをカットできます。

Delaunayにはもう1つの仕事があります。各穴には独自の橋が必要です。それを見つけるには境界全体をスキャンする必要があるため、数千の穴を持つポリゴンは数千回のスキャンに費用がかかります。当社の世界地図には11,675個の穴を持つポリゴンがあります。ear clippingはそれに606ミリ秒、Delaunayは17ミリ秒かかります。エンジンは現在、128個以上の穴を持つすべてのポリゴンをDelaunayに送信しています。マゾフシェデータでは、これは2つのポリゴンです。

ライブラリはこれらのフォールバックの両方を内部に保持します。earstripは、耳がなくなった場合、それらが残した三角形が間違った領域をカバーする場合、および128個以上のリングがある場合に、リングエッジの制約付きDelaunay三角形分割でポリゴンをカバーします。C++のrun()とJavaScriptのlastStatus()は、ポリゴンがどちらの方法で処理されたかを示します。当社の世界地図タイルの56,472個のポリゴン（1つのポリゴンに最大11,675個の穴）では、earstripは459ミリ秒、earcut.hppは2,892ミリ秒かかります。耳が正確にカットできないズームレベル13のポリゴンでは、earcutの三角形は、平均31%のピクセルでリングの偶奇塗りつぶしと異なります。earstripは平均0.15%です。

In the browser​

新しいクリッパーがアプリに組み込まれる前に、MapLibreとMapbox GLのバックエンドにあるear clipper、earcutのC++バージョンearcut.hpp 3.2.4と比較しました。ネイティブでは、両者は近く、それで終わりになる可能性もありました。しかし、私たちはWebAssemblyにコンパイルされた同じC++エンジンを実行するWebレンダラーに取り組んでいるため、さらに進んでブラウザ側も比較しました。両方のクリッパーをWebAssemblyにコンパイルし、今日Webマップが実行しているJavaScriptオリジナルであるearcut.js 3.2.4と比較しました。そのために、クリッパーをエンジンからearstripに抽出しました。当社のテストセットのすべての580万個のポリゴンで、エンジンと同じストリップを正確に書き出します。そのリポジトリには、このセクションのすべての数値の背後にあるポリゴンとコードがあります。ネイティブベンチマークは上記のチャートの背後にあるもので、ネイティブ行はearcut.hppに対して実行された最後のステップです。すべてのWebAssemblyビルドは、ズームレベル13と9の440万個のポリゴンすべてで、ネイティブバージョンとまったく同じ三角形を生成します。ズーム13の374万個のポリゴンに対するミリ秒、7回のベスト、1スレッド、NodeのV8およびBunのJavaScriptCore：

ネイティブではearcut.hppはearstripより14%遅いです。WebAssemblyにコンパイルされ、V8で実行されると、earstripは23%低下し、earcut.hppは67%低下します。そのほとんどは、earを探すメインループで発生しています。アロケーターではありません。mimallocとemmallocでのビルドは同じ時間をもたらします。さらに深く掘り下げていません。V8では、WebAssemblyのearcut.hppは、プレーンなearcut.jsよりもわずか19%高速になり、JavaScriptからの各ポリゴンへの個別の呼び出しでは10%速くなります。

パッケージ行はWebページが取得するものであり、そこでのearcut列はearcut.jsです。earcutを置き換える呼び出しは、座標をJavaScript数値として受け取り、毎回新しい三角形インデックスの配列を返します。V8ではearcut.jsの1.4倍、JavaScriptCoreでは1.6倍高速です。モジュール内のループに対する低下の約半分は、入力座標のコピーとストリップをリストに変換することに費やされ、残りはポリゴンごとに1回WebAssemblyに呼び出すことに費やされます。直接呼び出しはそれらをすべてスキップします。ポリゴンはモジュールのメモリに直接入り、WebGL 2またはWebGPU用の単一のインデックスバッファとして戻ってきます。earcut.jsの1.7倍から2.1倍高速です。

上記のアイデアの一部はearcutにも含まれています。earcutは、耳を妨げるのはへこみだけであるという事実に依存していますが、バウンディングボックスと三角形テストを通過したすべての頂点に対してそれを最後にチェックします。へこみのリストを保持するのではなく。80頂点を超えると、それらをz順にリンクし、三角形のバウンディングボックスがカバーするその曲線のストレッチのみをスキャンします。これはバケットの役割を果たします。earcutには、凸ポリゴンのショートカット、ストリップ、および分岐なしのカウントはありません。三角形ごとに3つのインデックスを書き込み、内側の最初の点で停止します。

earstripはGitHubとnpmで入手できます。この記事のポリゴンと、独自のコンピューターでearcutに対してタイミングを測定するためのコードが含まれています。