Houdini で複数本ユニークな最短経路を見つけ出す方法

こちらの記事は [Houdini Apprentice Advent Calendar 2023](https://qiita.com/advent-calendar/2023/happrentice) 6日目の記事です。

元々は [Houdini Advent Calendar 2023](https://qiita.com/advent-calendar/2023/houdini) の方で書かせた[記事](https://qiita.com/talkyren/items/9fd49cf550495254b462)の一部として書く予定ですが、内容はあまりにも離れすぎたので、独立とした記事をもう一枠を使わせて頂きます。

# はじめに
　この記事は Point 間の最短経路を探すノード `Find Shortest Path SOP` の基本機能である **Edges to Avoid** を利用し、重複しない経路を作り出す方法の一例です。

# 動作環境
Windows 10 Pro 64bit
Houdini FX 20.0.506 (Python 3.10)

# 1. Find Shortest Path SOP の基本セットアップ
　こちらに書かれたのは私が最短経路をプロジェクトに組み込む時によく使う手法ですが、ごく一般的な使い方なので、既に知っている方は読み飛ばしても構いません。

## 1.1. 全体の「地図」を準備

![Fig 01](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/43ddf819-548e-3ff9-afbe-2644fe163426.png)
<fig_01: grid_remesh_network>

　こちらは `Remesh SOP` 使って入力ジオメトリを再分割し Edge のネットワークを生成する方法です。最後の `Convert Line SOP` は Polygon を削除するために追加したノードですが、Polygon のままでも `Find Shortest Path SOP` は動作しますので、単なる見た目や私個人の習慣です。

　普段の利用上もちろん問題はありませんが、欠点として上げると：
- ジオメトリの表面上でしか生成できない
- 各頂点は近くの頂点としか繋がらない
- 入力ジオメトリのサイズが変えるたびにパラメータの数値が合わなくて、Houdini がフリーズしがち

がありますので、最短経路を作る時のベースとして私はあまり使いません。

　代わりに、`Scatter SOP` で Point を生成してから `nearpoints()` や `pcfind()` などの VEX 関数で Polyline のネットワークを繋ぐことが多いです。基本のノード ネットワークは下記のようになります。

![Fig 02](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/3a78fb6b-fda7-703f-5854-c0453759a1ae.png)
<fig_02: grid_scatter_vex_network>

```cpp:generate_polyline (Points)
// Run Over: Points

// Connect near points
float radius = chf("radius");
int maxpoints = chi("max_points");

int near_points[] = nearpoints(0, v@P, radius, maxpoints);
removevalue(near_points, i@ptnum); // remove ptnum itself from search result

foreach (int near_pt_num; near_points) {
    addprim(0, "polyline", i@ptnum, near_pt_num);
}
```

　`Attribute Wrangle SOP`(generate_polyline) の中身の VEX コードは先ず検索最大距離 `radius` と検索最大 Point 数 `maxpoints` のパラメータを引数として、`nearpoints()` 関数で Point の空間座標 `v@P` から近くの Point 番号の配列を取得します。
　なお、`nearpoints()` はあくまで**空間座標**からの距離で検索するので、座標が一致する Point 自身の Point 番号はもちろん返り値に含まれます。`removevalue()` 関数で配列から Point 自身の Point 番号である `i@ptnum` を削除します。
　`foreach` で結果である近くの Point 番号を順番に取り出して、`addprim()` 関数で Point と Polyline で繋ぎます。

　`radius` と `maxpoints` のパラメータをそれぞれ `10` に設定し、結果を確認すれば分かると思いますが、幾つかの Point の間にバイパスのような Polyline が作成されました。更に `Scatter SOP` の **Relax Iterations** パラメータをデフォルトの `10` から `1` に下げて、頂点の分布に少し疎密を付けます。

![Fig 03](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/f2b9aee1-eb9c-dbe0-8294-07c9927f7103.png)
<fig_03: scatter_relax_from_10_to_1>

## 1.2. 起点と終点に印を付ける
　最短経路を探すためにノードに起点と終点の Point を入力する必要があります。もちろん手動で Point の番号を指定することはできますが、Houdini のプロシージャ特性を活かすために Bounding Box (BBox) の最小座標値と最大座標値を使って、*start_points* と *end_points* の Point Group を設定したいと思います。

　なお、`Find Shortest Path SOP` は複数の起点と終点を入力できますが、結果を確認しやすいための起点と終点の数を1つずつにします。

![Fig 04](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/3c3b6b71-3312-baad-082f-0856e489bbcc.png)
<fig_04: set_start_points_and_end_points_group>

　`Group Expression SOP` と `Attribute Wrangle SOP` をどっちを使っても同じ結果になりますが、ノードの表示内容と VEX コードの分かりやすさの視点から見ると、私は `Attribute Wrangle SOP` の方が好みです。（実行効率的に言っても `Attribute Wrangle SOP` の方が微かに上）

　どの方法を使ってもコアとなるのは指定した空間座標から**一番近い Point を1つ取得する**ための `nearpoint()` 関数です。`nearpoints()` と関数名に **s** が1つ違うだけで引数と返り値が全く違うので、ご注意ください。
　`getbbox_min()` と `getbbox_max()` は関数名通りに Bounding Box の最小座標値と最大座標値を取得するので、特に加筆する必要がありません。

### 1.2.1. Group Expression SOP を使う場合
　**Group Type** を `Point` に変更し、**+** ボタンで項目を1つ追加します。
　それぞれの項目に Point Group 名と VEX コードを入力します：

|Group Name|VEXpression|
|:-:|:-|
|`start_points`|`i@ptnum == nearpoint(0, getbbox_min(0))`|
|`end_points`|`i@ptnum == nearpoint(0, getbbox_max(0))`|

### 1.2.2. Attribute Wrangle SOP を使う場合
　**Run Over** を `Detail (only once)` に変更し、下記の VEX コードを入力します。

```cpp:set_start_end_points_group (Detail)
// Run Over: Detail

vector bbox_min = getbbox_min(0);
vector bbox_max = getbbox_max(0);

int pt_start = nearpoint(0, bbox_min);
int pt_end = nearpoint(0, bbox_max);

setpointgroup(0, "start_points", pt_start, 1);
setpointgroup(0, "end_points", pt_end, 1);
```

## 1.3. セットアップ結果を確認

![Fig 05](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/0d790879-446f-260e-20a0-76a377356421.png)
<fig_05: find_shortest_path_setup_test>

　これで最短経路を見つけ出すためのベースとなる「地図」を準備できました。
　ノード ネットワークの最後に `Find Shortest Path SOP` を追加し、**Start Points** と **End Points** パラメータにそれぞれ *start_points* と *end_points* の Point Group 名を入力すれば（**▼** ボタンでも選択できます）、最短経路の結果が表示するようになりました。

## 1.4. 経路コストを加味
　`Find Shortest Path SOP` はデフォルトに Edge の長さを経路のコストとして計算し、その中でコストが一番低い経路を結果として出力するアルゴリズムを使っています。
　**Path Cost** タブに色んなコストに関する追加オプションを設定できますが、この記事では一番シンプルな **Point Cost Attribute** を使います。他のコストの応用例として、例え Edge の傾斜角度をコストとして加味すれば、地図上の急勾配の斜面が通らないようにすることもできます。

![Fig 06](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/6e20b3f4-2685-b7be-7b6e-d6ef8f4f16a4.png)
<fig_06: add_point_cost_attribute_noise>

　`Find Shortest Path SOP` の **Path Cost** タブの **Point Cost Attribute** を ON にして、アトリビュート名を `point_cost` を入力します。更に上流に `Attribute Noise SOP`(init_noise_point_cost) を追加し、同じアトリビュート名が `point_cost` とする Float 型のノイズを追加する。ある程度の初期値が欲しいので、**Range Values** を **Min/Max** に変更し、 **Min Value** を `0.5` に、**Max Value** を `1.0` に設定します。
　`Attribute Noise SOP` の **Offset** スライダーをドラッグすれば、最短経路の結果が変わることを確認できます。

------

　以上は `Find Shortest Path SOP` の基本的なセットアップになります。

# 2. 重複しない最短経路を作成

## 2.1. 全体のフローを考える
　`Find Shortest Path SOP` が出力された経路は Prim Group を付けることができるので、最短経路を複数回作成するには下記のような流れで作成した経路を分離させるのは一番シンプルだと思います。

```mermaid
graph LR
  Start([Edge Network]) --> FSP[Find Shortest Path]
  FSP --> Result((Paths))
  Start --> Mod[/Modify/]
  Result -.-> Mod
  Result --> Path([Shortest Paths])
  Mod --> End([Edge Network])
  End --Next Loop--> Start
```

　経路が重複しないように、反復処理の使って毎回作成された経路の情報をベースにして、元の Edge のネットワークを加工してから次のループに送ります。Houdini の For-Loop の仕様上、最短経路 Prim をループの途中から分離することが難しいので、ループごとの最後に一旦 `Merge SOP` でデータを合併してから、次のループの始まりに `Split SOP` で分離することにします。
　Houdini のノード ネットワークはこのようにセットアップします。（説明のため、加工の箇所は `Attribute Wrangle SOP` で代用します。）

![Fig 07](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/e2c9c6d6-0e61-7a03-2936-4bf2c935ad17.png)
<fig_07: example_of_loop_flow>

　`Block Begin SOP` と `Block End SOP` を手動で作成すると設定する必要な箇所が多くなるので、Tab メニューで `for-loop` を入力して、表示された **For-Loop with Feedback** を使います。

　`Find Shortest Path SOP` の **Output** タブの **Path Group** オプションを ON にすれば、出力された経路 Prim はデフォルトに *paths* という Prim Group を付けられます。この Prim Group を利用すれば、次のループで前ループの経路 Prim をメインの処理に入る前に分離することも簡単です。

## 2.2. 通らない Edge を指定
　次に通らせたくない経路を指定する方法を考えます。

　`Find Shortest Path SOP` の **Surface Constraints** タブに **Edges to Avoid** という項目があって、こちらに Edge 番号・Prim 番号・Edge Group 名・Prim Group 名のいずれかを入力すれば、ノードが最短経路を探索する時にその Edge（もしくは Prim に含まれる Edge）を通らないようにします。
　エレメント番号を入力として使うのが少し面倒なので、ここで Edge Group を使ってセットアップしたいと思います。

:::note info
　Prim Group を使わない理由として、Polygon と Polyline どちらの状態でも **Edges to Avoid** に指定できますが、処理時の挙動が少し違います。Polygon の場合 Prim を構成する全ての Edge が対象になることに対して、Polyline の場合 Prim 自体が Edge そのものになります。
　こうした入力データの状態によって、処理の出力結果が変わることを避けたいので、一義的な Group タイプを選びます。もちろんしっかり入力状態の検知や予備処理などの対策をすれば良いですが、却って不具合の要因になる可能性が高くなりますので、この記事では説明を省かせて頂きます。
:::

------

![Fig 08](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/1bef320c-7308-7db3-1fe1-0a363572134a.png)
<fig_08: A_add_avoid_edges_group>

　先ずは For-Loop に入る前に `Group Create SOP`(init_edge_group) で Edge Group を作成し *avoid_edges* と名前を付けます。Group の内容は空のままで大丈夫ですので、**Base Group** も OFF にします。こちらの Edge Group を作成しないと、初回のループの時 `Find Shortest Path SOP` は指定する Group が見つからなくてエラーになります。
　`Find Shortest Path SOP` の **Surface Constraints** タブの **Edges to Avoid** 項目も ON にして、*avoid_edges* の Group 名を入力します。

## 2.3. 経路の Edge を avoid_edges に追加

![Fig 09](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/1848d3e7-a2a2-5372-e018-38fdf51f1794.png)
<fig_09: A_transfer_to_edge_network>

　`Find Shortest Path SOP` から出力された *paths* Group は Prim Group なので、`Group Promote SOP` で Edge Group に変換します。同時に同名の違うタイプの Group を避けたいので、*to_avoid* の新しい名前に変更します。
　次に `Group Transfer SOP` で元の Edge のネットワークに転写します。**Distance Threshold** を小さい数値の `0.0001` に設定します。なお、初回のループ以降、前ループの *to_avoid* Group 結果が残っていたので、**Group Name Conflict** のオプションを **Overwrite** に変更します。（*to_avoid* はあくまで中間変数にあたる Group のため、`Group Delete SOP` でループごとに削除してもいいですが、上書きするだけで十分だと思います。）
　最後に `Group Combine SOP` で *to_avoid* の内容を既存の *avoid_edges* に合併（Union）します。

## 2.4. 最短経路の結果を確認

![Fig 10](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/70665de1-20c1-448f-9fa3-dfd99e427b8f.png)
<fig_10: A_path_result_visualization>

　For-Loop の最終結果は複数の経路と元の Edge のネットワークが一緒に出力されるので、`Blast SOP` で *paths* の Prim Group のみ取り出します。更に、同じ完全に重ねた Edge があるかどうかを確認するため、`Point Jitter SOP` と `Resample SOP` で結果を簡単に加工します。

------

　`Block End SOP`(repeat_end) の **Iterations** パラメータはデフォルトで `10` しかないので、試しに `20` に反復処理の回数を上げてみると、**出力された経路は一定数以上にならない**ことが分かりました。

　`Find Shortest Path SOP` に :warning: の警告が表示されたので Node Info で確認してみたら `No path exists between a start point and an end point.`（起点と終点を結ぶ経路が存在しません。）と教えてくれました。

![Fig 11](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/37d82af5-7c86-9bb5-b434-d3f3b851c7b5.png)
<fig_11: A_node_info_of_find_shortest_path_about_no_existing_path>

　For-Loop 内の `Group Combine SOP` に表示を切り替えて確認すれば分かると思いますが、*start_points* として指定された Point から出た全ての Edge が *avoid_edges* に指定されたので、使える経路は存在しないことが明白です。

![Fig 12](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/4f7b740f-118a-173d-f187-8bf6ca632bbf.png)
<fig_12: A_group_combine_result_of_over_10_iterations>

　確かに、この記事のタイトル通りに**完全に**ユニークな経路を作成されましたが、初期の Edge ネットワークによって、できる経路の本数があまりにも少なくなることが頻発することになります。このままではアルゴリズムとして使い勝手が悪すぎるので、ある程度条件を緩めてもっと汎用性の高い処理を目指したいです。
　この問題を回避するために必要のは *start_points* と *end_points* に近い Edge を *avoid_edges* に指定されないようにすれば解決になります。但し、除外エリアなどを設定することによって、却ってノード ネットワークが複雑になるので、逆に考えれば、ループごとに *avoid_edges* に指定するための使った *paths* Prim の両端を適切に切断すれば、切断された部分の Edge は *avoid_edges* に転写されなくなり、この問題が修正されると推定します。

# 3. avoid_edges の範囲を再調整

## 3.1. 転写のための経路を切断

![Fig 13](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/61b31c9c-5fd4-9400-9baf-205f4e387740.png)
<fig_13: B_cut_paths_prim_with_carve_sop>

　Polyline 型の Prim を両端から切断するために `Carve SOP`(cutoff_start_end) を使います。他のジオメトリを入力されることはないですが、念のために **Group** を *paths* を指定します。**Second U** を ON にして、**First U** と **Second U** のパラメータを適切な数値に設定します。
　なお、**Divisions** モードに設定すると Edge の途中から切断され、切断箇所から新しい頂点 Point が作成されますので、ここで **Breakpoints** モードに設定し、切断箇所の Edge を丸ごと削除するようにします。

 ------

![Fig 14](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/e1d36018-2033-9566-7e71-5c2992ca9fa1.png)
<fig_14: B_set_edge_group>

　`Carve SOP` で Prim を切断すると、Prim はノードの仕様上、元の Prim Group を保持できないので、ここで [***§2.3***](#23-経路の-edge-を-avoid_edges-に追加) に使われた `Group Promote SOP` は意味がなくなったので削除します。普通に `Group Create SOP` を使って全ての Edge を *to_avoid* の Edge Group に指定します。

## 3.2. 最短経路の結果を再確認

![Fig 15](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/33a14e5b-0aa0-1a0f-8936-607e9f81046d.png)
<fig_15: B_path_result_visualization>

　修正後の結果を再度確認すると、経路がしっかり **Iterations** に設定された通りに 20 本作成されたことが分かりました。しかも、*avoid_edges* 範囲の再調整によって除外された部分は複数本の経路が通ったことを確認できます。
　但し、更に `Block End SOP`(repeat_end) の **Iterations** パラメータの数値を上げて、試しに `100` に設定したら、結局作成された経路の数がある程度限られています。では、`Carve SOP`(cutoff_start_end) の除外範囲をもっと広げれば使える経路は増えるではないかと考えられますが、結局この下の画像のように、途中にある Edge を全部通った後に道が無くなります。

![Fig 16](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/5a165ff2-41e8-f15b-8004-30d7d321a78e.png)
<fig_16: B_upper_cut_range_to_0.3_and_0.7>

　この構造自体が毛細血管網のようになって、割と面白い結果になりましたが、やはり「道が足りない」問題を根本的に解決するには他の方法を考えなければなりません。

# 4. Edge を除外する方法を再考
　Edge を完成に除外すると使える経路が少なくなるので、経路を**通らせない**ではなく、経路が**通りにくくなる**に方針転換すれば、経路の数がある程度十分に確保できるようになります。
　ここで [***§1.4***](#14-経路コストを加味) で全ての頂点に追加された `point_cost` に注目します。一度に通った経路の `point_cost` を少し上げれば、今の道を塞ぐことなく、他の道を選ばせることに「誘導」するになります。**完全に**ユニークな経路でなくなりますが、**ある程度**重複しない経路を作ることができます。

## 4.1. 転写のための経路の Group を変更

![Fig 17](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/3ac20d0c-c9b7-88b3-fbc4-23f542c28385.png)
<fig_17: C_change_path_group_type_and_transfer_type>

　[***§3.1***](#31-転写のための経路を切断) で `Group Create SOP` を使って裁断された経路に *to_avoid* の Edge Group を付けましたが、今度は Edge ではなく各頂点をマークしたので、**Group Type** を **Points** に変更し、名前も用途に合わせて *to_adjust* に変わります。
　`Group Transfer SOP` のパラメータも合わせて Point Group 用に変更します。他のパラメータはそのままで大丈夫です。これで裁断された経路にあった Point が元の Edge のネットワークで *to_adjust* としてマークされました。

## 4.2. point_cost を調整

![Fig 18](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/0e98cd34-333d-eddf-cfe6-3fa5ecaf6729.png)
<fig_18: C_adjust_point_cost>

　[***§2.3***](#23-経路の-edge-を-avoid_edges-に追加) で追加された `Group Combine SOP` を使わなくなったので、`Group Combine SOP` ノードを削除して、`Attribute Adjust Float SOP`(adjust_point_cost) を新たに追加します。
　*to_adjust* の Point Group に対して、`point_cost` アトリビュートを調整します。では、数値をどのように調整すれば良いのでしょうか。一定の数値を加えるだけで効果がちょっと弱いので、ここで **Operation** を **Multiply** （乗算）に変更し、乗数の **Constant Value** を `5` と入力します。最後に **Enable Post-Process** を ON にして、**Maximum** の数値を `3.4e+35` （$3.4 \times 10 ^ {35}$）に設定します。

:::note info
　普段 Houdini を使う時にあまり気にしなくていい問題ですが、一応 `float` 型の数値の最大値が $3.4 \times 10 ^ {38}$ になるので、オーバーフローのことを配慮し、`point_cost` の最大値をその一歩手前に制限します。因みに $\text{log}_5 (3.4 \times 10 ^ {35}) \approx 50.834$ なので、割と現実的な数値です。
:::

## 4.3. 複数経路の結果を確認

![Fig 19](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/7b7eec70-c328-9073-736b-d8c7bec40a11.png)
<fig_19: C_iterations_with_300>

　`Block End SOP`(repeat_end) の **Iterations** パラメータの数値を `300` に設定し、少し処理の時間がかかりますが、きちんと 300 本の経路を作成されたことを確認できます。

------

![Fig 20](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/420722/c784aa54-c919-aa16-a5f7-612af484ac4f.png)
<fig_20: C_with_different_input_geo>

　試しに最初の入力に Gird ではなく Volume(VDB) や Pig Head を使っても、経路も期待通りに作成されました。

　ここから更に For-Loop に色んな種類のコスト要素を追加することができますが、これ以上にダラダラ書くと記事が長くなりすぎて年内投稿が間に合わなくなるので、皆さんにお任せします。

# おわりに

　いかがでしょうか。`Find Shortest Path SOP` の[公式ドキュメンテーション](https://www.sidefx.com/ja/docs/houdini/nodes/sop/findshortestpath.html)を読めば割とすぐに作れそうな内容ですが、まさかこんな長い記事になりました。私自身もこの記事を書くため、改めて色んなノードのパラメータを細部まで調べて、良い勉強の機会にもなりました。
　Group を操作するノードではなく、`Attribute Wrangle SOP` で VEX コードを並ぶバージョンも試作しましたので、時間があれば VEX バージョンも投稿する予定です。なお、私が [Houdini Advent Calendar 2023](https://qiita.com/advent-calendar/2023/houdini) の方で書いた[空間中の Point を補間する方法](https://qiita.com/talkyren/items/9fd49cf550495254b462)と結合し、もう少し複雑なものも作れますので、流石に年内には無理なので、来年の自分に託します。
　それでは、皆さん良いお年を。

　Talky Ren　2023/12/30