イラストを数式に

おはようございます。
今日は少し専門的な数学のお話です。
勉強ではなく、数学の美しさについて少し語っていこうかなと思います。

数学って聞くだけで計算やグラフを思い浮かべて頭が痛くなる方がいるかもしれませんが、
数学にはこんな美しさもあるんだなと思っていただけたら嬉しいです。

Contents

• 1.イラストを数式に
• 2.導入
• 3.準備
• 　1.画像の輪郭を取得
• 　2.輪郭を点群で表す
• 　3.点群にフーリエ変換を施す
• 4.まとめ

イラストを数式に

最初に結果から見ちゃいましょう。

3つの写真があります。

①はネット上で拾った2020東京オリンピックのロゴ
②は①にある処理をかけたもの
③は①にかけるある処理の精度を上げたもの

②から③になるに連れて、①に近づいているように見えます。
ある処理がなんなのかと言うと、フーリエ変換です。

フーリエ変換は大学で数学や物理を専門的に勉強している人しかわからないかもしれませんが、簡単に言うと「関数をSin関数の和で表す」ことです。
この変換元の"関数"は数式で表せない任意関数でも良いのです。

上の例で説明すると
①の輪郭を任意関数とし、フーリエ変換を施したのが②,③となっています。
②と③の違いは「フーリエ変換の精度の違い」ですが、
この精度と言うのはフーリエ変換後の関数の次数のことで、次数が高ければより元のイラストに近いものが生成されます。


つまり「あやふやな状態」から「完璧な状態」に近くことから下のようなとても面白い過程を見ることができます。

今回の記事はこんな感じのgifの作り方。

導入

今回この記事で使用するのは、Wolfram社が提供している数式処理システムMathematica

Mathematica

動作環境 : Wolfram Mathematica10 学生エディション for Mac OS X x86


これはPC上で動作する高性能電卓のようなソフトです。
四則演算はもちろん、中学校で習う連立方程式、高校で習うベクトルの計算、大学で習う微分方程式など、ありとあらゆる計算が可能。
また、計算だけではなく方程式を入力するだけでグラフまで可視化してくれるとても便利なソフトです。

この高性能電卓を用いれば、人間ではとても扱えないような膨大な量のデータを処理することが出来ます。


1.画像の輪郭を取得
2.輪郭を点群で表す
3.点群にフーリエ変換を施す

この順番で処理を見ていきます。

準備

コード簡略化の為に独自関数を5つ定義します。
5つとも先にMathematica上で実行しておくと良いです。

1は最も近い点と点を線で結ぶ関数
2から4はフーリエ変換(Fourier transform)やフーリエ級数(Fourier serie)の関数
5は数式をグラフに表す関数


1.
[mathematica collapse="1"]
definition1[pointList_, neighbothoodSize_: 6] :=
Module[{
L = DeleteDuplicates[pointList],
NF, lambda, lineBag, counter, seenQ, sLB,
nearest, nearest1, nextPoint, couldReverseQ, d, n, s
},
NF = Nearest[L];
lambda = Length[L];
Monitor[
lineBag = {};
counter = 0;
While[
counter < lambda, sLB = {RandomChoice[DeleteCases[L, _?seenQ]]}; seenQ[sLB[[1]]] = True; counter++; couldReverseQ = True; While[( nearest = NF[Last[sLB], {Infinity, neighbothoodSize}]; nearest1 = SortBy[DeleteCases[nearest, _?seenQ], 1. EuclideanDistance[Last[sLB], #] &]; nearest1 =!= {} || couldReverseQ ), If[ nearest1 === {}, sLB = Reverse[sLB]; couldReverseQ = False, nextPoint = If[Length[sLB] <= 3, nearest1[[1]], d = 1. Normalize[ (sLB[[-1]] - sLB[[-2]]) + 1/2 (sLB[[-2]] - sLB[[-3]]) ]; n = {-1, 1} Reverse[d]; s = Sort[{Sqrt[ (d.(# - sLB[[-1]]))^2 + 2 (n.(# - sLB[[-1]]))^2 ], #} & /@ nearest1]; s[[1, 2]]]; AppendTo[sLB, nextPoint]; seenQ[nextPoint] = True; counter++ ] ]; AppendTo[lineBag, sLB] ]; Reverse[SortBy[Select[lineBag, Length[#] > 12 &], Length]],
Grid[{
{
Text[Style["progress point joining", Darker[Green, 0.66]]],
ProgressIndicator[counter/lambda]
}, {
Text[Style["number of segments", Darker[Green, 0.66]]],
Length[lineBag] + 1
}
},
Alignment -> Left, Dividers -> Center
]
]
]
[/mathematica]

2.
[mathematica collapse="1"]
definition2[pointList_, nMax_, op_] :=
Module[{
epsilon = 10^-3, myu = 2^14, M = 10000, s, scale, delta, L,
nds, sMax, if, xyFunction, X, Y, XFT, YFT, type
},
scale = 1. Mean[Table[
Max[fl /@ pointList] - Min[fl /@ pointList],
{fl, {First, Last}}
]];
delta = EuclideanDistance[First[pointList], Last[pointList]];
L = Which[
op === "Closed", type = "Closed";
If[
First[pointList] === Last[pointList], pointList,
Append[pointList, First[pointList]]
],
op === "Open", type = "Open";
pointList, delta == 0., type = "Closed";
pointList, delta/scale < op, type = "Closed"; Append[pointList, First[pointList]], True, type = "Open"; Join[pointList, Rest[Reverse[pointList]]] ]; xyFunction = BSplineFunction[L, SplineDegree -> 4];
nds = NDSolve[
{s'[t] == Sqrt[xyFunction'[t].xyFunction'[t]], s[0] == 0},
s, {t, 0, 1}, MaxSteps -> 10^5, PrecisionGoal -> 4
];
sMax = s[1] /. nds[[1]];
if = Interpolation[
Table[{s[rho] /. nds[[1]], rho}, {rho, 0, 1, 1/M}]
];
X[t_Real] :=
BSplineFunction[L][Max[Min[1, if[(t + Pi)/(2 Pi) sMax]], 0]][[1]];
Y[t_Real] :=
BSplineFunction[L][Max[Min[1, if[(t + Pi)/(2 Pi) sMax]], 0]][[2]];
{XFT, YFT} = Fourier[Table[#[N@t], {
t, -Pi + epsilon, Pi - epsilon, (2 Pi - 2 epsilon)/myu
}]] & /@ {X, Y};
{type, 2 Pi/Sqrt[myu]*((
Transpose[Table[{Re[#], Im[#]} &[Exp[I k Pi] #[[k + 1]]],
{k, 0, nMax}]] & /@ {XFT, YFT}
))}
]
Options[definition3] = {"MaxOrder" -> num, "OpenClose" -> 0.025};
[/mathematica]

3.
[mm collapse="1"]
definition3[pointLists_, OptionsPattern[]] :=
Monitor[Table[
definition2[pointLists[[k]],
If[Head[#] === List, #[[k]], #] &[OptionValue["MaxOrder"]],
If[Head[#] === List, #[[k]], #] &[OptionValue["OpenClose"]]
], {k, Length[pointLists]}
],
Grid[
{{Text[
Style["progress calculating Fourier coefficients", Darker[Green, 0.66]]
],
ProgressIndicator[k/Length[pointLists]]
}},
Alignment -> Left, Dividers -> Center
]
] /;
Depth[pointLists] === 4
[/mm]

4.
[mm collapse="1"]
definition4[{"Closed" | "Open", {{cax_, sax_}, {cay_, say_}}}, t_, n_] :=
{Sum[
If[k == 0, 1/2, 1] cax[[k + 1]] Cos[k t] + sax[[k + 1]] Sin[k t], {k, 0, Min[n, Length[cax]]}
],
Sum[If[k == 0, 1/2, 1] cay[[k + 1]] Cos[k t] + say[[k + 1]] Sin[k t], {k, 0, Min[n, Length[cay]]}
]
}
[/mm]

5.
[mm collapse="1"]
definition5[n_] :=
Show[{
ParametricPlot[
Evaluate[definition4[#, t, n] & /@ fCs], {t, -Pi, Pi}, Axes -> False
]
}]
[/mm]

1.画像の輪郭を取得

まず、数式に変換したい画像を読み込みます。
[mm]img = Import["フォルダパス又はURL"][/mm]
これで任意の画像をimgと言う文字列に格納しました。
実行して画像が表示されるか確認してください。





次にこの画像の輪郭（edge）を取得します。
[mm]
edgeImage =
Thinning[EdgeDetect[
ColorConvert[
ImagePad[
Image[Map[Most, ImageData[img], {2}]], 20, White],
"Grayscale"
]
]
]
[/mm]
これを実行すると.....

このように、色の境目を抽出してくれます。この時点ではまだ数式で扱えるデータではありません。
次のステップで白黒画像に"座標"という概念を与え、境界線を"点群"で表します。

2.輪郭を点群で表す

下のコードを実行すると境界線が2次元直交座標系の点群で表されます。
[mm]edgePoints = {#2, -#1} & @@@ Position[ImageData[edgeImage], 1, {2}][/mm]

これがMathematicaの出力結果ですが使う画像が複雑であればあるほど点は多くなります。
今回使っているTokyo2020のロゴでは2万2千個以上の点が出力されました。




次に今求めた点と点を繋ぎ線にしていきます。
[mm]
SeedRandom[2];
hLines = definition1[edgePoints, 16];
Graphics[{ColorData["DarkRainbow"][RandomReal[]], Line[#]} & /@ hLines]
[/mm]
これを実行すると線で結ばれますが、数式ではありません。
＊definition1[]は独自関数なので準備のコードで定義してください。

これを任意曲線とし、次のステップでフーリエ変換を施します。

3.点群にフーリエ変換を施す

準備で定義した独自関数をバリバリ使ってフーリエ変換します。
[mm]
num = 200;
fCs = definition3[hLines];
definition5[num]
[/mm]

これで画像を数式で近似することができました。

ここで注目してもらいたいのは１行目のnumの値、今は200が代入されていますね。
numは"フーリエ変換後の数式の次数"を表しています。
このnumの値がフーリエ変換の精度になっており、大きいほど元の画像に近づきます。




numを1から増やしていくのをgif画にしたのが最初にもお見せしたこれ

最初は円や曲線なのにだんだんと形になっていく様子がとても美しいです。

まとめ

今回紹介したWolfram社の数式処理システム"Mathematica"
興味を持ち、私も使ってみたい！と思ってくれる方がいるかもしれません。

大学数学の複雑な計算ですらスラスラ解いてくれてとても便利なのですが
このソフト、実はライセンスがとても高価なのです。

一番安いプラン(個人の趣味用、商用利用不可)で年間５万円
企業向けの一番高いプランで年間100万円以上
興味があるからと言う理由だけでライセンス取得するのはすこし抵抗があるかもしれません。

もし、この画像を数式近似してgif画にしてほしい！と言うのがありましたら
コメントやメールなどで送ってください！
私がMathematicaで処理したgif画を添付して返信します！

mail : info@lonelynote.net

この美しさに共感してくれる方が少しでもいたら嬉しいです。
ありがとうございました。