関数アートジェネレーター　技術解説

# はじめに
「**関数アート**」と呼ばれるアートをご存知でしょうか。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/3d8e3c93-961b-7bc2-99ac-05d29338e3b4.png" width="40%">
こんなふうに**猫を数式で表現**してみたり$^{[1]}$、時には

https://www.youtube.com/watch?v=g3gowoSlOvI

こんなふうに**音楽のPV**そのものを数式で表してしまったりする、数学好きならではの**職人芸**です。

ところでこの関数アート、**自動生成**してしまえそうではないでしょうか。
ブラウザ上で動作する[「**関数アートジェネレーター**」](https://yodai49.github.io/functionArtGenerator/)を制作したので、このツールが自動的に関数アートを生成している**仕組みについて解説**します。[ **技育展2022**](https://talent.supporterz.jp/geekten/2022/)に出展しました。開発期間は**正味1週間ほど**です。

# 作ったもの

https://yodai49.github.io/functionArtGenerator/

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/461caead-8e79-6772-2394-067b92423252.png" width="70%">

お手持ちの画像から、**関数アート**を**自動的に生成**することができます。
検出精度は現状**まずまずの水準**$^{[2]}$で、単純なイラストなら**いい感じに**関数アートを生成してくれます。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/258eaa6e-75b4-5dea-8a72-742c194a6d9c.png" width="35%"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/29256552-948d-e95e-1b4a-01fb22a49ee8.jpeg" width="35%">
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/64f163d5-96f7-6900-3820-b98f8b977514.png" width="35%"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/f211e542-64d3-3332-28dd-83a0aeb4b7f2.jpeg" width="35%">
**写真**でもそこそこの検出精度です。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/c90a4938-a0a2-4aee-d32f-398b2bdc35d2.jpeg" width="35%"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/3908e62c-311b-ded3-10c0-a2b8b3d5525f.jpeg" width="35%">

**ブラウザ上で動作**するため、インストールやダウンロード不要で**誰でもすぐに**利用することができます。
なお、現時点で生成される関数の種類は「**線分**」「**円弧**」「**多項式**」の3種類です。
画像認識のライブラリや、機械学習のライブラリなどは用いず、**ゼロからフルスクラッチで開発**しています。
**フーリエ変換**などを用いれば非自明な形の関数アートが生成できそうですが、できる限り**人間が理解しやすい関数アートを生成**するため、今回はフーリエ変換を用いない仕様となっています。


> 現時点では、座標平面はx,yがそれぞれ-1〜1の範囲で固定です。
また、生成された関数アートの式を修正することはできません。
これらの機能は、そのうち**追加実装**するかもしれません。

# 概観
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/3e590a92-9513-77c3-d8b4-3eedc60bcc6e.png" width="90%">

今回製作した**関数アートジェネレーター**は、以下の**6つのステップ**で関数アートを生成します。
**1**. **エッジの検出**　画像の**エッジ**を検出します。
**2**. **勾配の計算**　　画像中の各点について、**局所的な傾き**を算出します。
**3**. **線分の検出**　　画像に含まれる**線分**を検出します。
**4**. **円弧の検出**　　画像に含まれる**円弧**を検出します。
**5**. **グラフ構造化**　ここまでで近似できなかった点を、**グラフ構造化**します。
**6**. **多項式で近似**　グラフ構造化した各連結成分を、**多項式で近似**します。

ここから、各ステップの動作について解説します。
なお、実装には**JavaScript**と、**Canvas要素**を用いています。

## ステップ１　エッジ検出
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/64f163d5-96f7-6900-3820-b98f8b977514.png" width="35%"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/c68249df-67c1-b37c-9ce4-af502ede28a5.png" width="35%>

　
ステップ1では、**エッジを検出**して画像を関数表現しやすくします。
**突貫工事の開発**だったこともあり、このステップはかなり**プリミティブなエッジ検出**を実装しています。


**1.** 画像を1ピクセルずつ**走査**する
**2.** 右隣のピクセルが存在するとき、そのピクセルとの「**差**」が閾値以上であれば、該当する点を中心に**一定の半径の円**を描く
**3.** 終端に到達したら終了

ここで、$2$ピクセル$\mathrm{P_1}(r_1,g_1,b_1),\ \mathrm{P_2}(r_2,g_2,b_2)$間の**差**とは、
$$d=\sqrt{(r_1-r_2)^2+(g_1-g_2)^2+(b_1-b_2)^2}$$
で定義しています。実際のプログラムでは、**平方根の計算は省略**しています。

この部分を、**輝度**を与える式を用いて
$$d=|0.3(r_1-r_2)+0.6(g_1-g_2)+0.1(b_1-b_2)|$$
などと改良すれば、より望ましい結果が得られるかもしれません。また、**垂直方向の差分**も考慮したり、[**キャニー法**](https://jp.mathworks.com/discovery/edge-detection.html)を用いるように改良したりすれば、よりよい結果が得られます。

とはいえ、関数アートジェネレーターの上では**十分望ましい結果**が得られているため、今回はこれでよしとします。

## ステップ２　勾配の計算
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/c68249df-67c1-b37c-9ce4-af502ede28a5.png" width="35%"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/34a1bb7d-d8b6-7c70-fd34-6825c063962f.png" width="35%">

エッジを検出できたら、次は検出されたエッジ上のピクセルそれぞれについて、付近の**局所的な傾きを算出**します。傾きを算出することで、次の処理である「**線分の検出**」と「**円弧の検出**」の**正確性**を高めることができます。

なおこの処理は**やや負荷が高い**ため、適当にエッジ上の**ピクセルを間引いて**処理を行なっています。

以下の図において、**局所的な傾き**を計算したい点が**赤い点**、**サンプリングする点**が**青枠の点**と**青い点**、**灰色の部分**が**検出されたエッジ**を表します。また、赤い点を$\mathrm{P}$、サンプリングする直線の角度を$\alpha$、サンプリングする各点を**図**に示した要領で$\mathrm{P_{-n},P_{-(n-1)},\cdots,P_{n-1},P_n}$とします。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/5374f271-6c0e-eaa9-e478-ebc1e347f2fc.png" width="35%">

このとき、**局所的な傾き**は以下のステップで算出します。

**1.** エッジ上の各点を順番に走査する
**2.** $\alpha:0\leq\alpha<\pi$を一つとる。
**3.** $\mathrm{P}$を通る傾き$\tan\alpha$の直線上の点を一定間隔で**サンプリング**する。
**4.** サンプリングした各点のうち、**エッジ上にあった点の割合**$r_\alpha$を求める。
**5.** $r$の**最大値**を与える$\alpha$を$\alpha_{max}$とする。
**6.** $r_{\alpha_{max}}$と$r$の平均$\overline{r}$の差$d_{\mathrm{P}}$を求め、これが閾値以上だったら$\alpha_{max}$を点$\mathrm{P}$の**傾きとして採用**する。

この処理を行った結果得られた局所的な傾きを**図示**したのが以下の図です。
※処理を分かりやすくするため、検出された勾配の**一部のみを描画**しています。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/09506a50-3842-4ade-4a96-14ebdb3bf11d.png" width="35%">

ここからさらに、**二次元配列を用意**して、得られた各点の傾きを**二次元配列上のデータに移し替え**ます$^{[3]}$。
この二次元配列を**可視化**したものが、本項目の**冒頭**で示した**右側の図**であり、これが**ステップ２のゴール**です。

**処理時間**が気になりますが、例で示している**キノコの画像**の場合、主要な**データ**は以下の通りです。
- **ピクセル数**：　$540\times540$
- 傾きを調査した**ピクセルの個数**：　$1165$　(縦横4ピクセル間隔で調査）
- サンプリングの**角度に対する分割数**： $90$（2度ずつ直線をずらす）
- サンプリングの**半径に対する分割数**： $30$（片方の方向あたり）
- $s_\mathrm{P}$の**閾値**： $0.25$
- 閾値を超えた傾きの**検知数**： $1070$

これより、**計算した回数の概算**は以下の通りです。
- 各点あたりの**傾きの計算**：$90\times(30\times2)\simeq5400$
- 全ての点の**傾きの計算**：$5400\times1165\simeq6.3\times10^6$

極端な話、仮に**全ての点がエッジ**だったとしても縦横$4$ピクセルずつ間引いて傾きを算出していけば、**概算計算量は$10^8$を超えない**ので、任意の画像について**現実的な時間で処理**できます。

## ステップ３　線分の検出
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/34a1bb7d-d8b6-7c70-fd34-6825c063962f.png" width="35%">
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/2619265f-cb61-64a0-2f29-256d6c8fb404.png" width="35%">

まずは**直線を検出**し、あとからその直線の上にある**線分を抽出**するという流れです。

はじめに、**拡張されたHough変換**を用いて直線を検出します。
拡張されたHough変換の解説については、**別途記事を書きました**。
Hough変換について知りたい方は、**以下の記事**をご覧ください。

https://qiita.com/trth/items/dbc18a7845e81d050886

その後、**直線の補正**を行います。直線の補正は以下のステップで行います。
$N_l$本の直線が検出されたとき、**検出された直線**を$l_n\ (0\leq n< N_l)$とし、その式が$$l_n:\ (\sin\alpha)x-(\cos\alpha)y+\beta=0$$で与えられるとします。

**1.** 検出された直線$l_k$を一つ選ぶ
**2.** $\alpha+i\delta_\alpha,\ \beta_k+j\delta_\beta$をパラメータに持つ直線$l_{ij}$について、$l_{ij}$上にある**エッジ部分の割合$r_{ij}$を調査**する。なお、$i,\ j$は適当な範囲内の整数であり、$\delta_\alpha,\ \delta_\beta$は正の適当な実数である。
**3.** $r_{ij}$が最も大きかった$i,\ j$の組を$i_{max},\ j_{max}$とし、$l_k$のパラメータをそれぞれ$\alpha_{i_{max}},\ \beta_{j_{max}}$で置き換える。
**4.** $\delta_\alpha,\ \delta_\beta$を既に調査した範囲と被りが出ないよう**定数倍**して、同様の処理を繰り返す。

ここまでが、直線検出のステップです。

その後、検出された**各直線$l_k$上の点をサンプリング**し、**エッジ上にあるかどうかを調査**します。
$2$つ以上連続して**エッジ上にある点の組**を検出し、該当する範囲を「**線分あり**」と判定します。

## ステップ４　円弧の検出
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/34a1bb7d-d8b6-7c70-fd34-6825c063962f.png" width="35%">
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/f0055333-456d-ce68-3121-e960e7210c53.png" width="35%">

線分の時と同様に、まずは**円を検出**し、あとからその円の上にある**円弧を抽出**します。

こちらも、[**拡張されたHough変換**](https://qiita.com/trth/items/dbc18a7845e81d050886)を用いてまずは**円の中心を抽出**します。
ステップ2の段階で各点の**局所的な特徴を一次式で取り出している**ため、ここでは**円の半径の情報はわからない**ことに注意が必要です。

その後、**円の半径を決定**します。
**円の半径**は、次のような流れで決定します。
$N_C$個の円の中心が検出されたとき、その中心を$\mathrm{P_n}\ (0\leq n<N_C)$と表します。

以下の図において、**灰色**は**画像のエッジ**、**赤い点**は$\mathrm{P}_n$を表します。
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/1db79617-3494-201a-b2f7-aa322a36f3b3.png" width="35%">

1. 検出された円の中心$\mathrm{P}_k$を一つ選ぶ。
2. $r=i\delta$を一定間隔で一つとる。ここで$i$は$0<i\leq i_R$を満たす自然数であり、$i_R$は適当な自然数、$\delta$は適当な正の実数である。
3. $\mathrm{P}_k$を中心とする半径$r$の円上から、**一定間隔で点をサンプリング**する。
4. サンプリングした各点のうち、**エッジ上にあった点の割合**$s_i$を求める。
5. $t_1=s_2-s_1,\ t_2=s_3-s_2,\ \cdots, t_{i_R-1}=s_{i_R}-s_{i_R-1}$を求める。
6. $t_1,\ t_2,\ \cdots t_{i_R-1}$のうち、$t_{n-1}t_{n}<0$となる$n$を見つける（**※**）。
7. 半径$n\delta$の円を検出する。

**※** $t_1,\ t_2,\ \cdots t_{i_R-1}$と順にみていったとき、**正負が逆転する点**を探しています。
$t$が連続関数の場合の**微分に相当**すると考えると分かりやすいです。

ここまでが、円検出のステップです。

その後、線分検出の場合と同様、**検出された各円上の点をサンプリング**し、**エッジ上にあるかどうかを調査**します。
$2$つ以上連続して**エッジ上にある点の組**を検出して「**円弧あり**」と判定します。
この際、$x$軸と交差するような円弧を別物と認識してしまわないように注意が必要です。
今回の実装では、**エッジ上にない円弧上の点**をまず見つけ、**そこから一周**して調査しています。

## ステップ５　グラフ構造化
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/8367b9f9-cda6-de83-73c7-1c5877644fe5.png" width="35%"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/6978fc2d-98ec-b2e8-0284-d50f4e3b152a.png" width="35%">

このステップでは、**多項式で近似するための準備**を行います。

まず、**ステップ3**と**ステップ4**で検出された部分に相当するベクトル場の大きさを**強制的に$0$とします**。

その後、**ベクトル場の値が0となっていない点をエッジ通りに結ぶグラフ**を以下の要領で作成します。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/1ada557b-a31a-b6b6-7f97-084270a649e8.png" width="35%">

**1.** ステップ2で傾きが検出されていた点で、なおかつまだその付近のベクトル場の大きさが$0$ではないものを一つ選ぶ。
**2.** その場の傾き方向の周辺で、半径$r_0+i\delta$の円周上の点に、**エッジ上の点**がないかどうかを調べる。ここで、$r_0,\ \delta$は適当な正の実数、$i$は$0$から適当な整数の閾値未満の整数である。
**3.** エッジ上の点が見つかったらそこを**新たな頂点**とし、$2$点間に連結を生成したあとで同様の探索を行う。
**4.** エッジ上の点が見つからなかったら、$i$を**インクリメント**して同様の探索を行う。
**5.** $i$が閾値を超えたら**探索を終了**し、その頂点は**グラフの端点**とする。

上部の図の**ピンク色の円弧**が**2**.で探索をおこなったラインです。ピンク色の円弧の間隔が狭いところはグラフの**連結成分**が生成されています。

この探索を行うと、線分や円弧で取りきれなかった**エッジ上の点を滑らか$^{[4]}$に結ぶ連続成分**をいくつかもった**グラフ構造**が得られます。
ステップ5の**冒頭に示した右図**をよく見ると、線分と円弧で取り去れなかったベクトル場の付近に、**連結成分ごとに色分けされたグラフ**があるのがわかるかと思います。

## ステップ６　多項式で近似
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/6978fc2d-98ec-b2e8-0284-d50f4e3b152a.png" width="35%"><img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2519394/2f22476e-5583-ae93-3406-590ae7cfd3a5.png" width="35%">

最後のステップで、連続成分をそれぞれ**多項式で近似**します。

$x$と$y$の間の式を得る手順は、以下の通りです。

**1.** 連結成分を一つ選ぶ。

**2.** [**確率的勾配降下法**](https://ja.wikipedia.org/wiki/%E7%A2%BA%E7%8E%87%E7%9A%84%E5%8B%BE%E9%85%8D%E9%99%8D%E4%B8%8B%E6%B3%95)を用いて、連結成分に含まれる**各点の$y$座標**に対して、**各点の$x$座標の累乗**による重回帰を行う。

**3.** **2**.と同様に、**各点の$x$座標**に対して、**各点の$y$座標の累乗**による重回帰を行う。

**4.** **2**.と**3**.で得られた近似式のうち、**残差の二乗和がより小さい方**をその**連結成分の近似式**として採用する。

**1**.と**2**.においては、今回は正則化は行わずに**残差の二乗和の最小化**を行なっています。重回帰の次数は、「**事前に設定した閾値**」と、「**連結成分に含まれる点の個数**」のうち小さい方です。

なお、次数として後者が採用される場合はわざわざ確率的勾配降下法を用いずとも、[LU分解](https://ja.wikipedia.org/wiki/LU%E5%88%86%E8%A7%A3)を用いて**連立一次方程式**を解けばよいです。
しかし、**極端なパラメータ**を持つ式が生成されることは**関数アート的には望ましくない**ので、今回は確率的勾配降下法に統一しています。

**1**.と**2**.で得られた回帰式における残差の二乗和が**いずれも閾値を上回る場合**は、連結成分を**強制的に中央付近で分割**して再び重回帰を行います。
どちらも望ましくない場合に$x$と$y$それぞれを連結成分内の通し番号$t$で重回帰するようにすれば、パラメータ表示による近似曲線が得られそうです。

# おわりに
関数アートを生成するにあたり、最終的に採用しなかった**いくつかのアプローチ**も考えていましたが、なかなかうまくいかないものも多かったです。
また今回は**時間の制約**上、**実装を諦めたアイディア**もかなりあります。
**技育展2022**でのフィードバックも活かしつつ、さらに**パワーアップ**したものをいつかリリースしたい所存です。
最後までお読みくださってありがとうございました。

### 注釈
[1] 表示されている数式は猫のかかと部分（緑色の部分）の数式です。
[2] かなり曖昧ですが、見た目が不自然ではない、ということです。
[3] 傾きが検出された点それぞれについて、その近くに相当する部分の$2$次元配列の要素にその傾きを影響させる処理を実装しました。
[4] 数学的な意味ではなく、日常生活での意味における「滑らか」です。