Procedural Hard Surface ModelingのためのWorley noiseの実装

この記事は[Houdiniアドベントカレンダー2021](https://qiita.com/advent-calendar/2021/houdini)  25日目の記事です。

# はじめに
[Procedural Design](https://www.procedural-design.com/)の[Akira Saito](https://twitter.com/a_saito)です。
Houdini等を使用しロボットなどのメカデザインを自動で行う研究をしています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/02d2e6bc-2d5c-eafa-0791-a7b0fcf0bea2.png)

今回は入力されたメッシュのシルエットに抑揚をつける為のWorley noiseを使った手法とその実装方法を紹介したいと思います。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/25d1a339-78de-a8e4-e184-75c75735f524.png)
細分化したメッシュにWorley noise(Worley Cellular F2-F1)を加え、左右対称にすると、メカニカルな印象が強まります。
### 動作環境
Houdini 19.0.436(Python3)
macOS Big Sur バージョン11.6
MacBook Pro (15-inch, 2017)
プロセッサ 3.1 GHz クアッドコアIntel Core i7
メモリ 16 GB 2133 MHz LPDDR3
グラフィックス Radeon Pro 560 4 GB





# Worley noiseについて
Worley noiseは、1996年にSteven Worleyによって発表されたノイズ関数です。
空間内のランダムな複数のポイントをから、空間内のすべての場所について、n番目に近いポイント（たとえば、2番目に近いポイント）までの距離を調べ、それらの組み合わせを使用して色などを制御します。


# Houdiniでは様々な方法でWorley noiseを扱うことが可能です。

### Mountain(Attribute Noise) SOP
入力されたメッシュの各頂点を法線方向に様々なノイズで押し出すSOPノードです。
その中の、Worley cellular(F2-F1)を使うことで、工業製品のようなソリッドな印象の形状を生成することが可能です。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/2b930e4d-809f-aae4-3d96-099223a4adf9.png)
Houdini 19へのバージョンアップでMountain SOPは大きく仕様が変わった(Attribute Noise　SOPに置き換わった)ので注意
### Point VOP
![スクリーンショット 2021-12-24 22.54.49.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/2816da25-2951-c094-c611-ce963b4b779f.png)
以下のようなVOPネットワークを組めばMountain SOPと同様のWorley cellular(F2-F1)になります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/28957a28-c2cd-ce76-b275-0d5842fc9b1d.png)

### VEX
Point Wrangle(RunOverがPointのAttribute Wrangle)でも以下のコードを実行すればWorley cellular(F2-F1)が表現可能です。
![スクリーンショット 2021-12-24 23.19.44.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/d73c9c72-56f1-9bf1-3cd2-fc7110978cc8.png)

```c
float f1,f2;
int seed;
wnoise(@P,seed,f1,f2);
@P += (f2 - f1)*1*@N;
i@seed = seed;
```


# 標準のWorley noiseでは、なぜダメなのか？

前述のようにHoudiniでは標準機能でWorley noiseを使うことができますが、自作する理由は２点ほどあります。
Mountain SOPやVEX関数のwnoise()は、ノイズの元となる母点が定型のポワソン分布による一様なものしか使えず、配置そのものや粗密などのコントロールができません。
もう一つの自作する理由は、面情報です。
Worley noiseを使用したメッシュの変形を行なった後に、更なる形状の変更やディティールの追加などを行うために、ポリゴンリダクション等のトポロジーの整理を行う必要があります。
それには、平坦な面とその境界のエッジの情報が必要です。VEXやVOPのWorley noiseならば、それに近いデータを出力することが可能ですが面情報にはなりません。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/7b122e2e-2165-41cf-1e08-783a34cba72f.png)
この様な細分化さらたメッシュでも、面情報が用意できればポリゴンリダクションやトポロジーの最適化も容易になります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/f2cff4e2-21a8-bf30-2b39-c47ef2df6b13.png)
ノイズの折り目で区切られた領域に同じ値を設定しそれを面情報とします。(可視化のために色分けをしていますが、情報としては整数型のアトリビュートです)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/0bad75ed-0f0d-44cf-bba6-e0a88398386d.png)
この面情報を元にポリゴンリダクションを効率よく行うことが可能です。

#Worley noiseの仕組みと実装
## 仕組み
##### まずはシンプルなグリッド状のポリゴンメッシュでWorley noiseの仕組みを解説していきます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/f824ab1a-d373-f5cb-bacf-bfc35200044e.png)
##### 次にノイズの元となる「母点」をランダムに配置します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/426f288f-7db6-a5cc-a072-c29f23df8d9b.png)
##### グリッドの全ての頂点から一番近い母点を求めます。
一番近い母点ごとに色分けするとこの様な見た目になります。
![スクリーンショット 2021-12-25 0.16.02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/9c03ae35-8664-1421-83ac-d7d0b298a98d.png)
##### グリッドの全ての頂点から一番近い母点までの距離の二乗を求めます。
グリッドの各頂点から一番近い母点までの距離を母点に近い方が青、遠い方が赤として可視化するとこの様な見た目になります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/9e7b2699-9fa5-cba0-f535-cac0a336c256.png)
この模様こそが一番シンプルなWorley noise(Worley CellularF1)の見た目になります。
##### この距離を法線方向に押し出しているのがMountain SOPになります
各頂点から一番近い母点までの距離の二乗を各頂点のY座標に加算します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/4a566460-2d95-1814-1b17-244d1678ec77.png)
グリッドの解像度を上げると形がはっきり見えてきます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/739b2708-40e4-b9ca-39f9-75a7b75e4638.png)
さらにグリッドの解像度を上げ母点の数を増やします。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/a54ea0bc-7768-b2d6-ac80-d20392e15a0b.png)
## 実装

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/d3f17dd7-b759-d991-9959-76fffc0a8f4f.png)
この例では母点の配置にscatter SOPを用いていますが、頂点を配置することが出来ればどのような方法でも構いません。
まさに、このことが自作する最初の理由。母点を任意に変更できるという事になります。
配置した母点に母点番号をPoint Wrangleで設定します。

```c
//母点の番号を設定
i@__seed = @ptnum;
```
今回は母点の頂点番号をそのまま母点番号i@__seedとしています

編集対象メッシュ(GRID)を１番、母点を２番目の入力に繋いだpoint wrangleに記述したコードは以下の通り

```c
//@P(現在の処理している頂点の座標)から一番近い頂点番号を入力2から検索する。
int nps[] = nearpoints(1,@P,10000,1);

//入力2の一番近い頂点のi@__seedを求める。
int seed1 = point(1,"__seed",nps[0]);
//入力2の一番近い頂点までの距離の二乗を求める。
float f1 = distance2(@P,point(1,"P",nps[0]));
 
//__seedを現在処理している頂点のi@__seed1に格納する
i@__seed1 = seed1;
//距離を現在処理している頂点のf@__f1に格納する
f@__f1 = f1;
```
i@__seed1に一番近い母点番号

f@__f1に一番近い母点までの距離の二乗が格納されます





# F2-F1 とは？
メカニカルな印象を与えるWorley Cellular　F2-F1とは、どの様な物なのでしょうか？
実装を見ながら説明します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/412c2b9b-4cc9-54bb-07e9-63c7b00fca2d.png)
先ほど、一番近い母点までの距離の二乗を求めたpoint wrangleに数行コードが追加されています。

```c
//@P(現在の処理している頂点の座標)から1番目2番目に近い頂点番号を入力2から検索する。
int nps[] = nearpoints(1,@P,10000,2);

//入力2の一番近い頂点のi@__seedを求める。
int seed1 = point(1,"__seed",nps[0]);
//入力2の一番近い頂点までの距離の二乗を求める。
float f1 = distance2(@P,point(1,"P",nps[0]));

//入力2の二番目に近い頂点のi@__seedを求める。
int seed2 = point(1,"__seed",nps[1]);
//入力2の二番目に近い頂点までの距離の二乗を求める。
float f2 = distance2(@P,point(1,"P",nps[1]));

i@__seed1 = seed1;
i@__seed2 = seed2;

f@__f1 = f1;
f@__f2 = f2;

//seed1とseed2を組み合わせた面情報i@__normalClusterを生成
i@__normalCluster = seed1*npoints(1) + seed2;
```
この様なコードにより、
1番目に近い母点、２番目に近い母点の両方の距離の二乗と、各母点番号を求めています。

### F1
最後のpoint wrangleのコードを以下の様に記載すれば。

```c
@P.y += f@__f1 * chf("heightScale");
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/15eb67fd-4469-7a22-733c-7bb82d352716.png)
一番近い母点までの距離が加味された __Worley Cellular　F1__

### F２
最後のpoint wrangleのコードを以下の様に記載すれば。

```c
@P.y += f@__f2 * chf("heightScale");
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/ba2266c7-c4bd-c34e-af9e-7e92868a07eb.png)
二番目に近い母点までの距離が加味された __Worley Cellular F2__
折り目が増え幾何学的な印象が増しました。

### F2-F1

最後のpoint wrangleのコードを以下の様に記載すれば。

```c
@P.y += (f@__f2 - f@__f1) * chf("heightScale");
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/ec2a06cd-90f1-5600-1c80-0cec112de00b.png)
二番目に近い母点までの距離から一番近い母点までの距離を引く事により、F2の曲面をF1の曲面が打ち消し、直線と平面でできたメカニカルな印象を与えるノイズになります。
これが、__Worley Cellular　F2-F1__
## Worley　Cellular　F2-F1の面情報
形状は求められましたが、各面を定義する面情報はどのように求められるでしょうか？
一番近い母点番号i@__seed1と
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/f33b9ba3-8518-24d1-d201-8ccad5cde670.png)

二番目に近い母点番号i@__seed2

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/e4b0c57c-a4d2-23a2-e442-8f1bf3539c70.png)


この二つを組み合わせることで求められます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/ef7aba33-440e-807f-bb34-1983e58bec32.png)

母点との距離を求めるpoint wrangleの一番最後

```c
//seed1とseed2を組み合わせた面情報i@__normalClusterを生成
i@__normalCluster = seed1*npoints(1) + seed2;
```
この行で、npoints()関数で求めた総母点数と一番近い母点番号を掛け合わせた数に、二番目に近い母点番号を足して求められます。
これにより、自作Worley Noiseの目的の二つ目。面情報i@__normalClusterの取得も可能になりました。



# 面情報を使ったトポロジーの整理
入力メッシュをSphereに置き換え
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/96c5751f-d67f-437c-b7d1-becb58241418.png)

最後のpoint wrangleのコードを以下のように書き換えれば、

```c
@P+= @N * (f@__f2 - f@__f1) * chf("heightScale");
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/25dee7f0-add2-b62a-95f5-3cdc9c8d66a6.png)
Worley　Cellular　F2-F1で法線方向に押し出すことができます。

しかし、この面構成では後の行程が困難です。

そこで、
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/68af8fc8-0fff-d1c0-badc-be7dccf5089d.png)
attribute promote SOPによりpoint attributeである面情報i@__normalClusterを、mode(最頻値)でprim attributeに変換。

変換したprim attribute i@__normalClusterを元にGroup from Attribute Border SOPでedge groupを生成。
最後に、edge smooth SOPでedge groupを指定すれば、より編集しやすいメッシュに修正することが可能です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/2b7566d5-12df-815a-0423-0cce0d7cabc0.png)






# 様々な距離
今まで紹介してきたWorley noiseは各母点までの距離をユークリッド距離という一般的な距離の概念で求めてきましたが、別の概念の距離を使うことでさらに面白い表情を与えることができます。それぞれの距離の概念を紹介していきます。

### ユークリッド距離
![2021-11-09_16h01_15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/4527c7e2-4949-3a4d-c302-d7e1ea7c6aeb.png)

```c
distance = sqrt(A^2 + B^2)
sqrt(4^2+5^2)  =  6.40312423743
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/e518cfe6-480f-c155-cbbd-2130d1b2247a.png)
__ユークリッド距離の面情報__
母点番号の組み合わせが面情報になります

### マンハッタン距離
![2021-11-09_16h01_26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/88541474-4aa6-0c25-8189-8bbd29c9fb1f.png)

各座標の差（の絶対値）の総和を2点間の距離とする。

```c
distance = abs(A) + abs(B)

4 + 5 = 9
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/e28f8505-eb7e-8e8c-cdcf-926f060e2acd.png)
__マンハッタン距離の面情報__
母点番号と母点との位置の差の各要素の符号の組み合わせ８通りにより面の切れ目が発生します。

```c
//一番近い母点の位置
vector p0 = point(1,'P',i1);
//それぞれの位置の差分を求める
vector t1 = p0 - @P;
nc1 = i@__seed1 * 8;
if(t1.x<0) nc1 += 1;
if(t1.y<0) nc1 += 2;
if(t1.z<0) nc1 += 4;
```
### チェビシェフ距離
![2021-11-09_16h01_26.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/ddb30561-fe13-a026-176f-630472febd0a.png)


各座標の差（の絶対値）の最大値の距離とする。

```c
distance = max(abs(A),abs(B))

max(4,5) = 5
```
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/918826/0647d7cb-28c6-6112-244d-fc3c1fa72452.png)

__チェビシェフ距離の面情報__
母点番号と母点との位置の差のでどの要素の絶対値が最大だったかとその要素の符号の組み合わせ６通りにより面の切れ目が発生します。

```c
//一番近い母点の位置
vector p0 = point(1,'P',i1);
//それぞれの位置の差分を求める
vector t1 = p0 - @P;
vector absT1 = abs(t1);
float xyz[] = array(t1.x,t1.y,t1.z);
float absXyz[] = array(absT1.x,absT1.y,absT1.z);
int sortIndex[] = argsort(absXyz);
int maxIndex = sortIndex[2];
float maxDist = xyz[maxIndex];

if(maxDist<0)
{
    maxIndex+=3;
}

nc1 =  i@__seed1 * 6 + maxIndex;
```

### 少し大変なところ
近い母点を求める際に利用したnearpoints()などは基本的に、ユークリッド距離を使用するので、マンハッタン、チェビシェフそれぞれに対応したnearpoints()と同等の関数を用意する必要があります。


# まとめ
プロシージャルモデリングの起点となる事の多い様々なノイズ。
各ノイズの仕組みを理解することで、更なる機能追加や、応用など多くの可能性につながると思います。
気に入ったノイズや、プリセットで用意されている機能の根本まで調べ尽くしてみる事の大切さを、この記事で伝えることができたら幸いです。質問やご指摘などがあれば遠慮なく、[Akira Saito](https://twitter.com/a_saito)までお願いします。

### すみません・・・
公私ともども忙しく締め切りが一番遠い最終日の投稿を選んでしまいましたが、連日の皆様の素晴らしい投稿の最後を飾るには、かなり用途が限定的な話になってしまいました・・・。

__それでは、みなさま。良いクリスマスを！__