推しを数式化するアレ

#はじめに
この記事はTokyo Submerged UniversityTokyo City University Advent Calendar 2019 24日目の記事です.

Qiitaはどんなクソを垂れ流しても"それっぽく"見えるようになると聞いたので, 今回はこのプラットフォームを利用させてもらいました.
まるでトイレを題材にした現代アートですね.
#推しを数式で表すとは

まつりちゃんを数式で表した夏色まつり近似曲線を作ったぞ＼＼\\٩( 'ω' )و //／／#祭絵 pic.twitter.com/1zXN8HLdzc

— Nulldehyde (@mt_1682) December 27, 2018

こう言うことですね. 因みにこの子は元気いっぱい! チア部女子高生の夏色まつりちゃんです.

#方法
抽出した輪郭線をフーリエ級数展開で数式化します.
今回はTCUバードの翼くんを数式の世界に堕としていきます.

Q. そもそもフーリエ級数展開ってなんぞ?
A. 「スゥ...ジョゼフフーリエです...」

フーリエ級数展開 pic.twitter.com/nz101gwKqQ

— Nulldehyde (@mt_1682) May 14, 2019

これはフーリエ級数展開のイメージで, 感覚としては「時間と大きさ」の世界を「周波数と大きさ」の世界に変換する感じです. 別の言い方をすると周期的な複雑な波を正弦波の組み合わせに分解することです. 詳しくはggってね.

手順

大まかな手順は

1. 画像のエッジ検出

1. 検出したエッジを良い感じに加工
2. 作成した線をフーリエ級数展開

となります.

コードは下の方にまとめてあるので, 記事中では基本的に日本語での説明のみになります.
#####1. 画像のエッジ検出
pythonとOpenCVで輪郭を検出しようとするとこうなります

import cv2

input = "TCU_Bird_0.png"
mid = "gray.png"
output = "edge.png"

graying = cv2.imread(input, cv2.IMREAD_GRAYSCALE)
detedge = cv2.Canny(graying, 50, 110)

cv2.imwrite(output, detedge)

良い感じに抽出されてます.

因みにWolframに魂を売るとMathematicaの守護霊を憑けてもらうことができ, 画像読み込み, エッジ検出がこれだけで済むようになります. 便利ですね.

img = Import["画像ファイルのアドレス"]
EdgeDetect[img]

これでめでたく画像中のエッジを抽出できるようになり, 「推しを数式化する」という我々の野望に一歩近づいたわけです. 主はWolfram信者なので以降はMathematicaで話を進めて行きます. #####2. 検出したエッジを良い感じに加工する

MathematicaのEdgeDetect関数では例によって例の如く輪郭を線ではなく点のリストとして出力しているので, フーリエ級数展開するためにも何とかしてその点郡を線として扱いたいわけです.
そこで適当な点を始点として, 近い点同士を繋いでEdgeDetectで検出した輪郭の点郡をいくつかの線に分割し, 1本の線として組み合わせた点郡をそれぞれ線として扱えるようにグループ化します.
...
...
文字に起こしてもよくわからないので図を用意しました.

なんとなく理解できれば十分です.

これは実際に点郡を線化した上で, 各線を色別に表示した画像です.

輪郭検出しただけの白黒画像と見比べると, ある程度離れている部分は別の線として扱われているのがわかるかと思います.
ただ, このままでは「線として見做せるようにした点郡」に過ぎないので, Bスプライン曲線を作成して関数化してやります.

これでフーリエ級数展開を行う準備は全て整いました.

#####3. 作成した線をフーリエ級数展開
曲線が閉じているか開いているかをちょいちょい気にしながらフーリエ係数を求めます.

Q. なぜ開閉を気にするのか?
A. ギブスの現象を回避するため.

WolframのHPにいい感じのプログラムがあったので借りてきました.
まずこれは閉じた曲線

綺麗に星形を再現できてます.
次に閉じてない曲線

第3, 4象限でｳﾆｮｳﾆｮしてるのがわかります.
このｳﾆｮｳﾆｮがいわゆるギブスの現象と言われる物で, 余計な成分を持ってくるので排除したかったのです.

あとは適当に数式を表示させれば終了です.
これが推しを堕とし込んだ数式ですね.

結果

{{-((1759 cos(t))/6)+13/19 cos(2 t)+<<31>>+2311/210 sin(12 t)+11437/25,<<1>>},
{(606 cos(t))/19+3388/23 cos(2 t)+<<32>>+17701/29,<<1>>},
{-((892 cos(t))/19)-184/25 cos(2 t)+<<32>>+29/21 sin(12 t)+9881/26,-(<<1>>/18)+<<36>>},
{-((418 cos(t))/25)+<<35>>+11/19 sin(12 t)+6856/21,-((17 <<1>>)/13)+<<34>>},
{-((45 cos(t))/16)+973/19 cos(2 t)+617/88 cos(3 t)+<<30>>+4537/17,<<1>>/19+<<42>>},
{(110 cos(t))/21+<<31>>+17/50 sin(12 t)+15526/31,<<1>>},
{-((883 cos(t))/19)+173/16 cos(2 t)-226/23 cos(3 t)+<<30>>+2391/7,(57<<1>><<1>>)/25+<<37>>},
{-((1751 cos(t))/44)+<<34>>+31/22 sin(12 t)+1818/11,<<1>>},
{-((1288 cos(t))/25)-33/32 cos(2 t)+<<30>>+11/13 sin(12 t)+2638/17,<<1>>},
{-((2204 cos(t))/23)+15/14 cos(2 t)-210/29 cos(3 t)+<<29>>+3839/10,<<1>>/15+<<36>>},
{(886 cos(t))/29+145/17 cos(2 t)-66/17 cos(3 t)+<<32>>+415,-((76 <<1>>)/11)+<<39>>},
{(823 cos(t))/22-25/38 cos(2 t)-129/22 cos(3 t)+<<32>>+15608/29,<<1>>/13+<<38>>},
{-((409 cos(t))/27)-71/19 cos(2 t)+<<35>>+7/29 sin(12 t)+7191/14,<<1>>},
{(64 cos(t))/23+434/19 cos(2 t)+<<30>>+3/25 sin(12 t)+6869/21,-(<<1>>/19)+<<37>>},
{(1395 cos(t))/82+400/39 cos(2 t)+<<32>>+4549/18,<<1>>},
{-((168 cos(t))/169)+1/6 cos(2 t)+<<32>>+5/11 sin(12 t)+10300/23,<<1>>/63+<<36>>},
{-((227 cos(t))/18)+9/13 cos(2 t)+<<37>>+7663/20,<<1>>},
{<<42>>+8807/23,<<1>>},
{<<38>>+2351/19,(466 cos(t))/15+59/22 cos(2 t)+<<33>>},
{-((82 cos(t))/11)-51/23 cos(2 t)+<<31>>+8/39 sin(12 t)+6131/16,<<1>>},
{(608 cos(t))/15-17/27 cos(2 t)+13/3 cos(3 t)+<<23>>+5/24 sin(12 t)+5962/13,<<1>>},
{-((13 cos(t))/40)+44/35 cos(2 t)+<<29>>+14523/28,(1028 <<1>>)/27+<<28>>},
{<<41>>+6613/21,<<1>>},
{-((29 cos(t))/10)-9/22 cos(2 t)+<<30>>+1/20 sin(12 t)+7013/15,<<1>>},
{-((18 cos(t))/7)+<<33>>+1/49 sin(12 t)+5077/17,<<1>>},
{-((297 cos(t))/14)+<<31>>+1/33 sin(12 t)+2141/9,<<1>>},
{(109 cos(t))/12+<<37>>+9157/15,(64 cos(t))/17-8/17 cos(2 t)+<<37>>},
{(238 cos(t))/13+cos(2 t)+48/25 cos(3 t)+<<21>>+1/25 sin(12 t)+8419/21,<<1>>},
{(113 cos(t))/16+5/23 cos(2 t)+7/20 cos(3 t)+<<35>>+4256/25,<<1>>},
{<<36>>+1/23 sin(12 t)+20256/25,(248 cos(t))/21-1/17 cos(2 t)+<<29>>},
{-((69 cos(t))/16)+<<32>>+1/180 sin(12 t)+2703/7,-((39<<1>>cos<<1>><<1>>))/20)+<<32>>}}

また, 近似する次数を変化させていくと図形がどのように変化していくのか, を見るためにManipulate関数を利用してこんな感じのアニメーションを作成できます.

こうして見るとどのように近似されていくのかがわかって面白いですよね.

まとめたコードはこっちに上げときます.

晴れて推しを数学の世界に堕とせるようになった我々に最早敵はいません.
法律の範囲内で楽しみましょう.
#終わりに
最近はラズベリーパイのraspbianに無料のmathematicaが搭載されているのでそれを使うのもアリですね.
ラズベリーパイの性能的に相当時間はかかりそうですが, 並列計算でなんとかなるかな?
と思って, ラズパイ6枚でクラスタ組んで適当に並列計算させたのですが, よくわかりませんでした.

完成した pic.twitter.com/AprBp4CBqP

— Nulldehyde (@mt_1682) December 27, 2018

そもそもラズパイ6枚でクラスタを組むお金があったら, 学生版のMathematica買った方が安いですね. 実機はコストがアレだけど仮想環境にraspbianをインストールすればもしかしたら...? とやってみました. > 「_ちくしょう! お前はいつもそうだ. このエラーはまるでお前の人生そのものだ. 何かやるにしても基本のシステムがぶっ壊れてやがる. 誰もお前を愛さない_」

よくよく考えてみれば今年の5月にWolfram Engineが無料公開されてましたね.
「Mathematica」や“Wolfram|Alpha”の核「Wolfram Engine」が開発者向けに無償開放
早速Jupyter notebookで弄ってみました.

おおー
Jupyter notebookでWolfram languageが使えてる pic.twitter.com/Sxtykrw0ep

— Nulldehyde (@mt_1682) December 22, 2019

EdgeDetectもちゃんと使えるね pic.twitter.com/BaM6nfj3CP

— Nulldehyde (@mt_1682) December 22, 2019

たのしい.

####アドカレ情報
明日はTokyo City University Advent Calendar 2019最終日, らぴーと氏です.

######参考文献
wolfram blog