Wolframコンファレンスに登壇しました！

2018年12月18日(火)、Wolframコンファレンスが行われ、僕が登壇しました。
そのスライドを記事化して公開します。

はじめに

Mathematicaは ListPlot, ArrayPlot, Graphics, ListPlay, Export など、画像処理や音声処理の関数が充実しています。
それに加えて数式処理の強みを活かし、画像→音声などの変換にも強いのが特徴です。
本講演では5つの事例をプログラムを通して紹介します。

お品書き

• 物智数
• bion
• 数式曲線
• 解析信号
• フラクタル音楽

自己紹介

僕は岩淵 勇樹 (IWABUCHI Yu(u)ki)と申します。

• 金沢大学自然科学研究科修了
• 博士（工学）
• 面白法人カヤック HTMLファイ部・人事部
• 博士論文は「図形と音声の変換方法とその応用に関する研究」

（このプロフィール画像もMathematicaを使って合成しました。）

数システム → 画像（その1）

Mathematicaによる可視化

Wolfram Mathematicaは数学的対象をプロットしたりグラフィックス表示するのが簡単です。
それにより、研究の効率を上げることができます。

物智数

2進数や10進数などの「記数法」のひとつです。
発表者が2006年に発見しました。

参考:
物智数（Butchi number, 平面的2進数）
岩淵勇樹 卒業論文「2進数の平面的表現に関する考察」 (2007)

2進数変換の再発明

2進数は IntegerDigit でリストとして得ることがきますが、このように多項式を使って実装することもできます。

binaryFromPolynomial[n_] := Reverse[
  CoefficientList[
    FixedPoint[
      Expand[(# - PolynomialMod[#, 2])/2 p + PolynomialMod[#, 2]] &,
      n
    ],
    p
  ]
]

In[1]:= binaryFromPolynomial[5]
Out[1]= {1, 0, 1}

In[2]:= binaryFromPolynomial[p^2 + 1]
Out[2]= {1, 0, 1}

2進数の「2」を「p」に置き換えて多項式表現（ex: 5 = Subscript[101, 2]= 2^2+1 = p^2+1）しています。
これは数式処理ソフトならではの実装といえるでしょう。

2進数の平面拡張

物智数は「左に繰り上がる」という2進数を拡張し、「左と上に繰り上がる」という繰り上がり則にした記数法です。

先ほど再発明した関数の、pの多項式をp, qの多項式に置き換えることにより、物智数を得ることができます。

butchiNumberFromPolynomial[n_?PolynomialQ] := butchiNumber[
  CoefficientList[
    FixedPoint[
      Expand[(# - PolynomialMod[#, 2])/2 (p + q) + PolynomialMod[#, 2]] &,
      n
    ],
    {q, p}
  ]
]

物智数は配列をグリッド配置して右と下に線を入れる決まりとしているため、Formatを使ってすぐに表示できるようにします。

Format[butchiNumber[li_]] := Grid[Reverse[li, {1, 2}], Dividers -> -1 -> True]

In[3]:= butchiNumberFromPolynomial[4]
In[4]= First[%]

Out[3]=
  0 1 1
  1 0 0
  1 0 0

Out[4]= {{0, 0, 1}, {0, 0, 1}, {1, 1, 0}}

1→黒、0→白としてプロットする butchiNumberPlot も定義します。

butchiNumberPlot[butchiNumber[li_]] := ArrayPlot[Reverse[li, {1, 2}]]

In[5]:= Table[butchiNumberPlot[butchiNumberFromPolynomial[n]], {n, 10}]

数を増やすとこうなります。最初これを見たときは滅茶苦茶興奮しました。（当時はC言語による実装でしたが。）

In[6]:= Manipulate[butchiNumberPlot[butchiNumberFromPolynomial[n]], {{n, 10}, 1, 1000}]

セルオートマトンを用いた実装

butchiNumberCellularAutomaton[n_] := butchiNumber[
  FixedPoint[
    CellularAutomaton[{
      Mod[#[[1]], 2] + Quotient[#[[2]], 2] + Quotient[#[[3]], 2] &,
      {},
      {{0, 0}, {-1, 0}, {0, -1}}}, {#, 0}, {{{2}}}
    ] &,
    {{n}}
  ]
]

お蔵入りとはなりましたが、セルオートマトンを利用した実装なんかはCやJavaなどよりは短いコードで書けるので、こういったコードを気軽に(?)使えるのもMathematicaの強みでしょう。
数式処理と豊富な関数により、実装の選択肢が増えるので、適材適所のプログラミングができます。

Mathematicaを使った可視化による研究のメリット

対話的なプログラミングにより、簡潔な関数定義と1行レベルの命令で実験可能となります。
即座に可視化結果を確認できるので試行錯誤しやすいところが特長です。

このように、物智数そのものやそれに関するグラフをプロットしたりして、研究が捗ります。

数システム → 画像（その2）

bion

音楽における「倍音」をイメージした数システムです。
螺旋状に1ずつ加算され、放射状に2倍されます。

複素数でグラフィックス処理

x座標を実部、y座標を虚部として扱います。リスト操作や幾何学的変換が{x, y}のリストに比べて簡便で、個人的にはほとんどこの方式を用いています。

In[1]:= ReIm[3 + 4 I]
Out[1]= {3, 4}

In[2]:= Complex @@ {3, 4}
Out[2]= 3 + 4 I

In[3]:= (-1)^(4/5 Range[0, 5]);
In[4]:= ReIm /@ %;
In[5]:= Graphics[Line[%]]

詳しくは、複素数の魔法のエントリーをご覧ください。

bionの作品化

bionは研究の一環でありつつ、CG作品としても出力しました。
Exportによるグラフィックス出力はラスター形式もベクトル形式もサポートしています。

In[6]:= n = 2048;
In[7]:= gBion = Graphics[Table[
  {
    Hue[0.57, (1 - Log[2, k]/Log[2, n])^2, 1],
    Disk[ReIm[bion[k]], 0.5 k^-0.25]
  },
  {k, n, 1, -1}
]]

In[8]:= Export["bion.eps", gBion]
Out[8]= "bion.eps"

waon

bionの関連作品です。RGB, HSB, CMYKなど表色系も充実しています。

waon[n_] := E^(2 \[Pi] I Log[2, n])

In[9]:= n = 256;
In[10]:= gWaon = Graphics[Table[
  {
    Hue[0.58 + Log[2, k], 1 - (Log[2, k]/Log[2, n])^2, 1], 
    Disk[ReIm[waon[k]], 0.3 k^-0.2]
  },
  {k, n, 1, -2}
]]

In[11]:= Export["waon.eps", gWaon]
Out[11]= "waon.eps"

映像 → 数式

数式曲線

Wolfram|Alphaの “Hatsune Miku curve” がネットで話題になりました。

Pikachu curveなんかもあります。
（余談ですが、ポケモンキャラの知識データベースまであることをつい先日知りました。）

イラストでできるならアニメでもできる！と思い立ち、3年前、アニメーション版の『Bad Apple Curve!!』を制作しました。
スライド

こちらのブログに作り方が載っていますが、あまり参考にせず全部自力で実装しました。

そして、モノクロのアニメーション版ができたので今度はカラー版もと、今回のコンファレンスの目玉として制作しました。
（11月後半からこのデモ制作にかなり時間を費やしてしまい、スライド制作時間がかなり圧迫されてしまいました…。）

でびどる曲線！

今期放送中のアニメ『でびどる！』OPがメインカラーが絞られていて色を抽出しやすく、さらに画像の利用を快く許可していただけたので、コンファレンスでも堂々と上映させていただきました。
『でびどる！』制作のコチンPa!製作委員会ならびにまさたかPに大感謝！
（ちなみに、コンファレンスの週に最終回を迎え、本記事投稿の今週に最終放送の総集編があります。）

（↑ 静止画キャプチャ、GIFアニメーション（音声、BGMはなし）をSNS等で無断アップロードをぜんぜん禁じません。ですって！！！！）

例えばこのコマが・・・

以下のように変身します。

動画で見ると・・・

↑はフル尺ですが、記事に貼る都合上カクカクです。12fpsのもうちょっとヌルヌルした動画は↓のツイートで見られます。

アニメを1フレームずつ処理して、領域ごとに曲線を抽出し、それをsinだけでできた数式に変換しました。

フーリエ記述子

「でびどる曲線！」の原理は、曲線の座標をフーリエ変換し、低周波成分のみを拾う「フーリエ記述子」の理論に基づいています。
今回は何種類かあるフーリエ記述子のうち、G型のフーリエ記述子を採用しました（↑の例だとP型なのでちょっと違います）。

参考文献:
岩淵勇樹・秋田純一 「非周期的な画素配置を持つCMOSイメージセンサの基礎検討」, 電子情報通信学会技術研究報告 (2007)

でびどる曲線！ 制作過程

実際使ったコードを整理し、一連の流れを抽出しました。

画像の読み込み

In[1]:= img = Import["aira.png"]

二値化

binaryImageAira[image_] := ColorNegate[Binarize[ColorDistance[image, Red], 0.5]]

In[2]:= binImg = binaryImageAira[img]

ノイズの除去

In[3]:= gcImg = GeodesicClosing[binImg, 1]

領域分割

In[4]: mcImg = MorphologicalComponents[gcImg];
In[5]:= mcImg // Colorize

レイヤー化

In[6]:= layerLi = Image /@ Table[Unitize[mcImg - n], {n, 1, Max[mcImg]}]

細かい領域を排除

In[7]:= selLayerLi = Select[layerLi, Total[Flatten[1 - ImageData[#]]] > 200 &]

境界追跡

findBoundary[img_] := 
Module[{imgData = 
  ImageData[Binarize[ImagePad[ImagePad[img, -2], 2, 1]]],
  d = ( {
    {0, 1},
    {1, 1},
    {1, 0},
    {1, -1},
    {0, -1},
    {-1, -1},
    {-1, 0},
    {-1, 1}
  } ),
  vec = 0,
  nextPointLi, offset, boundary},
  boundary = Position[imgData, 0, 2, 1];
  While[Length[boundary] == 1 || First[boundary] != Last[boundary],
    nextPointLi = 
    Select[
      RotateLeft[d, vec - 3],
      Extract[imgData, Last[boundary] + #] == 0 &
    ];
    If[nextPointLi == {}, Break[]];
    offset = First[nextPointLi];
    AppendTo[boundary, Last[boundary] + offset];
    vec = First[First[Position[d, offset]]] - 1;
  ];
  Return[Most[Reverse /@ boundary]];
]

線分の座標列を算出

In[8]:= lineLi = Cases[findBoundary /@ selLayerLi, Except[{}]]
In[9]:= Graphics[Line /@ lineLi]

Out[8]= {{
  {247, 3}, {248, 3}, {249, 3}, {248, 4}, {247, 5},
  {247, 6}, {246, 5}, {245, 6}, {244, 7}, {243, 8}, ...
}}

離散フーリエ変換

fft[li_] := Fourier[Complex @@ # & /@ li, FourierParameters -> {-1, 1}]

In[10]:= fftLi = fft /@ lineLi

Out[10] = {{
  337.995 + 151.436 I, -36.9354 + 8.77659 I, 35.2291 - 11.0688 I, 
  8.47604 - 3.34839 I, -22.2869 + 7.85072 I, 11.6755 - 11.2943 I, ...
}}

sinの係数・位相のリストを算出

atan[0, b_] := Pi/2
atan[0., b_] := N[Pi/2]
atan[a_, b_] := ArcTan[b/a] + Piecewise[{{0, a >= 0}}, Pi]

(* reSin・imSinの係数と位相だけのリスト *)
reSinComponent[c1_, c2_] := With[
  {p = -Im[c1] + Im[c2], q = Re[c1] + Re[c2]},
  {Sqrt[p^2 + q^2], atan[p, q]}
]
imSinComponent[c1_, c2_] := With[
  {p = Re[c1] - Re[c2], q = Im[c1] + Im[c2]},
  {Sqrt[p^2 + q^2], atan[p, q]}
]

(* 係数・位相リストのリスト *)
reComponents[components_] := With[
  {n = Length[components]/2}, 
  Join[
    {Re[First[components]]},
    Table[reSinComponent[components[[k + 1]], components[[-k]]], {k, 1, n}]
  ]
]
imComponents[components_] := With[
  {n = Length[components]/2},
  Join[
    {Im[First[components]]}, 
    Table[imSinComponent[components[[k + 1]], components[[-k]]], {k, 1, n}]
  ]
]

In[11]:= reComp = reComponents /@ fftLi;
In[12]:= imComp = imComponents /@ fftLi;
In[13]:= Short[First[reComp]]
In[14]:= Short[First[imComp]]

Out[13] = {337.995, {130.314, 3.87657}, ... , {0.00670017, 1.5708}}
Out[14] = {151.436, {80.2655, -1.40162}, ... , {7.25102*10^-16, 1.5708}}

↑最初の要素だけ定数なので位相なしです。

数式に変換

sinFromComponent[{k_, θ_}, n_] := k Sin[n t + θ]
reExpressionList[n_] := First[#] + Sum[
  sinFromComponent[PadRight[#, n, {{0, 0}}][[m]], m - 1],
  {m, 2, n}
] & /@ reComp;
imExpressionList[n_] := First[#] + Sum[
  sinFromComponent[PadRight[#, n, {{0, 0}}][[m]], m - 1],
  {m, 2, n}
] & /@ imComp;

In[15]:= reExprLi = reExpressionList[10];
In[16]:= imExprLi = imExpressionList[10];
In[17]:= First[reExprLi]
In[18]:= First[imExprLi]

Out[17] = 337.995 + 130.314 Sin[3.87657 + t] +
  13.6335 Sin[4.67236 + 2 t] + 8.17294 Sin[2.51724 + 3 t] +
  16.0637 Sin[3.44023 + 4 t] + 2.02671 Sin[2.66355 + 5 t] +
  12.9663 Sin[0.334349 + 6 t] + 11.9775 Sin[4.08646 + 7 t] +
  5.31037 Sin[1.20734 + 8 t] + 5.20582 Sin[3.13833 + 9 t]
Out[18] = 151.436 - 33.1266 Sin[0.827451 - 5 t] -
  18.0512 Sin[0.830595 - 3 t] - 87.0867 Sin[0.263506 - 2 t] -
  80.2655 Sin[1.40162 - t] + 39.8491 Sin[3.13284 + 4 t] +
  8.84287 Sin[1.45426 + 6 t] + 16.8957 Sin[4.56961 + 7 t] +
  6.99099 Sin[3.02382 + 8 t] + 14.3262 Sin[2.39458 + 9 t]

プロットしてみる

In[19]:= ParametricPlot[
  Transpose[{reExprLi, 360 - imExprLi}],
  {t, 0, 2 Pi},
  PlotRange -> {{0, 640}, {0, 360}}, 
  AspectRatio -> Automatic
]

次数を変えながらコマをプロット

In[20]:= Manipulate[
  ParametricPlot[
    Evaluate[Transpose[{reExpressionList[n], 360 - imExpressionList[n]}],
    {t, 0, 2 Pi}, 
    PlotRange -> {{0, 640}, {0, 360}}
  ],
  {{n, 50}, 2, 150, 1}
]

これをさらに、動画シーケンスに適用して大量処理するコードを実装して、書き出された連番画像を結合してアニメーション化しました。

音声 → 曲線

解析信号

解析信号とは、簡単に言えば音声信号を複素信号化してプロットしたものです。
周期信号に対して閉曲線が描かれます。

解析信号フィルタ

負の周波数成分をカットするフィルタをFFT後のスペクトルに対して適用することで、解析信号が得られます。

参考文献:
岩淵勇樹，秋田純一，北川章夫 「閉曲線を利用した音色の操作法」 芸術科学会論文誌 (2011)
岩淵勇樹 博士論文 「図形と音声の変換方法とその応用に関する研究」 (2012)

解析信号の実演

解析信号フィルタの正体は、ヘビサイドのステップ関数です。

step[n_Integer?Positive] := 
 1 + Sign[Sin[\[Pi] Range[1, -1, -2/n] // Most]]

In[1]:= step[10]
Out[1] = {1, 2, 2, 2, 2, 1, 0, 0, 0, 0}

「2」の部分が正の周波数成分を拾う部分（都合上1倍ではなく2倍になってます）、「0」の部分が負の周波数成分をカットする部分です。

analyticSignalTransform[li_List] := InverseFourier[Fourier[li]*step[Length[li]]]

離散フーリエ変換された信号は、前半が正の周波数成分、後半が負の周波数成分であることに注意してください。

例として、フルートの音を使います。

In[2]:= snd = ExampleData[{"Sound", "Flute"}]

音声のデータ部分、サンプルレート、サンプル数を抽出します。

In[3]:= sndData = snd[[1]][[1]]
In[4]:= sampleRate = snd[[1]][[2]]
In[5]:= len = Length[sndData]

Out[3] = {0.00244148, 0.00332652, 0.00283822, 0.00357067 ... }
Out[4] = 22050
Out[5] = 57920

解析信号化します。

In[6]:= asLi = analyticSignalTransform[sndData]

Out[6] = {0.00244148 + 0.0137746 I, 0.00332652 + 0.0105611 I, 
  0.00283822 + 0.0104776 I, 0.00357067 + 0.00923803 I, ... }

この複素数をよく見ると、実部がsndDataになってます。

全体をプロットすると、塗りつぶされたような形が出てきます。

In[7]:= ListLinePlot[ReIm /@ asLi, PlotRange -> All, AspectRatio -> Automatic]

20000サンプル目から500サンプルだけを抽出してプロットします。

In[8]:= ListLinePlot[ReIm /@ Take[asLi, {20000, 20500}], PlotRange -> All,  AspectRatio -> Automatic]

これがフルートの音色の形です。

解析信号のアニメーション

ListAnimate でアニメーションさせます。

In[9]:= fps = 30;
In[10]:= totalSecond = N[len/sampleRate];
In[11]:= totalFrame = Floor[len*fps/sampleRate];
In[12]:= samplePerFrame = Floor[sampleRate/fps];
In[13]:= ListAnimate[
  Table[ListPlot[
    Take[
      Transpose[{Re[asLi], Im[asLi]}],
      {frame*samplePerFrame + 1, (frame + 1)*samplePerFrame}
    ],
    Joined -> True, 
    PlotRange -> {{-1, 1}, {-1, 1}}, AspectRatio -> Automatic, Axes -> False, 
    PlotStyle -> Black
  ],
  {frame, 0, totalFrame - 1}],
  fps
]

さまざまな解析信号

楽器や声質など、音色に特有の形が現れます。

参考:
解析信号ギャラリー

図形 → 音楽

フラクタル音楽

「フラクタル音楽」とは、フラクタル図形を可聴化するアート活動です。
TwitterとInstagramで約1年間毎日投稿、計459作品のフラクタル音楽を作りました。

カーペットフラクタル(@carpet_fractal) - Instagram

フラクタル音楽制作の流れ

1. フラクタル図形の描画
2. 音声信号に変換
3. 解析信号で可視化

参考:
岩淵勇樹 「カーペットフラクタルの可聴化と音楽生成」, 情報処理学会第105回音楽情報科学研究会 (2014)

フラクタル生成から可聴化まで

簡単な例として、「シェルピンスキーのカーペット」を可聴化します。

iterate[x_] := ArrayFlatten[x /. {
  0 -> ( {
    {0, 0, 0},
    {0, 0, 0},
    {0, 0, 0}
  } ),
  1 -> ( {
    {1, 1, 1},
    {1, 0, 1},
    {1, 1, 1}
  } )
}]

1回のイテレーションで3×3の配列が得られます。

In[1]:= iterate[1] // Grid

Out[1] =
1 1 1
1 0 1
1 1 1

2回のイテレーションで9×9の配列が得られます。

In[2]:= iterate[iterate[1]] // Grid

Out[2] = 
1 1 1 1 1 1 1 1 1
1 0 1 1 0 1 1 0 1
1 1 1 1 1 1 1 1 1
1 1 1 0 0 0 1 1 1
1 0 1 0 0 0 1 0 1
1 1 1 0 0 0 1 0 1
1 1 1 1 1 1 1 1 1
1 0 1 1 0 1 1 0 1
1 1 1 1 1 1 1 1 1

イテレーションを繰り返し行うために Nest を使います。

In[3]:= Nest[iterate, 1, 2];

Out[3] = 
1 1 1 1 1 1 1 1 1
1 0 1 1 0 1 1 0 1
1 1 1 1 1 1 1 1 1
1 1 1 0 0 0 1 1 1
1 0 1 0 0 0 1 0 1
1 1 1 0 0 0 1 1 1
1 1 1 1 1 1 1 1 1
1 0 1 1 0 1 1 0 1
1 1 1 1 1 1 1 1 1

白黒で描画するとこうなります。

In[4]:= Nest[iterate, 1, 2] // ArrayPlot

イテレーションを数回繰り返すとこのように多数の穴が空いた図形が得られます。

In[5]:= Manipulate[ArrayPlot[Nest[iterate, 1, n]], {n, 1, 6, 1}]

可聴化に、5回のイテレーションを使います。

In[6]:= carpet = Nest[iterate, 1, 5];
In[7]:= ArrayPlot[carpet, Frame -> False, PixelConstrained -> 1]

可聴化は非常にシンプルです。2次元配列を単純に1次元にする Flatten を通した後、リストを音声として聞く ListPlay を適用するだけです。

In[8]:= ListPlay[Flatten[carpet]]

バリエーションを増やすために、乱数を用います。

randomIterate[x_] := ArrayFlatten[x /. {
  0 -> ( {
    {0, 0, 0},
    {0, 0, 0},
    {0, 0, 0}
  } ),
  1 -> iterationTbl
}]

In[9]:= SeedRandom[0]
In[10]:= iterationTbl = RandomInteger[{0, 1}, {3, 3}];
In[11]:= Grid[iterationTbl]

Out[11] =
1 0 1
0 0 1
1 1 0

In[12]:= randomCarpet = Nest[randomIterate, 1, 5];
In[13]:= ArrayPlot[randomCarpet, Frame -> False, PixelConstrained -> 1]

In[14]:= ListPlay[Flatten[randomCarpet]]

さらにバリエーションを増やすために、平行移動してビット演算をしてみます。

In[15]:= shiftCarpet = BitXor[carpet, RotateLeft[carpet, {20, 45}]];
In[16]:= ArrayPlot[
  shiftCarpet,
  PixelConstrained -> 1,
  Frame -> False, 
  ColorFunction -> (Blend[{LCHColor[0.2, 1, 1.8], White}, #] &)
]

In[17]:= ListPlay[Flatten[shiftCarpet]]

このようにして週替りでプログラムを改変し、無数のフラクタル音楽を作ることができました。

まとめ

• 「数システム → 画像」「映像 → 数式」「音声 → 曲線」「図形 → 音楽」の変換について解説
• 数式を介して画像・音声間を柔軟に行き来可能
• Wolfram Mathematicaは研究にも芸術活動にも最適

自分にとってはMathematica人生こそ研究人生といっても過言ではないので、今回の登壇で５つの研究を紹介できたことは卒業研究以降の数学研究のひとつの集大成となりました。
Wolfram社ならびにウルフラムリサーチアジアリミティッドの皆さんに、貴重な機会をいただけたことを感謝いたします！