ハーフエッジデータ構造を理解する

# はじめに

本記事ではハーフエッジデータ構造をなるべく根本から理解しつつ、なるべく楽してハーフエッジデータ構造を実装します。

前提・レギュレーション

- 入力は境界付き多様体三角形メッシュ
- 入力メッシュの全ての面は表向き（本記事では頂点順が時計回り）であると仮定する
- C++20 を使う
- Halfedge クラス・構造体を作らない
- パフォーマンスは重視しない

要約すると、そこまで複雑ではない形状に対して使える簡易的なデータ構造を作るということです。  

また、タイトルではわかりやすさの観点からハーフエッジという名称を使っていますが、以降の本記事内ではハーフエッジを<ruby>半辺<rt>はんぺん</rt></ruby>と訳して記述します。


# メッシュについて
半辺に触れる前に、まずはその前提から理解していきます。

本記事におけるメッシュとは複数の三角形要素によって表された形状のことです。次の図は平面におけるメッシュの例です。

![mesh_image11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/18561d1a-d143-e332-d949-c5a57c839539.png)

これはメッシュの一部を拡大しているためわかりづらいですが、実際のメッシュは三角形要素の集合によって何らかの形状を表現していると思ってください。

メッシュは頂点集合 `vertices` と面集合 `faces` を持っています。各頂点は座標を持ち、各面は面を構築する 3 つの頂点番号を持っています。  
平面のメッシュの要素を構造体で作ると大体次のようになります。

``` c++
struct Vertex{
    double x;
    double y;
};

struct Face{
    int v0;
    int v1;
    int v2;
};

struct SimpleMesh{
    std::vector<Vertex> vertices;
    std::vector<Face> faces;
}
```

ところで、あるメッシュに対して、データ量削減のために面数をもっと少なくしたいとか、角ばっている部分を滑らかにしたいとか、そのようなメッシュの変形の操作を考えることはよくあります。あるいは変形ではなく、単にメッシュの面や頂点の隣接・接続情報を得たいといったこともよくあります。

このような操作を実装するために便利なのが **半辺（ハーフエッジ）データ構造** です。


# 半辺について
先に半辺データ構造の理論と構築について説明します。半辺を使った具体的な操作や半辺の用途は後で記述します。

そもそも **半辺（ハーフエッジ）** とは、雑に言えば 1 つの三角形面につき 3 つある有向辺のことです。  
しかし、通常の辺も言ってしまえば 1 つの三角形に 3 つあるわけで、言葉だけでは半辺と通常の辺の違いはよくわかりません。メッシュを図にして考えると半辺と通常の辺の違いがわかりやすいです。

![image22.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/1853eef9-b8fa-2e93-8d9a-efffc1283d6a.png)

例えば、上の図のように 2 つの面からなるメッシュを考えます。
このメッシュの辺の数は見てわかるように 5 つですが、半辺の数は 6 つとして考えます。
1 つの面につき 3 つの半辺があると考えるためです。  

![image23.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/2373ae69-95ce-c1c7-9aea-ab47c6a351bb.png)

半辺を矢印で示したものが上の図になります。 1 つの面につき矢印、つまり半辺が 3 本あり、合計で 6 本の半辺があるのがわかります。通常の辺の数は 5 本ですが、これは図の青と赤の半辺をまとめて 1 つの辺と数えているからです。

一旦、この図からわかる半辺の要所を箇条書きでまとめます。  

- 1 つの面に 3 つの半辺がある
- 半辺は向きを持つ
- 面の 3 つの半辺は時計回りになる
- 境界辺は 1 つ、内部辺は 2 つの半辺を持つ
- 内部辺に属する半辺は、始点と終点を逆にした **反対辺（opposite）** を持つ（図の青と赤の半辺）。境界辺に属する半辺は反対辺を持たない
- ある半辺の反対辺は、その半辺とは異なる面に属する

ここで、**境界辺** と **内部辺** という概念が登場しました。境界辺とは 1 つの面のみに接続する辺で、内部辺は 2 つの面に接続する辺のことです。上の図では、真ん中の縦線が内部辺で、そのほかの 4 つの辺が境界辺となります。

内部辺は 2 つの半辺を持ち、これが半辺の名前の由来でもあります。上の図にある内部辺は、青と赤の 2 つの半辺を持っています。また、青の半辺の始点と終点を反転させると赤の半辺が得られるため、青の半辺と赤の半辺は反対辺の関係にあると言えます。一方、境界辺にあたる黒の半辺は反対辺がありません。この反対辺という情報が半辺データ構造では重要です。

半辺データ構造では、全ての半辺の反対辺を記憶しておきます。問題は反対辺をどう求めるかという点です。


# 半辺データ構造の構築
全ての半辺について反対辺を探す処理を作りますが、その前に半辺の番号付けの規則を決めておきます。

![mesh_image21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/243b275e-8036-7da1-15f2-64fcf4252530.png)

上の図は面 $f_i$ の半辺の番号付けです。  
面 $f_i$ には順に $(v_n, v_m )$, $(v_m, v_l )$, $(v_l, v_n )$ の計 3 つの半辺が存在しますが、これらの半辺の番号をそれぞれ $3i$, $3i+1$, $3i+2$ とします。こうしておくと半辺を扱いやすいです。

それでは、各半辺に対して反対辺を探します。次のようにします。

1. 頂点添え字 2 つのペアをキーとするマップ $M$ を用意する。
2. メッシュの各面の各半辺 $h_i$ （始点を $v_b$, 終点を $v_e$ とする）に対して次の処理を行う。
3. $M[(v_e, v_b )]$ に既に値が入っているかどうかで分岐する。
    - 入っている場合： $M[(v_b, v_e )] = -(M[(v_e, v_b )] + 1)$, $M[(v_e, v_b )] = h_i$ とする。
    - 入っていない場合： $M[(v_b, v_e )] = -(h_i + 1)$ とする。 

コードはこのようになります。

``` c++
using halfedge_type = std::pair<int, int>;
std::map<halfedge_type, int> opps;

// 各半辺の反対辺を探して opps に記録する処理
for (int i = 0; i < (int)faces.size(); ++i) {
    for (int j = 0; j < 3; ++j) {
        halfedge_type edge = get_halfedge(i, j);
        halfedge_type sw_edge = swapped_halfedge(edge);
        if (opps.contains(sw_edge)) {
            // 既に反対辺が含まれているならばお互いにその半辺番号を登録
            opps[edge] = -(opps[sw_edge] + 1);
            opps[sw_edge] = i * 3 + j;
        }
        else {
            // まだその半辺の反対辺が無い場合、自分の半辺番号を負にして登録しておく
            opps[edge] = -(i * 3 + j + 1);
        }
    }
}
```

このように、全ての半辺ごとに反対辺がマップに登録されているか確認し、その都度半辺の番号をマップに格納していけばよいです。
`get_halfedge()` と `swapped_halfedge()` の実装は以下です。どちらも難しいものではないです。

``` c++
// face_idx 番目の面の halfedge_idx 番目の半辺を取得する
halfedge_type get_halfedge(int const face_idx, int const halfedge_idx){
    halfedge_type result;
    switch (halfedge_idx) {
    case 0:
        // v0 -> v1 の半辺
        result.first = faces[face_idx].v0;
        result.second = faces[face_idx].v1;
        break;
    case 1:
        // v1 -> v2 の半辺
        result.first = faces[face_idx].v1;
        result.second = faces[face_idx].v2;
        break;
    case 2:
        // v2 -> v0 の半辺
        result.first = faces[face_idx].v2;
        result.second = faces[face_idx].v0;
        break;
    }
    return result;
}

// 辺の始点と終点をひっくり返す
halfedge_type swapped_halfedge(halfedge_type const& halfedge){
    halfedge_type result = halfedge;
    std::swap(result.first, result.second);
    return result;
}
```

このやり方の肝は、反対辺がマップにまだ登録されていない場合自身の半辺の番号を負にして登録することです。こうしておくとことで、反対辺が存在しない半辺の `opps` の値は自動的に負になるので、境界辺の判定が楽になります。境界辺の判定はよく登場するため、楽であることに越したことはありません。

これにより各半辺の反対辺が `opps` というマップに記録され、半辺データ構造が大体出来上がりました。


# 半辺の操作
実はこの時点で、既に次のようなメッシュ上の操作が高速に行えるようになっています。

- ある面からその面に属する半辺を得る（face_to_halfedge | f_to_h）
- 半辺番号からその半辺が属する面番号を得る（halfedge_to_face | h_to_f）
- 半辺から反対辺にあたる半辺番号を得る（halfedge_to_halfedge | h_opps）
- 半辺からその始点・終点の頂点番号を得る（halfedge_to_vertex | h_src, h_dst）

f_to_h は上述した `get_halfedge()` と同じです。  
h_to_f は、半辺の番号付けの規則より、半辺の番号を 3 で割れば面番号を得られます。  
h_opps はマップ `opps` を使えばよいです。  
h_src と h_dst は、半辺の始点 `vb` または終点 `ve` を参照すればよいです。  

また、次の操作も可能です。

- 半辺番号から半辺（始点・終点）を得る（h_make）
- ある半辺が境界辺かを判定する（is_boundary）
- ある半辺の半辺番号を得る（h_idx）

半辺番号を x とすると、 `get_halfedge(x/3, x%3)` で始点・終点を持つ半辺を得られます。  
半辺の始点を `vb`, 終点を `ve` とすると、 `opps[(vb, ve)] < 0` が true ならその半辺は境界辺です。  
また、 `opps[(ve, vb)]` でその半辺の番号を得られます。半辺が境界辺である場合は、 `-(opps[(vb, ve)] + 1)` で半辺番号を得られます。

実際に半辺の操作の例を見てみます。

![image24.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/762ea790-d0cf-5351-beed-bd7dd823f1e8.png)

この図は半辺データ構造の具体的な例です。あるメッシュから 2 つの面を取り出して考えた状態だと思ってください。  
それぞれの面・頂点・半辺に通し番号がついています。頂点と面の番号は適当です。

面 $f_3$ について、 11 番目の半辺を跨いだ先にある面の番号を計算する操作を考えます。  
まず h_opps で 11 番目の半辺の反対辺、つまり 21 番目の半辺を取得します。
次に h_to_f で、 21 / 3 = 7 を計算し、これで目標としていた $f_7$ を得られます。  

半辺を利用することで、隣の面の番号を得ることができました。この処理を繰り返せば、ある面に隣接する面を全て列挙することもできます。このような局所的なメッシュ操作は半辺データ構造の得意分野です。

ここで、新たな半辺の機能を追加します。

- ある半辺について、その面内における次の半辺を得る（halfedge_to_halfedge | h_next）
- ある半辺について、その面内における前の半辺を得る（halfedge_to_halfedge | h_prev）
- ある半辺からその半辺の対角頂点の番号を得る（halfedge_to_vertex | h_to_v_opps）

![mesh_image22.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/f5be95c2-1468-2297-65d8-d9e7b82fa6c3.png)

上の図は、よくある半辺の解説図です。今注目している半辺が図中の赤い矢印であり、その半辺に対する h_next や h_prev, opposite や v_opps などの位置関係が記されています。赤い半辺一つから、黒い文字で書かれている情報にすべてアクセスできます。

h_next を強引に実装し、 h_prev と h_to_v_opps を h_next を使って実装します。

``` cpp
halfedge_type h_next(halfedge_type const & halfedge) {
    int const face_idx = h_to_f(h_idx(halfedge));
    if (h_dst(halfedge) == faces[face_idx].v0) {
        return f_to_h(face_idx, 0);
    }
    else if (h_dst(halfedge) == faces[face_idx].v1) {
        return f_to_h(face_idx, 1);
    }
    else {
        return f_to_h(face_idx, 2);
    }
}

halfedge_type h_prev(halfedge_type const& halfedge) {
    return h_next(h_next(halfedge));
}

int h_to_v_opps(halfedge_type const halfedge) {
    return h_dst(h_next(halfedge));
};
```

三角形メッシュでは一つ前の半辺は二つ後の半辺と等しいので、 h_prev は h_next 二回で実装できます。  
h_to_v_opps は、頂点の位置的に一つ後の半辺の終点を見ればよいので、 h_next 一回で実装できます。

h_next と h_prev があれば半辺を介した面の操作は更にやりやすくなります。これらの操作は実際に後で使います。

以上により、楽してメッシュの半辺データ構造を構築できました。次に具体的なメッシュ操作の中で半辺データ構造がどのように利用されるかを見ていきます。


# メッシュの基本操作
メッシュにおける幾何処理の 3 つの基本操作というものがあります。

- 辺縮約（edge collapse）
- 対角変形（edge flip）
- 辺細分（edge split）

これらを実装して半辺データ構造の切れ味を確かめたいところですが、これらの操作は実装するのは意外と大変です。ここでは辺縮約と対角変形を省き、辺細分のみを実装してみて、大体どういう感じの流れになるのかを確かめます。


## 辺細分について
辺細分は、メッシュへ新たに頂点を 1 つ追加することにより、 1 つの辺を 2 つの辺に細分します。頂点を追加する位置は辺上とは限りません。

![mesh_image32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/f9c10368-b278-581b-3669-441b3a6ef45f.png)

上の図では、赤い辺が辺細分の対象であり、新頂点は $\boldsymbol{v}_{new}$ となっています。  
辺細分により、頂点や辺だけでなく、面の数も増えていることがわかります。


## 辺細分の実装
実装で難しいのは半辺の反対辺の正しい再構築です。まずは具体的に面を見てみます。

![mesh_image34.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/60e1ca4e-67d5-cdb2-067f-31d55ce35df4.png)

この図は、左上が辺細分前のメッシュの状態、右側が辺細分後のメッシュの状態を示しています。赤い辺は辺細分の対象となった辺です。

$f_a$, $f_b$ は既存の面です。 $f_{new0}$, $f_{new1}$ が新たに発生する面です。辺細分後のこれらの面は、それぞれ次のような頂点の構成になります。

- $f_a$ : $(v_0, v_1, v_{new})$
- $f_b$ : $(v_2, v_0, v_{new})$
- $f_{new0}$ : $(v_1, v_3, v_{new})$
- $f_{new1}$ : $(v_3, v_2, v_{new})$

また、細分前の二つの面の外側に位置する半辺 $h_0$, $h_1$, $h_2$, $h_3$ は細分後も残ります。これらの値も後で使うため、細分をする前に記憶しておきます。ここで h_next や h_prev を使っていきます。

この図に従い、まずは半辺を介していくつかの情報を取得し、その後実際に頂点を追加し面の構成を更新していきます。

``` cpp
void edge_split(halfedge_type const& halfedge, Vertex const pos) {
    int fa, fb, fnew0, fnew1;
    int h0, h1, h2, h3;
    int v0, v1, v2, v3, vnew;

    // 面番号を取得
    fa = h_to_f(halfedge);
    fb = h_to_f(h_opps(halfedge));
    fnew0 = (int)faces.size();
    fnew1 = (int)faces.size() + 1;

    // 半辺の番号を取得
    h0 = h_opps(h_prev(halfedge));
    h1 = h_opps(h_next(halfedge));
    h2 = h_opps(h_prev(swapped_halfedge(halfedge)));
    h3 = h_opps(h_next(swapped_halfedge(halfedge)));

    // 頂点番号を取得
    v0 = h_dst(halfedge);
    v1 = h_to_v_opps(halfedge);
    v2 = h_to_v_opps(swapped_halfedge(halfedge));
    v3 = h_src(halfedge);
    vnew = (int)vertices.size();

    // 頂点を追加
    vertices.emplace_back(pos);

    // 面の更新
    // fa
    faces[fa].v0 = v0;
    faces[fa].v1 = v1;
    faces[fa].v2 = vnew;

    // fb
    faces[fb].v0 = v2;
    faces[fb].v1 = v0;
    faces[fb].v2 = vnew;

    // fnew0
    faces.emplace_back(v1, v3, vnew);

    // fnew1
    faces.emplace_back(v3, v2, vnew);

    // 後半に続く
}
```

次に、更新された面の頂点順序に従い、半辺の反対辺を更新していきます。 $f_{new0}$, $f_{new1}$ の新たに発生する半辺はもちろん、頂点の順序が変わる可能性のある $f_a$, $f_b$ の半辺の反対辺も更新する必要があります。

![mesh_image35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/b322dbd0-a7d4-bd40-86a4-4b6be498325c.png)

この図は反対辺の更新が必要な 16 本の半辺を示しています。青で書かれた数字は、その半辺の面内における相対的な番号付けです。

この番号付けを元に、半辺の反対辺を更新していきます。

```cpp
void edge_split(halfedge_type const& halfedge, Vertex const pos) {
    // 前半は省略

    // 半辺の反対辺を更新
    // fa
    opps[{v0, v1}] = h1 < 0 ? -(fa * 3 + 0 + 1) : h1;
    opps[{v1, vnew}] = fnew0 * 3 + 2;
    opps[{vnew, v0}] = fb * 3 + 1;

    // fb
    opps[{v2, v0}] = h2 < 0 ? -(fb * 3 + 0 + 1) : h2;
    opps[{v0, vnew}] = fa * 3 + 2;
    opps[{vnew, v2}] = fnew1 * 3 + 1;

    // fnew0
    opps[{v1, v3}] = h0 < 0 ? -(fnew0 * 3 + 0 + 1) : h0;
    opps[{v3, vnew}] = fnew1 * 3 + 2;
    opps[{vnew, v1}] = fa * 3 + 1;

    // fnew1
    opps[{v3, v2}] = h3 < 0 ? -(fnew1 * 3 + 0 + 1) : h3;
    opps[{v2, vnew}] = fb * 3 + 2;
    opps[{vnew, v3}] = fnew0 * 3 + 1;

    // 外周
    if (opps.contains({ v1, v0 })) {
        opps[{v1, v0}] = fa * 3 + 0;
    }

    if (opps.contains({ v0, v2 })) {
        opps[{v0, v2}] = fb * 3 + 0;
    }

    if (opps.contains({ v2, v3 })) {
        opps[{v2, v3}] = fnew1 * 3 + 0;
    }

    if (opps.contains({ v3, v1 })) {
        opps[{v3, v1}] = fnew0 * 3 + 0;
    }
}
```

地道な作業ですが、これで完成です。  
基本的には図を見ながら反対辺の対応を作っていきますが、 $h_0$, $h_1$, $h_2$, $h_3$ にあたる半辺 $( 
v_1, v_0)$, $(v_0, v_2)$, $(v_2, v_3)$, $(v_3, v_1)$ は注意が必要です。これらはその辺が境界辺である可能性があり、つまりこれらの半辺は存在しない可能性があります。境界辺の半辺の反対辺には自身の半辺番号を負にした値を入れていたので、もし境界辺であれば負の半辺番号を代入しなおす必要があります。コード中の三項演算子の部分がこの処理にあたります。  
また、存在しない半辺の反対辺を更新しないよう、 if 文で条件分岐を書いています。

これで辺細分した際の半辺情報の再構築は上手くできました。  
ここまで実装した辺細分は、内部辺についての辺細分です。面倒なことに、境界辺を細分した際の動きも作らなければなりません。

![mesh_image36.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/bd1e23d8-198a-afd5-355c-e696e5012b91.png)

と言っても、境界辺の辺細分は内部辺のときと比べ新たに発生する面が 1 つだけなので、面や反対辺の変更点も少ないです。上の図は境界辺の辺細分を表しています。

大体内部辺のときと同じように実装できるので、投げやりですが頑張って作ります。


## 辺細分のデモ
半辺のナンバリングが正しくできているか確かめるため、実際のメッシュ上で辺細分の動きを見てみます。

![mesh_image41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/74e5cc35-0c0b-f51b-acbb-10b6e12e5bac.png)

これは正六角形のメッシュの図です。数字は半辺の番号です。青は内部辺、緑は境界辺を示しています。  
1 番の半辺にあたる辺を辺細分してみます。

![mesh_image42.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/b5a0bbf0-4de3-0548-f62e-c5447bc69406.png)

こんな感じになりました。 2 や 8 の半辺の反対辺が変わっていることが確認できます。  
次に 4 の半辺を辺細分します。

![mesh_image43.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2984217/783e94ef-fa21-4096-8dc3-ff990e741d1d.png)

こんな感じになります。 1 の半辺の反対辺が 3 になっており、新たに発生した 24 が境界辺になりました。

この調子でひたすら辺細分し続ける動画です。

<iframe width="560" height="315" src="https://www.youtube.com/embed/xgeK3quesV0?si=L6XVvPuYv2lGba2v" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture; web-share" allowfullscreen></iframe>

半辺データ構造によって正しく面の接続情報を更新しているため、繰り返して辺細分の操作を行うことができます。新頂点の位置を辺の中点以外にしたり、面積の大きい面を優先して細分するなど、いろいろなバリエーションも考えられます。


## ソースコード
半辺データ構造の構築と辺細分の例のコードを載せています。

<details><summary>ソースコード</summary><div>

```cpp
struct Vertex {
	double x;
	double y;
};

struct Face {
	int v0;
	int v1;
	int v2;
};

struct SimpleMesh {
	using halfedge_type = std::pair<int, int>;

	std::vector<Vertex> vertices;
	std::vector<Face> faces;
	std::map<halfedge_type, int> opps;

	void input_mesh(std::vector<Vertex> const & v, std::vector<Face> const & f) {
		vertices = v;
		faces = f;
	}

	// face_idx 番目の面の halfedge_idx 番目の半辺を取得する
	halfedge_type get_halfedge(int32 const face_idx, int32 const halfedge_idx) {
		halfedge_type result;
		switch (halfedge_idx) {
		case 0:
			// v0 -> v1 の半辺
			result.first = faces[face_idx].v0;
			result.second = faces[face_idx].v1;
			break;
		case 1:
			// v1 -> v2 の半辺
			result.first = faces[face_idx].v1;
			result.second = faces[face_idx].v2;
			break;
		case 2:
			// v2 -> v0 の半辺
			result.first = faces[face_idx].v2;
			result.second = faces[face_idx].v0;
			break;
		}
		return result;
	}

	// 半辺の始点と終点をひっくり返す
	halfedge_type swapped_halfedge(halfedge_type const& halfedge) {
		halfedge_type result = halfedge;
		std::swap(result.first, result.second);
		return result;
	}

	// 半辺データ構造構築
	void halfedge_construction() {
		// 半辺により境界辺を探す
		opps.clear();

		for (int i = 0; i < (int)faces.size(); ++i) {
			for (int j = 0; j < 3; ++j) {
				halfedge_type edge = get_halfedge(i, j);
				if (opps.contains(swapped_halfedge(edge))) {
					// 既に反対辺が含まれているならばお互いにその辺番号を登録
					opps[edge] = -(opps[swapped_halfedge(edge)] + 1);
					opps[swapped_halfedge(edge)] = i * 3 + j;
				}
				else {
					// まだその辺の反対辺が無い場合、自分の辺番号を負にして登録しておく
					opps[edge] = -(i * 3 + j + 1);
				}
			}
		}
	}

	// 半辺が境界か判定
	bool is_boundary(halfedge_type const& halfedge) {
		return opps[halfedge] < 0;
	}

	// 半辺の番号を計算
	int h_idx(halfedge_type const& halfedge) {
		halfedge_type h = swapped_halfedge(halfedge);
		if (opps.contains(h)) {
			return opps[h];
		}
		else {
			return -(opps[halfedge] + 1);
		}
	}

	// 面から半辺生成
	halfedge_type f_to_h(int const face_idx, int const halfedge_idx) {
		return get_halfedge(face_idx, halfedge_idx);
	}

	// 半辺の始点
	int h_src(halfedge_type const& halfedge) {
		return halfedge.first;
	}
    // 半辺の終点
	int h_dst(halfedge_type const& halfedge) {
		return halfedge.second;
	}

	// 半辺番号から面番号取得
	int h_to_f(int const halfedge_idx) {
		return halfedge_idx / 3;
	}
    // 半辺から面番号取得
	int h_to_f(halfedge_type const& halfedge) {
		return h_idx(halfedge) / 3;
	}

	// 半辺の反対辺の半辺番号取得
	int h_opps(halfedge_type const & halfedge) {
		return opps[halfedge];
	}
    // 半辺番号から半辺生成
	halfedge_type h_make(int const halfedge_idx) {
		return get_halfedge(halfedge_idx / 3, halfedge_idx % 3);
	}
    // 半辺の次の半辺取得
	halfedge_type h_next(halfedge_type const & halfedge) {
		int const face_idx = h_to_f(h_idx(halfedge));
		if (h_dst(halfedge) == faces[face_idx].v0) {
			return f_to_h(face_idx, 0);
		}
		else if (h_dst(halfedge) == faces[face_idx].v1) {
			return f_to_h(face_idx, 1);
		}
		else {
			return f_to_h(face_idx, 2);
		}
	}
    // 半辺の前の半辺取得
	halfedge_type h_prev(halfedge_type const& halfedge) {
		return h_next(h_next(halfedge));
	}

	// ある面の辺の対角頂点番号を取得
	int h_to_v_opps(halfedge_type const halfedge) {
		return h_dst(h_next(halfedge));
	};

    // 辺細分
	void edge_split(halfedge_type const& halfedge, Vertex const& pos) {
		if (is_boundary(halfedge)) {
			// 境界辺の辺細分
			int fa{}, fnew{};
			int h0{}, h1{}, h2{};
			int v0{}, v1{}, v2{}, vnew{};

			// 面番号を取得
			fa = h_to_f(halfedge);
			fnew = (int)faces.size();

			// 半辺の番号を取得
			h0 = h_opps(h_prev(halfedge));
			h1 = h_opps(h_next(halfedge));
			h2 = h_opps(halfedge);

			// 頂点番号を取得
			v0 = h_dst(halfedge);
			v1 = h_to_v_opps(halfedge);
			v2 = h_src(halfedge);
			vnew = (int)vertices.size();

			// 頂点を追加
			vertices.emplace_back(pos);

			// 面の更新
			// fa
			faces[fa].v0 = v0;
			faces[fa].v1 = v1;
			faces[fa].v2 = vnew;

			// fnew
			faces.emplace_back(v1, v2, vnew);

			// 半辺の反対辺を構築
			// fa
			opps[{v0, v1}] = h1 < 0 ? -(fa * 3 + 0 + 1) : h1;
			opps[{v1, vnew}] = fnew * 3 + 2;
			opps[{vnew, v0}] = h2 < 0 ? -(fa * 3 + 2 + 1) : h2;

			// fnew
			opps[{v1, v2}] = h0 < 0 ? -(fnew * 3 + 0 + 1) : h0;
			opps[{v2, vnew}] = -(fnew * 3 + 1 + 1);
			opps[{vnew, v1}] = fa * 3 + 1;

			// 外周
			if (opps.contains({ v1, v0 })) {
				opps[{v1, v0}] = fa * 3 + 0;
			}
			if (opps.contains({ v2, v1 })) {
				opps[{v2, v1}] = fnew * 3 + 0;
			}
		}
		else {
			// 内部辺の辺細分
			int fa{}, fb{}, fnew0{}, fnew1{};
			int h0{}, h1, h2, h3{};
			int v0{}, v1{}, v2{}, v3{}, vnew{};

			// 面番号を取得
			fa = h_to_f(halfedge);
			fb = h_to_f(h_opps(halfedge));
			fnew0 = (int)faces.size();
			fnew1 = (int)faces.size() + 1;

			// 半辺の番号を取得
			h0 = h_opps(h_prev(halfedge));
			h1 = h_opps(h_next(halfedge));
			h2 = h_opps(h_prev(swapped_halfedge(halfedge)));
			h3 = h_opps(h_next(swapped_halfedge(halfedge)));

			// 頂点番号を取得
			v0 = h_dst(halfedge);
			v1 = h_to_v_opps(halfedge);
			v2 = h_to_v_opps(swapped_halfedge(halfedge));
			v3 = h_src(halfedge);
			vnew = (int)vertices.size();

			// 頂点を追加
			vertices.emplace_back(pos);

			// 面の更新
			// fa
			faces[fa].v0 = v0;
			faces[fa].v1 = v1;
			faces[fa].v2 = vnew;

			// fb
			faces[fb].v0 = v2;
			faces[fb].v1 = v0;
			faces[fb].v2 = vnew;

			// fnew0
			faces.emplace_back(v1, v3, vnew);

			// fnew1
			faces.emplace_back(v3, v2, vnew);

			// 半辺の反対辺を更新
			// fa
			opps[{v0, v1}] = h1 < 0 ? -(fa * 3 + 0 + 1) : h1;
			opps[{v1, vnew}] = fnew0 * 3 + 2;
			opps[{vnew, v0}] = fb * 3 + 1;

			// fb
			opps[{v2, v0}] = h2 < 0 ? -(fb * 3 + 0 + 1) : h2;
			opps[{v0, vnew}] = fa * 3 + 2;
			opps[{vnew, v2}] = fnew1 * 3 + 1;

			// fnew0
			opps[{v1, v3}] = h0 < 0 ? -(fnew0 * 3 + 0 + 1) : h0;
			opps[{v3, vnew}] = fnew1 * 3 + 2;
			opps[{vnew, v1}] = fa * 3 + 1;

			// fnew1
			opps[{v3, v2}] = h3 < 0 ? -(fnew1 * 3 + 0 + 1) : h3;
			opps[{v2, vnew}] = fb * 3 + 2;
			opps[{vnew, v3}] = fnew0 * 3 + 1;

			// 外周
			if (opps.contains({ v1, v0 })) {
				opps[{v1, v0}] = fa * 3 + 0;
			}

			if (opps.contains({ v0, v2 })) {
				opps[{v0, v2}] = fb * 3 + 0;
			}

			if (opps.contains({ v2, v3 })) {
				opps[{v2, v3}] = fnew1 * 3 + 0;
			}

			if (opps.contains({ v3, v1 })) {
				opps[{v3, v1}] = fnew0 * 3 + 0;
			}
		}
	}
};
```

</div></details>

# まとめ
メッシュの操作を行うための半辺データ構造を理解しつつ実装しました。半辺を通じて、メッシュ上の様々な情報にアクセスすることができます。

ここで紹介した実装は色々改善の余地があると思います。遊んでみてください。