【Pythonコードあり】たった一枚のシルエット画像から一瞬で三次元形状を復元する、機械学習なしで

これが
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/250024/9b8fab99-1c31-4f55-eb08-ed9ef0272381.png"  width="256">
一瞬でぷくぷくシールみたいに膨らんでこうなる。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/250024/f045f294-0658-4267-af40-ace65766f2a3.gif"  width="512">

ぷくぷくシール
https://www.google.com/search?q=%E3%81%B7%E3%81%8F%E3%81%B7%E3%81%8F%E3%82%B7%E3%83%BC%E3%83%AB

# 背景
## 前置き　Computer Vision（CV）の分野の話
シルエット画像が三次元形状を復元する上で重要なキーであることはよく知られています。
Principal Pointからシルエットのエッジに飛ばしたレイは対象物の三次元形状をProjectiveに拘束するのでたくさんのシルエット画像（とカメラパラメータ）があればwatertightな対象物の3次元形状を幾何的に推定することが可能です（Visual Hull）。
こうしたShape-from-Silhouetteと呼ばれる分野は古くからあり、精度よくキャリブレーションされた複数視点の入力を前提としていました。
Patch-based MVS（PMVS）で有名な日本人の3D Computer Vision研究者である古川先生も昔はこういうのをやってたみたいですね。
Furukawa, Yasutaka, and Jean Ponce. "Carved visual hulls for image-based modeling." European Conference on Computer Vision. Springer, Berlin, Heidelberg, 2006.
https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.519.3426&rep=rep1&type=pdf
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/250024/5f6403b7-97ed-ea95-797a-b39742ec51a6.png" width="512">

一枚（以上）の画像から人体の3次元形状と色を推定する手法として一躍有名になったPiFUもシルエットを入力としていますね。これはDeep Learningをバリバリ使って幾何的には絶対に解けない問題を解いています。陰関数は最高。
Saito, Shunsuke, et al. "Pifu: Pixel-aligned implicit function for high-resolution clothed human digitization." Proceedings of the IEEE/CVF International Conference on Computer Vision. 2019.
https://shunsukesaito.github.io/PIFu/
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/250024/46d35f13-01d2-8aab-eda2-2e14a64f207b.png" width="512">


最近アツいDiffentiable Rendering（DR）界隈でもマスク画像を入力して初期形状（例えば球）を変形させることで3次元形状を推定するデモが多いですね。
DR界隈では単なる技術デモの意味合いが強いようです（CGだけで実画像でテストされていない、古き良きShape-from-Silhouetteの知見がまっっったく応用されていないなど、よく見ると突っ込みどころが多数ありブチ切れそうになります）。
Liu, Shichen, et al. "Soft rasterizer: Differentiable rendering for unsupervised single-view mesh reconstruction." arXiv preprint arXiv:1901.05567 (2019).
https://github.com/ShichenLiu/SoftRas
<img src="https://raw.githubusercontent.com/ShichenLiu/SoftRas/master/data/media/demo/deform/sphere_to_car.gif">

以上は主にComputer Vision (CV)の分野の話でした。
**CVでは現実世界の形状を精度よく復元することが目的**になっています。

## 本命　Computer Graphics（CG）の分野の話
一方でComputer Graphics (CG)の分野では違う文脈でシルエットからの3次元復元の研究が行われていました。
**Human Computer Interaction（HCI）寄りのCGではユーザーが入力した2次元のスケッチから3次元（あるいは2.5次元）の形状をつくるという研究分野があります。
つまり現実には存在しない形状の制作支援をするという技術ですね。**
シルエット画像から高さマップ（height map）やデプスマップ（depth map）を求め、さらにそれに面を張ったメッシュを出力することから**Inflation**と呼ばれているようです。
（シルエットを）膨らませるという意味で、経済のインフレと同じ単語ですね。

**Repoussé**という論文では3Dフォント作成が応用例として挙がっています。
Joshi, Pushkar, and Nathan A. Carr. "Repoussé: Automatic Inflation of 2D Artwork." SBM. 2008.
http://pushkarjoshi.org/pdf/JC_SBIM2008.pdf
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/250024/ef17cdfe-4889-dc52-9242-25f70e11d1be.png)


こういうのが進化してインタラクションが入ると**五十嵐先生のTeddy**になったり
Igarashi, Takeo, Satoshi Matsuoka, and Hidehiko Tanaka. "Teddy: a sketching interface for 3D freeform design." ACM SIGGRAPH 2006 Courses. 2006. 11-es.
https://www-ui.is.s.u-tokyo.ac.jp/~takeo/teddy/teddy-j.htm
(Unityプラグインがあるの今知った)
https://assetstore.unity.com/packages/tools/modeling/teddy-99075
<iframe width="560" height="315" src="https://www.youtube.com/embed/r6X21Py2l2g" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

アニメーション生成が入ると**Monster Mash**などになっていくわけですね。
Dvorožňák, Marek, et al. "Monster mash: a single-view approach to casual 3D modeling and animation." ACM Transactions on Graphics (TOG) 39.6 (2020): 1-12.
https://monstermash.zone/
<iframe width="560" height="315" src="https://www.youtube.com/embed/qQ-8GZ-MvrQ" title="YouTube video player" frameborder="0" allow="accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture" allowfullscreen></iframe>

この記事ではこっちのCGよりの技術のうちでかなり基礎的な手法の解説をします。

# 解説
## 本命　ポアソン方程式に基づく手法
ポアソン方程式（Poisson's equation、Poisson equation）は物理学で流体や電磁気などいろんな分野ででてきます。
https://ja.wikipedia.org/wiki/%E3%83%9D%E3%82%A2%E3%82%BD%E3%83%B3%E6%96%B9%E7%A8%8B%E5%BC%8F

重要なことは**二階微分=ラプラシアン=滑らかさ**　が入ってる方程式だということです。
**2次元のシルエット画像から推定する3次元形状は滑らかに変化してほしいですよね？？？？？？（超天下り的）**
というわけでポアソン方程式を解くのがこの問題, Inflationの常套手段なわけです。

みんな大好きジオメトリ界隈の若き彗星アレックヤコブソン（Alec Jacobson）大先生謹製のMatlab実装があります。
https://github.com/alecjacobson/gptoolbox/blob/master/mesh/inflate.m

ポアソン方程式にさらにどういう設定を入れるか、という点でいくつかバリエーションがあるようです。
もともとのポアソン方程式の原則はシルエットから計算したデプスや高さが滑らかで連続的に変化するというというだけで、絶対値がどうなるべきなのか、という制約はありません。
そこに色々工夫の余地があるわけですね。

残念ながらC++のヘッダオンリーライブラリであるlibiglには実装されていないようです。

この問題は最終的に**疎（スパース）な線型方程式を解くという問題に帰着**します。
いいですか、もう一回言います、**スパースな線型方程式を解くという問題に帰着します**。
つまり、**省メモリかつ高速に解けます**。
楽勝ということです。

### 関係性がなくもない話　ポアソンブレンディング
CV分野の老害がポアソン方程式と聞くとあ、進研ゼミでやったところだ！と思うことでしょう。
そう、コピペした画像を馴染ませる**ポアソンブレンディング**です。
Pérez, Patrick, Michel Gangnet, and Andrew Blake. "Poisson image editing." ACM SIGGRAPH 2003 Papers. 2003. 313-318.
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/250024/14963be0-f812-f1c7-d147-e2bff1a19073.png)


OpenCVにも`seamlessClone`と名前で実装されていますし、
https://docs.opencv.org/4.5.2/df/da0/group__photo__clone.html#ga2bf426e4c93a6b1f21705513dfeca49d

解説記事やMATLAB、Python、C++など実装もいたるところに転がっています。
まあまあ上手くいくし実装も簡単なので一世を風靡した手法なんですね。
実装はこの記事で紹介しているよりInflationより若干複雑ですがかなり似ています。
PRML勉強会が流行る前後にはみんなこれやってたんじゃないでしょうか。


## 誰もが思いつく距離変換に基づくやり方
画像処理をちょっと齧ったことがある人であれば、これは距離変換（Distance Transform）でいいんじゃね？という発想になるはずです。
距離変換とはシルエット画像を入力してシルエット内部のピクセルにおけるエッジからの2次元的な距離を効率的に計算するアルゴリズムの総称です。
エッジからの2次元の距離を高さにすれば目的は達成されるじゃん、という話です。
しかしこれは実は想定通りの結果にならないことが多いです。
なぜかというと**距離変換には勾配の滑らかさ（二次微分、ラプラシアン）の制約がない**からです。
距離変換では各ピクセルを近傍点以上の距離で更新していきますが、雑に言うと右から更新していった距離はどこかで左から更新していった距離とぶつかり、そこで微分不能な尖った点や峰が発生します。

## 比較

冒頭でも紹介した"あ"のシルエットを入力した場合の結果です。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/250024/9cf3eebb-f631-273e-6113-7c8ee6e26098.png)

|ポアソン方程式|距離変換|
|:---:|:---:|
|<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/250024/917ed3ee-2d5b-9930-cc27-1adf77851dd7.png"  width="512">|<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/250024/c5d1c8c3-247d-6f49-0e73-0c45fa87507f.png"  width="512">|

距離変換のほうは尖ってますがポアソン方程式のほうは丸まってますね。
滑らかさ制約のおかげです。

もっとわかりやすい例は円のシルエットを入力したときですかね。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/250024/1c620ffc-e97b-4a27-6881-355307f7815b.png)

|ポアソン方程式|距離変換|
|:---:|:---:|
|<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/250024/32641561-a590-e333-3481-d4f9e72296ce.png"  width="512">|<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/250024/4622fbc9-2fc0-ed6d-968b-61c63c866890.png"  width="512">|

ポアソン方程式は滑らかな半球になっていますが、距離変換は円錐で尖っていることがわかります。


# Pythonで実装
ポアソン方程式に円のシルエットが半球になる制約をいれて実装してみました（参考文献を参照）。

```math
\nabla^2h(x,y)=-4,  \quad \text{subject to } h(\partialΩ) = 0\\
\text{set } z = \sqrt h 
```

これはどういうこっちゃというと円の方程式$x^2+y^2=r^2$を変形した$h(x,y)$を考えて

```math
h(x,y)=-x^2-y^2+r^2=0
```

これにラプラシアン$\nabla^2$を作用させると

```math
\nabla^2h(x,y)=\frac{\partial^2 h}{\partial x^2}+\frac{\partial^2 h}{\partial y^2}=-2-2=-4
```

となるからです。最後にルートをとるのは元々ついている二乗を打ち消すためです。
なお、‐4を4にすると膨らみと凹みが逆になります。

numpyとscipy（sparse matrix用）で実装したコードをGitHubに置いてあります。
https://github.com/unclearness/inflation_py

一部を抜粋します。
愚直に実装しているだけですね。
離散的な画像に対するラプラシアンが4近傍のカーネルを畳み込むことに相当するということを知っていれば読めると思います。

``` python

    ...

    # 4 neighbor laplacian
    for y in range(1, h-1):
        for x in range(1, w-1):
            c = mask[y, x]
            if c == 0:
                continue
            triplets.append([cur_row, img2param_idx[get_idx(x, y)], -4.0])
            kernels = [(y, x - 1), (y, x + 1), (y - 1, x), (y + 1, x)]
            for kernel in kernels:
                jj, ii = kernel
                if mask[jj, ii] != 0:
                    triplets.append([cur_row, img2param_idx[get_idx(ii, jj)], 1.0])
            cur_row += 1  # Go to the next equation
    # Prepare right hand side
    b = np.zeros((num_param, 1))
    rhs = -4.0
    cur_row = 0
    for y in range(1, h-1):
        for x in range(1, w-1):
            c = mask[y, x]
            if c == 0:
                continue
            b[cur_row] = rhs
            cur_row += 1
    if use_sparse:
        # Sparse matrix version
        data, row, col = [], [], []
        for tri in triplets:
            row.append(tri[0])
            col.append(tri[1])
            data.append(tri[2])
        data = np.array(data)
        row = np.array(row, dtype=np.int)
        col = np.array(col, dtype=np.int)
        A = coo_matrix((data, (row, col)), shape=(num_param, num_param))
        x = linalg.spsolve(A, b)

   ...

    for j in range(1, h-1):
        for i in range(1, w-1):
            c = mask[j, i]
            if c == 0:
                continue
            idx = img2param_idx[get_idx(i, j)]
            # setting z = √ h
            depth[j, i] = np.sqrt(x[idx])
```

# まとめ

- いい感じに滑らかに一枚のシルエットを膨らませて3次元（2.5次元）形状にするにはポアソン方程式を解くのが常套手段
- この問題はポアソンブレンディングなどと同様にスパースな線型方程式を解くことに帰着するので省メモリで高速に解ける

# 参考文献
"Notes on Inflating Curves" [Baran and Lehtinen 2009](http://alecjacobson.com/weblog/media/notes-on-inflating-curves-2009-baran.pdf).