プログラミングの備忘録

プログラムをつくる過程を残すもの

タッパーの自己言及式の謎を解く

こんにちは。
今回は「タッパーの自己言及式 (Tupper's self-referential formula)」を取りあげようと思います。
(記事タイトルを AI が生成してくれるようになったので、さっそく使ってみました。)


例のごとく、少し前にこんなツイートを見かけました。

(Fermat's Library は理系的な雑学のツイートが多く、知見が広がるのでよく見ています。)

どこかで見たことがあるという方もいるかもしれません。
数式を図にすると数式そのものが浮かび上がる。
これは掘り下げない訳にはいきません。


目次


まずは実装してみる

さっそくですが、その真偽を確かめるべく実装してみます。


先程のツイートには載っていませんでしたが、$k$ という値が必要になります。

詳しくはまた後で触れますが、$k$ の値は…

$$ k = 960939379918958884971672962127852754715004339660129306651505519271702802395266424689642842174350718121267153782770623355993237280874144307891325963941337723487857735749823926629715517173716995165232890538221612403238855866184013235585136048828693337902491454229288667081096184496091705183454067827731551705405381627380967602565625016981482083418783163849115590225610003652351370343874461848378737238198224849863465033159410054974700593138339226497249461751545728366702369745461014655997933798537483143786841806593422227898388722980000748404719 $$

543桁にもなる大きな数です。


図は、ある座標 $(x, y)$ について右辺を計算し、その計算結果が $\displaystyle \frac{1}{2}$ と比べて大きいか小さいかで色分けをすることによって描くことができます。
(そのまま図示すると上下左右が逆になった図になるため、少し調整が必要になります。)


また、$\lfloor a \rfloor$ は $a$ を超えない最大の整数、$\mathrm{mod} (a, b)$ は $a$ を $b$ で割った余りを指します。


ここでは、Processing と Python で実装してみようと思います。

Processing

巨大数を扱うにあたって

実は、素直に実装しようとしても、いつも通りやっていてはうまくいきません。

普段は全く気にしませんが、扱える数値には下限上限があります。
int型での上限は $2147483647$ とされており、到底543桁には及びません。
(int型だ、と宣言するとある程度のデータ量をその数値に割くようにしてくれて、その「ある程度」が関係するとか)

long型という、もっと大きな (小さな) 数を扱えるものもありますが、それでも上限は $9223372036854775807$ と足りません。


そこで使うのが、「BigDecimal型 (BigDecimal (Java Platform SE 8))」です。
上限は $2147483647$ くらいだとされているようで、とんでもなく大きな数まで扱えることがわかると思います。

そんなBigDecimal型は Java のライブラリで、

import java.math.BigDecimal;

で読み込むことができます。

実装

あとは (関数は BigDecimal のものを使わなければいけませんが) 普通に計算していけば実装できます。

import java.math.BigDecimal;

String strk = "960939379918958884971672962127852754715004339660129306651505519271702802395266424689642842174350718121267153782770623355993237280874144307891325963941337723487857735749823926629715517173716995165232890538221612403238855866184013235585136048828693337902491454229288667081096184496091705183454067827731551705405381627380967602565625016981482083418783163849115590225610003652351370343874461848378737238198224849863465033159410054974700593138339226497249461751545728366702369745461014655997933798537483143786841806593422227898388722980000748404719";
BigDecimal k = new BigDecimal(strk);
int a = 17; // 縦
int b = 106; // 横

void setup(){
  size(800, 200);
  background(255);
  
  float l = width / b; // セル1辺の長さ
  
  for(int i = 0; i < b; i ++){
  for(int j = 0; j < a; j ++){
    BigDecimal x = new BigDecimal(i);
    BigDecimal y = new BigDecimal(j).add(k);
    
    BigDecimal v1 = new BigDecimal(a).multiply(x);
    BigDecimal v2 = y.remainder(new BigDecimal(a));
    BigDecimal v3 = v1.add(v2);
    BigDecimal v4 = new BigDecimal(2).pow(v3.intValue());
    BigDecimal v5 = new BigDecimal(1).divide(v4);
    BigDecimal v6 = y.divideToIntegralValue(new BigDecimal(a));
    BigDecimal v7 = v6.multiply(v5);
    BigDecimal v8 = v7.remainder(new BigDecimal(2));
    
    if(v8.intValue() > 0.5){
      fill(0);
    }
    else{
      fill(255);
    }
    rect(width-i*l-2*l, j*l+2*l, l, l);
  }
  }
  print("end"); // 計算が終わったら「end」と出力
}


変数に v シリーズがいくつも出てきてややこしかったり、描画位置が微妙な位置になっていて気持ち悪かったりしますが、実行すると一応描画できました。
(終了までに20秒程度かかりました。)

Python

Python では、int型の上限はなんとなし
しかし、float型は上限 $308$ 桁くらいと足りません。
(それでも十分大きいですが…)

結局、「decimal (decimal --- 十進固定及び浮動小数点数の算術演算 — Python 3.12.1 ドキュメント)」をインポートして巨大数を扱うことになります。

from decimal import Decimal, getcontext
import matplotlib.pyplot as plt

def tupper(k, a, i, j): # タッパーの自己言及式
    x = i
    y = j + k

    t = int(((Decimal(y) // 17) * (2 ** Decimal(- a * x - y % a))) % 2)

    if t > 0.5:
        return 1
    else:
        return 0

getcontext().prec = 1000 # 精度

k = 960939379918958884971672962127852754715004339660129306651505519271702802395266424689642842174350718121267153782770623355993237280874144307891325963941337723487857735749823926629715517173716995165232890538221612403238855866184013235585136048828693337902491454229288667081096184496091705183454067827731551705405381627380967602565625016981482083418783163849115590225610003652351370343874461848378737238198224849863465033159410054974700593138339226497249461751545728366702369745461014655997933798537483143786841806593422227898388722980000748404719
a = 17 # 縦
b = 106 # 横

print('a = {}, b = {}'.format(a, b))

xx = []
yy = []
for x in range(0, b):
    for y in range(0, a):
        if tupper(k, a, x, y) == 1:
            xx.append(-x)
            yy.append(-y)

# 散布図の描画
fig = plt.figure()
ax = fig.add_subplot()
ax.scatter(xx, yy, s=25, c="black", marker=",")
plt.show()


ネタばらし

改めて、タッパーの自己言及式を見てみます。

$$ \frac {1}{2} < \left\lfloor \mathrm{mod} \left( \left\lfloor {\frac {y}{17}} \right\rfloor 2^{-17 \lfloor x \rfloor - \mathrm{mod} \left( \lfloor y \rfloor ,17 \right)}, 2 \right) \right\rfloor $$

この数式では何が起こっているのでしょうか。
(いろいろありまして、先程つくったコードと以下の説明で使っている文字が一部対応していません。悪しからず。)


ざっくり言えば、この数式は $k$ をグラフ上に図示する機械で、$k$ は図示するとその数式そのものの形になるように調整されたもの、です。
つまり数式が先にあって、$k$ をうまいことつくって自己言及的になるようにした、という感じです。


さて、ある2進数 $a$ を $2^{0}$ の位から座標平面に $(0, 0), (0, 1), (0, 2), \dots$ と $h$ 個並べ、$(1, 0)$ に戻ってきてまた縦に並べていく、という操作を繰り返し、縦 $h$、横 $w$ の四角形をつくることを考えます。
($a$ をうまく調整して、この四角形 (ビットマップ) が何かしらの絵 (上下左右が逆) になるようにしておきます。)

$a$ の $n$ 桁目 ($2^{0}$ の位を $n = 0$ とする) の数は

$$ \left\lfloor \mathrm{mod} \left( a \times 2^{-n}, 2 \right) \right\rfloor $$

で求められ、これが $0$ か $1$ かを判定したいので $\displaystyle \frac{1}{2}$ を基準として、

$$ \frac{1}{2} < \left\lfloor \mathrm{mod} \left( a \times 2^{-n}, 2 \right) \right\rfloor $$

ここで、四角形における位置 $(i, j)$ が $a$ の何桁目にあたるかを求めるには、$n = h \times i + j$ を計算すれば良いので、

$$ \frac{1}{2} < \left\lfloor \mathrm{mod} \left( a \times 2^{- h \times i - j}, 2 \right) \right\rfloor $$


座標からその位置でのビットを計算したいので、数式から $a$ を消去します。
これは数学的なトリック (?) なのですが、$(i, j)$ から $(x, y) = (i, j + k)$ ($k = h \times a$) に平行移動させます。
すると、$\displaystyle a = \left\lfloor \frac{y}{h} \right\rfloor$、$j = \mathrm{mod} (y, h)$ と表せるので、

$$ \frac{1}{2} < \left\lfloor \mathrm{mod} \left( \left\lfloor \frac{y}{h} \right\rfloor \times 2^{- h \times x - \mathrm{mod} (y, h)}, 2 \right) \right\rfloor $$


$x$、$y$ は実数なので、四角形に描画するために整数に変換して、

$$ \frac{1}{2} < \left\lfloor \mathrm{mod} \left( \left\lfloor \frac{y}{h} \right\rfloor \times 2^{- h \times \lfloor x \rfloor - \mathrm{mod} (\lfloor y \rfloor, h)}, 2 \right) \right\rfloor $$


以上で、一般化した式が求められました。
タッパーの自己言及式では縦 $h$ を $17$ にしているため、初めに挙げたような式となります。
横 $w$ は式に含まれず、2進数 $a$ の桁数を $d$ として、$\displaystyle \left\lceil \frac{d}{h} \right\rceil$ で求めることができます。


ということは…?

$$ \frac {1}{2} < \left\lfloor \mathrm{mod} \left( \left\lfloor {\frac {y}{17}} \right\rfloor 2^{-17 \lfloor x \rfloor - \mathrm{mod} \left( \lfloor y \rfloor ,17 \right)}, 2 \right) \right\rfloor $$

この数式は $k$ をグラフ上に図示する機械で、$k$ は図示するとその数式そのものの形になるように調整されたものだった、ということでした。

つまり、$k$ を変えることで、好きなように絵を描くことができるということになります。


これに関して、便利なサイトがありました。

Tupper's Formula Tools

ドット絵を $k$ に変えてくれます。(逆も可)


例えば、k = 37067602316440724456404475121574976829800827292700312151932062052691617906507927381280234914856650067959341707704037256202764341388719745773546243699883854050921017111451676703126563717998971858347996166958438528749534028200461401835215861451434273715380484400478716115151572056421758253407955986722944213998525692367560978044283475116945590335229170642234077735071680796317197511645338572277873548247607374220806113516110376740106816215430532077738073620942940121330994950836541296197937646843827484372199351437787467973352862443913609216 とすると、


こうしてしまうと、自己言及でも何でもない、ただのお絵かきマシーンになってしまいますが。


まとめ

以上、「タッパーの自己言及式」を Processing と Python で実装し、そのカラクリの解明をしてみました。


参考

プログラミング言語「Uiua」を触ってみる

こんにちは。
今回は、プログラミング言語「Uiua」を紹介します。


こんなツイートを見た。


謎の文字の羅列。なぜこれで絵が描けるのか。
こんなん見せられたらやるしかありません。

というわけでやります。


目次


「Uiua」とは

Uiua (Uiua) はプログラミング言語のひとつで、今年出たばかりのできたてほやほや言語です。
発音は「ウィーワー (?)」という感じらしい。

公式サイトには「単純で、美しく、そして暗黙なコードを書くことに重きを置いた」や「可読性を残しつつもできるだけ短いコードを書ける」とありました。


ちょうど最近「Piet」を触ってみましたが、Uiua も同様にスタックベースで、後入れ後出しで値を取り扱います。

taq.hatenadiary.jp

ただ、数字や文字列だけでなく、画像や音の出力も可能なようです。


はサイン関数、 は配列の長さを返す関数というように特定の文字 (Glyph) に関数が組み込まれており、パッと見は難解プログラミング言語ともいえそうです。

組み込み関数は Documentation (https://www.uiua.org/docs) の「Functions」の章に羅列されています。
基本的には、上記サイトの「Tutorial」や「Other Docs」を読めば基礎的な力はつくと思われます。


インストールと実行環境

インストールしなくてもいい

初っ端から矛盾してますが、インストールしなくても公式サイトにエディタ (https://www.uiua.org/pad?src=0_4_1__) があるのでそちらを使うと楽です。

コードが URL に埋め込まれるのか、つくったプログラムを共有することもできるようです。


Uiua interpreter のインストール

基本的には公式サイト (https://www.uiua.org/docs/install) に記載の方法に従えばよいと思います。

しかし、私は Windows ユーザーなのですが、なぜか GitHub からダウンロードした exe ファイルを起動してもインストールされず、結局 Windows ユーザー以外の方法に従ってインストールしました。
(ダウンロードしてきた exe ファイルをダブルクリックすればいいんですよね・・・?)


ざっくり言うと、Rust をインストールして、コマンドプロンプトcargo install uiua を実行すれば完了です。
ただし私の環境では「Visual Studio C++ Build tools」なるものも足りなかったようで、これもダウンロードしました。

うまくいかなかったら、パソコンを再起動するといけるかもしれません。


VSCode での開発環境の構築

「Extensions」のタブで「uiua」と検索すると、VSCode で Uiua のシンタックスハイライトやフォーマッティングができる拡張機能が見つかります。
また、Uiua における関数は特殊文字もあるので、簡単に打ち込めるようにしたキーパッドが使える拡張機能もあります。


プログラムの実行方法

Python などと同じように、コードファイルのあるフォルダでコマンドプロンプトを起動し、

uiua run

uiua ファイル名.ua

と打ち込むことで実行結果が表示されます。


コーディング

コードは基本、右から左に実行されます。


例えば、

+ 2 ⇡ 5


は range 関数で、0 から n-1 までの数を入れた配列をつくります。
つまり、⇡ 5[0 1 2 3 4] となります。

そして + 2 で各要素に 2 を足していきます。

結果、[2 3 4 5 6] となります、


Uiua における組み込み関数を簡単に打ち込めるようになるキーパッドもあると書きましたが、関数の名前自体を書くことでも代用できます。

つまり、先程のコードは

add 1 range 5

と書いても問題ありません。
(実際は、実行時に勝手に Glyph に書き換えられる、という形になっています。)

一意に決められれば、関数名は多少省略してもよいようです。


ここからは、実際にいろいろつくってみます。


フィボナッチ数列

⍥(⊂/+↙2.)10[1 1]
# 結果: [144 89 55 34 21 13 8 5 3 2 1 1]

初めに配列 [1 1] をつくり、() 内の処理を ⍥()10 で 10回繰り返します。

() 内では、. で配列の複製、↙2 で複製した配列の先頭2つを取り出す、/+ で取り出した2つの和を計算する、 で元の配列に計算結果を入れる、という感じです。


素数リストの作成

▽¬∊:♭⊞×...+1↘1⇡20
# 結果: [2 3 5 7 11 13 17 19]

まず、+1↘1⇡20 で 2 から 20 までの自然数を入れた配列をつくり、... で同じ配列を3つ増やします。
は先頭から指定した個数だけ要素を配列から除きます。

⊞× で2つの配列の各要素について積を計算し、2次元の配列をつくります。いわば合成数を集めた配列ですね。
(100マス計算みたいな感じ。)

で2次元配列を1次元にします。扱いやすくしただけです。
(フラットだけに。)

: でスタックの上2つを入れ替え、¬∊ で初めの配列と合成数を集めた配列を比較し、被っている要素は 0、被っていない要素は 1 に変えます。
つまり、初めの配列で合成数は 0、素数は 1 に置き換えられた配列ができあがります。

最後に で、初めの配列と 0 と 1 の配列を比較し、0 の部分は消し、1 の部分は残します。

これで晴れて合成数は除かれ、素数だけが残る、ということになります。


ちなみに、関数をつくることもでき、

f ← ▽¬∊:♭⊞×...+1↘1⇡
f20

で、先程と同じ意味になります。


ウルフラムコード

Uiua では画像もつくることができます。

画像に変換する関数があるわけではなく、30 × 30 以上の2次元配列がつくられると勝手に変換される、という感じのようです。

⇌[⍥(∊:⊚⋯90⍘⋯⇌◫3⊂0⇌⊂0.)⌊÷2⧻.]=⌊÷2:⇡.200

ここまでくると解説も大変です。ざっくりとだけ。

=⌊÷2:⇡.200 で最初の盤面をつくります。
⍥()⌊÷2⧻.() 内の処理を、横幅を2で割った数 (⌊÷2⧻.) 回だけ繰り返します。
() 内では、⊂0⇌⊂0. で両端に 0 を追加 (境界条件)、⍘⋯⇌◫3 で3つずつ切り出して10進数に変換、∊:⊚⋯90 でルール90を適用、という流れです。

ウルフラムコードは過去にセル・オートマトンの回で触れていますが、そのときのルールをそのまま適用した感じなので、良ければそちらも読んでみてください。

taq.hatenadiary.jp


セル・オートマトンに関しては の扱いが重要になってきそうです。


その他、マンデルブロ集合やライフゲームを実装している例もありました。
実は Uiua のロゴも Uiua でつくられているとかいないとか・・・。


まとめ

以上、Uiua を触ってみました。

パッと見は難解ですが、実際やってみると意外と読めるので楽しいです。

gif や音声もつくることができますが、それはまた機会があればということにしておきます。


最後まで読んでいただいてありがとうございました。また次回。

processingで粉遊び的なものをつくってみる

こんにちは。
今回は、processing で「粉遊び」的なものをつくってみます。


「粉遊び」というゲーム自体は、以前にセル・オートマトンを取り上げた際に「ラングトンのアリ」の部分で少し引き合いに出しました。

dan-ball.jp taq.hatenadiary.jp


ダンボールのゲーム自体は昔から好きで、シンプルながら奥深いというか、「棒レンジャー」なんかは割とやりこんだりしていました。
このサイトに載っているようなものをつくってみたいと思って processing に手を伸ばした節もあります。


そんなとき、こんなツイートを見かけました。

まさに「粉遊び」だと思うとともに、以前にセル・オートマトンを少し触ったこともあるので、自分でやってみたいと思いました。

ということで、やります。


目次


前準備

メインの計算部分をつくる前に、いろいろ準備しておきます。

描画画面を整える

格子状に線を引き、「粉」のある部分を白で塗ります。

int n = 100; // 横 or 縦に並べるセルの数
int[][] cell = new int[n][n]; // セルの情報を入れておく配列 (0: 無、1: 粉)
float l;

void setup(){
  size(500, 500);

  l = width/n; // セルの一辺の長さ
}

void draw(){
  background(0);
 
  // グリッド線
  stroke(100);
  for(int i = 0; i < n; i ++){
    line(i*l, 0, i*l, height);
  }
  for(int j = 0; j < n; j ++){
    line(0, j*l, width, j*l);
  }
  
  // 粉の描画
  for(int i = 0; i < n; i ++){
  for(int j = 0; j < n; j ++){
    if(sand[i][j] == 1){ // 粉のある部分を白で塗る
      fill(255);
      rect(i*l, j*l, l, l);
    }
  }
  }

  // メインの計算部分
  // cell = calc(cell);
}


粉を追加できるようにする

draw() 内に以下のコードを書き加えます。
このコードでは、マウスの重なった位置のセルを「粉」のあるセルに変えます。

// 粉の追加
if(mousePressed == true){
  for(int i = 0; i < n; i ++){
  for(int j = 0; j < n; j ++){
    if(i*l < mouseX && mouseX <= (i+1)*l){
    if(j*l < mouseY && mouseY <= (j+1)*l){
      sand[i][j] = 1;
    }
    }
  }
  }
}


まだこの段階では「粉」は動きません。
以降で動かし方をつくっていきます。


素直に実装してみる

それでは、計算部分の作成に移ります。


イントロで挙げたツイートから、16種類の変化をそのままコードに打ち込めばひとまず実装できそうです。

計算部分の処理として関数 calc() を定義して、draw() 内で呼び出すことで実行できるようにしてみました。

他にやり方はあるかもしれませんが、参考までに。

int[][] calc(int[][] cell){
  for(int j = n-2; j > -1; j --){
  for(int i = 0; i < n-1; i ++){
    int a = cell[i][j]; // 左上
    int b = cell[i+1][j]; // 右上
    int c = cell[i][j+1]; // 左下
    int d = cell[i+1][j+1]; // 右下
    
    if(a == 0 && b == 0 && c == 0 && d == 0){
      cell[i][j] = 0;
      cell[i+1][j] = 0;
      cell[i][j+1] = 0;
      cell[i+1][j+1] = 0;
    }
    if(a == 0 && b == 0 && c == 0 && d == 1){
      cell[i][j] = 0;
      cell[i+1][j] = 0;
      cell[i][j+1] = 0;
      cell[i+1][j+1] = 1;
    }
    if(a == 0 && b == 0 && c == 1 && d == 0){
      cell[i][j] = 0;
      cell[i+1][j] = 0;
      cell[i][j+1] = 1;
      cell[i+1][j+1] = 0;
    }
    if(a == 0 && b == 0 && c == 1 && d == 1 ||
       a == 0 && b == 1 && c == 0 && d == 1 ||
       a == 0 && b == 1 && c == 1 && d == 0 ||
       a == 1 && b == 0 && c == 1 && d == 0 ||
       a == 1 && b == 1 && c == 0 && d == 0 ||
       a == 1 && b == 0 && c == 0 && d == 1){
      cell[i][j] = 0;
      cell[i+1][j] = 0;
      cell[i][j+1] = 1;
      cell[i+1][j+1] = 1;
    }
    if(a == 0 && b == 1 && c == 0 && d == 0){
      cell[i][j] = 0;
      cell[i+1][j] = 0;
      cell[i][j+1] = 0;
      cell[i+1][j+1] = 1;
    }
    if(a == 0 && b == 1 && c == 1 && d == 1 ||
       a == 1 && b == 1 && c == 0 && d == 1){
      cell[i][j] = 0;
      cell[i+1][j] = 1;
      cell[i][j+1] = 1;
      cell[i+1][j+1] = 1;
    }
    if(a == 1 && b == 0 && c == 0 && d == 0){
      cell[i][j] = 0;
      cell[i+1][j] = 0;
      cell[i][j+1] = 1;
      cell[i+1][j+1] = 0;
    }
    if(a == 1 && b == 0 && c == 1 && d == 1 ||
       a == 1 && b == 1 && c == 1 && d == 0){
      cell[i][j] = 1;
      cell[i+1][j] = 0;
      cell[i][j+1] = 1;
      cell[i+1][j+1] = 1;
    }
    if(a == 1 && b == 1 && c == 1 && d == 1){
      cell[i][j] = 1;
      cell[i+1][j] = 1;
      cell[i][j+1] = 1;
      cell[i+1][j+1] = 1;
    }
  }
  }
  
  return(cell);
}


まあ煩雑ですね。

2 × 2 のセルを、左上、右上、左下、右下の順に abcd に対応させ、各場合について変化前後での状態をそのまま書き下していったような感じになっています。

これを画面全体について繰り返し適用していくことで、次のステップでの「粉」の位置を決める、という方法です。


実行すると、


ここで、「画面全体に適用する」というときに、

for(int j = n-2; j > -1; j --){
for(int i = 0; i < n-1; i ++){

という、下から考えていくようにしないとうまく動きません。


for(int i = 0; i < n-1; i ++){
for(int j = 0; j < n-1; j ++){

このようにすると、動きはしますが過程が無く、「粉」を追加した途端に一番下まで落ちている、というような形になってしまいます。

(これを解決するのに調べ物をしていたら、自分がやってみたいことを詰め込んだようなものを見つけてしまいました (砂遊びシミュレーション|クリエイティブコーディングの教科書)。)


少し工夫

先程のコードでは煩雑だと思い、改良 (?) してみました。
(あまり変わってない感じもしますが・・・)

int[][] calc(int[][] cell){
  int[][] cn = {{0, 0}, {0, 0}}; // 変換後の 2×2 セル
  
  // 16種類の 2×2 セル
  int[][] c0000 = {{0, 0}, {0, 0}};
  int[][] c0001 = {{0, 0}, {0, 1}};
  int[][] c0010 = {{0, 0}, {1, 0}};
  int[][] c0011 = {{0, 0}, {1, 1}};
  int[][] c0100 = {{0, 1}, {0, 0}};
  int[][] c0101 = {{0, 1}, {0, 1}};
  int[][] c0110 = {{0, 1}, {1, 0}};
  int[][] c0111 = {{0, 1}, {1, 1}};
  int[][] c1000 = {{1, 0}, {0, 0}};
  int[][] c1001 = {{1, 0}, {0, 1}};
  int[][] c1010 = {{1, 0}, {1, 0}};
  int[][] c1011 = {{1, 0}, {1, 1}};
  int[][] c1100 = {{1, 1}, {0, 0}};
  int[][] c1101 = {{1, 1}, {0, 1}};
  int[][] c1110 = {{1, 1}, {1, 0}};
  int[][] c1111 = {{1, 1}, {1, 1}};

  for(int j = n-2; j > -1; j --){
  for(int i = 0; i < n-1; i ++){
    int[][] c = {{cell[i][j], cell[i+1][j]}, {cell[i][j+1], cell[i+1][j+1]}};
    
    // c の種類に応じて cn を変える
    if(check(c, c0001) == true){cn = c0001;}
    if(check(c, c0010) == true){cn = c0010;}
    if(check(c, c0011) == true){cn = c0011;}
    if(check(c, c0100) == true){cn = c0001;}
    if(check(c, c0101) == true){cn = c0011;}
    if(check(c, c0110) == true){cn = c0011;}
    if(check(c, c0111) == true){cn = c0111;}
    if(check(c, c1000) == true){cn = c0010;}
    if(check(c, c1001) == true){cn = c0011;}
    if(check(c, c1010) == true){cn = c0011;}
    if(check(c, c1011) == true){cn = c1011;}
    if(check(c, c1100) == true){cn = c0011;}
    if(check(c, c1101) == true){cn = c0111;}
    if(check(c, c1110) == true){cn = c1011;}
    if(check(c, c1111) == true){cn = c1111;}
    if(check(c, c0000) == true){cn = c0000;}

    // cn を cell にあてはめる
    cell[i][j] = cn[0][0];
    cell[i+1][j] = cn[0][1];
    cell[i][j+1] = cn[1][0];
    cell[i+1][j+1] = cn[1][1];
  }
  }
  
  return(cell);
}

boolean check(int[][] a, int[][] b){ // 配列 a と b を比較
  if(a[0][0] == b[0][0] &&
     a[0][1] == b[0][1] &&
     a[1][0] == b[1][0] &&
     a[1][1] == b[1][1]){
    return(true);
  }
  else{
    return(false);
  }
}


先に16種類の 2 × 2 のセル c0000 ~ c1111 をつくっておき、今注目している部分はどれに当てはまるかを check() で調べます。
それに応じて変化させ、その結果を cn に入れます。
最後に cn を元の cell に戻して完了、という感じです。


やり方を見直す

素直に実装したりちょっと工夫してみたりしましたが、どちらにせよコードがやたら長くなってしまって分かりづらいです。
また、この後の話になりますが、「粉」の種類を増やす点でもやりづらい部分があります。


ということで、やり方を見直さなければなりません。

初めの16種類の変化でやりたいことは、「粉」が下に落ちることと、「粉」が重なっていたら上のものが左右どちらかに落ちること、と考えられます。
というわけで、これを実装すると、

int[][] calc(int[][] cell){
  for(int j = n-2; j > -1; j --){
  for(int i = 0; i < n-1; i ++){
    if(cell[i][j] == 1){
      if(cell[i][j+1] == 0){ // 下に何もなかったら
        cell[i][j+1] = 1; // 下に動かす
        cell[i][j] = 0; // 元々あったところから消す
      }
      else if(cell[i][j+1] == 1){ // 下に「粉」があったら
        int a = 2;
        int p = int(random(-a, a)); // 移動量を決める
        if(i+p < 0 || n-1 < i+p){ // 画面外に出そうだったら
          continue; // 何もしない
        }
        else if(cell[i+p][j] == 0){ // 移動予定地に何もなかったら
          cell[i+p][j] = 1; // 動かす
          cell[i][j] = 0; // 元の位置のは消す
        }
      }
    }
  }
  }
  
  return(cell);
}

やり方の性質上、「粉」が常に動くような感じにはなってしまいますが、一応それなりに単純に書けるようになりました。


下に「粉」があったときの処理を、

else if(cell[i][j+1] == 1 || cell[i][j+1] == 2){
  float p = random(0, 100);
  if(p < 1){
    if(i-1 < 0){
      continue;
    }
    else if(cell[i-1][j] == 0){
      cell[i-1][j] = 1;
      cell[i][j] = 0;
    }
  }
  else if(p >= 99){
    if(n-1 < i+1){
      continue;
    }
    else if(cell[i+1][j] == 0){
      cell[i+1][j] = 1;
      cell[i][j] = 0;
    }
  }
}

のようにすると、p の大きさや移動方向の分岐をどの程度の p までにするかを調整することで、左右の動きが小さい「粉」をつくることもできるかと思います。


「粉」の種類を増やす

これまでで「粉遊び」的な基本となる部分はできたかと思います。
ここからはおまけ的な感じで、「粉」の種類を増やしてみます。


セル・オートマトンの良いところですが、セルの情報は割り当てる数によって好きなように設定できます。
ここでは、1 を水、2 を油、3 を壁に対応させてみます。

int[][] calc(int[][] cell){
  for(int j = n-2; j > -1; j --){
  for(int i = 0; i < n-1; i ++){
    if(cell[i][j] == 1){ // 水に関する処理
      if(cell[i][j+1] == 0){
        cell[i][j+1] = 1;
        cell[i][j] = 0;
      }
      else if(cell[i][j+1] == 1 || cell[i][j+1] == 3){
        int a = 2;
        int p = int(random(-a, a));
        if(i+p < 0 || n-1 < i+p){
          continue;          
        }
        else if(cell[i+p][j] == 0 || cell[i+p][j] == 2){
          cell[i][j] = cell[i+p][j];
          cell[i+p][j] = 1;
        }
      }
      else if(cell[i][j+1] == 2){
        float p = random(0, 1);
        if(p < 0.5){
          cell[i][j+1] = 1;
          cell[i][j] = 2;
        }
      }
    }
    if(cell[i][j] == 2){ // 油に関する処理
      if(cell[i][j+1] == 0){
        cell[i][j+1] = 2;
        cell[i][j] = 0;
      }
      else if(cell[i][j+1] == 1 || cell[i][j+1] == 2 || cell[i][j+1] == 3){
        int a = 2;
        int p = int(random(-a, a));
        if(i+p < 0 || n-1 < i+p){
          continue;          
        }
        else if(cell[i+p][j] == 0){
          cell[i+p][j] = 2;
          cell[i][j] = 0;
        }
      }
    }
  }
  }

  return(cell);
}

突貫工事感もありますが、油 (灰色) は水 (白色) に浮く一方で、水は油の下に沈んでいくという表現ができたのではないでしょうか。
なんとなくですが、水が沈んでいくときに抵抗があるような感じにしてみました。


これを初めの方法でやるとなると、{{1, 0}, {2, 2}}{{2, 0}, {1, 2}} に変える、などとやっていけばいいのかもしれませんが、面倒です。


まとめ

以上、processing で「粉遊び」的なものをつくってみました。

皆さんも、自分なりのルールで、多種多様な「粉」をつくって遊んでみてください。


最後まで読んでいただいてありがとうございました。また次回。

processingの備忘録 -音声-

こんにちは。
今回は「processing で音声を扱ってみよう」です。


かなり間が空きましたが、過去に画像とテキストの生成・読み込みを取り上げています。

processingの備忘録 -png、gifの生成- - プログラミングの備忘録
processingの備忘録 -txt、csvの生成- - プログラミングの備忘録
processingの備忘録 -画像、テキストの読み込み- - プログラミングの備忘録

今回は「ファイルの扱いシリーズ」として、さっくりとですが音声ファイルに手を伸ばしてみます。


目次


ライブラリ

processing における音声の取り扱いには、主に2つのライブラリがあります。
「Sound」(Sound / Libraries / Processing.org) と
「Minim」(https://code.compartmental.net/minim/Overview) です。

せっかくなので、この記事ではどちらも取り上げます。


インポートする際には、

  • Sound
import processing.sound.*;


  • Minim
import ddf.minim.*;


また、言わずもがなですが、なんらかの音声ファイルもしくはマイクが必要になります。
音声ファイルはスケッチフォルダに入れておいてください。


音声ファイルの読み込みと再生

このコードでは「test.mp3」という音声ファイルをスケッチフォルダに入れておく必要があります。


  • Sound
import processing.sound.*;

SoundFile file;

void setup() {
  file = new SoundFile(this, "test.mp3"); // ファイルの読み込み
  file.play(); // 再生
}      

void draw() {
}


  • Minim
import ddf.minim.*;

Minim minim;
AudioPlayer player;

void setup(){
  minim = new Minim(this); // 初期化
  player = minim.loadFile("test.mp3"); // ファイルの読み込み
  player.play(); // 再生
}

void draw(){
}


マイクからの入力

  • Sound

マイクから受け取った音声の音量を、円の大きさで表します。

import processing.sound.*;
AudioIn in;
Amplitude amp;

void setup() {
  size(500, 500);
  stroke(255);

  in = new AudioIn(this); // マイクの音声を受け取る
  //in.play(); // 受け取った音声を再生

  amp = new Amplitude(this);
  amp.input(in); // 受け取った音声の音量
}      

void draw(){
  background(0);
  
  float d = 300 * amp.analyze(); // 音量を円の直径に
  ellipse(width/2, height/2, d, d);
}


  • Minim

マイクから受け取った音声を、波形として画面に図示します。

import ddf.minim.*;
 
Minim minim;
AudioInput in;

void setup(){
  size(500, 500);
  stroke(255);

  minim = new Minim(this);
  in = minim.getLineIn(); // マイクの音声を受け取る
  //in.enableMonitoring(); // 受け取った音声を再生
}

void draw(){
  background(0);

  for(int i = 0; i < in.bufferSize()-1; i++){
    line(i, 100 + 100*in.left.get(i), i+1, 100 + 100*in.left.get(i+1)); // 左マイクの音
    line(i, 400 + 100*in.right.get(i), i+1, 400 + 100*in.right.get(i+1)); // 右マイクの音
  }
}


音の再生

これまでは音声ファイルやマイクからの入力音声を再生していましたが、波形から音を再生することもできます。


  • Sound
import processing.sound.*;

SinOsc sine;

void setup() {
  sine = new SinOsc(this); // サイン波の生成
  sine.freq(500); // 周波数を 500 Hz に
  sine.play(); // 再生
}

void draw() {
}


  • Minm
import ddf.minim.*;
import ddf.minim.signals.*; //波をつくるためのライブラリ

Minim minim;
AudioOutput out;
SineWave sine;

void setup(){
  size(500, 500);
  minim = new Minim(this);

  out = minim.getLineOut();

  sine = new SineWave(500, 0.5, out.sampleRate()); // サイン波の生成
  
  out.addSignal(sine); // 出力
}

void draw(){
  background(0);
  stroke(255);

  for(int i = 0; i < out.bufferSize()-1; i++){
    line(i, 100 + 100*out.left.get(i), i+1, 100 + 100*out.left.get(i+1));
    line(i, 400 + 100*out.left.get(i), i+1, 400 + 100*out.left.get(i+1));
  }
}


サイン波の他、矩形波三角波、のこぎり波などもあります。

また合成波に関しては、sine1sine2 など複数の波を生成させ、それぞれについて play()addSignal() を行えば良いということになります。


おまけ

少し前にフーリエ変換について記事にしました。

processingの備忘録 -離散フーリエ変換- - プログラミングの備忘録
processingの備忘録 -高速フーリエ変換- - プログラミングの備忘録
フーリエ変換を利用して絵を描く - プログラミングの備忘録


フーリエ変換ができるようになったら、音を変換して周波数成分を見てみたくなりますよね?


離散フーリエ変換できたてほやほやの頃にやってみたりしてました。
(ちなみに、音声としては木村わいPの「高音厨音域テスト」を使ってみました。今年8月には10周年バージョンが出てましたね。)

(ライブラリ使うんだったら組み込まれてる FFT() 使えば良くない?というコメントは受け付けていません。)


変換用のシグナルとしてマイクからの入力音声を使えば自分の声などもフーリエ変換できます。

今回はやりませんでしたが、フーリエ変換後のスペクトルからどこかの周波数成分を除いて波形に戻す、という操作で音がどう変わるかなど見てみるのも面白いかと思います。


まとめ

以上、processing で音声を扱えるようになりました。
全体として「Sound」の方が簡単な音声の取り扱いができるような気がします。
が、調べたいと思ったときに「sound」という語が一般的すぎてうまく検索できないなと感じました。

ここに挙げた機能が全てはないので、レファレンスを読みながらいろいろ開拓してみたいですね。
音声ファイルの書き出しについても機会があれば記事にしたいと思います。


最後まで読んでいただいてありがとうございました。また次回。

プログラミング言語「Piet」を触ってみる

こんにちは。
今回は、プログラミング言語「Piet」を紹介します。


実は最近、こんなツイートをしました。


これを見て"自己紹介"を感じていただけたならそれはそれで芸術の面白さを感じますが、芸術に疎い私にはこの画像に自己を乗せるというような高尚な表現はできないので、別の意味を持たせてあります。

それが今回紹介するプログラミング言語「Piet」に関わってきます。


目次


「Piet」とは

Piet は難解プログラミング言語のひとつです。
難解とよばれる所以は、この言語の特徴である画像ソースコードとする点にあります。


公式サイト
DM's Esoteric Programming Languages - Piet

wiki
Piet - Wikipedia


名前は画家のピエト・モンドリアンに由来し、「Mondrian」という名前にしようと思ったが既にあったので「Piet」にした、といういきさつがあるようです。

ピエト・モンドリアンは『赤・青・黄のコンポジション』『ブロードウェイ・ブギウギ』などの抽象画で有名ですが、これをソースコードにしようという発想が面白い。
実際に書いてみると、はじめに挙げたようにモザイクに近いような見た目になることが多いです。

寄せて作ってみるとこんな感じになります。


仕様

Piet は画像をソースコードとしています。
実行時は基本、左上から時計回りにコーデルをたどっていき、移動の前後での色相 (Hue) と明度 (Lightness) の変化によって実行する命令を指定できます。

コーデルは色を置いてあるマスのことを指し、同じ色が接している部分ではひとまとまりとみなされます。

数値は同じ色のコーデルが何個接しているかによって表されます。


コーデルを移動する"インタプリタ"は、DP (Direction Pointer) とCC (Codel Chooser) を持っています。
DPとCCはそれぞれ、車と乗っている人それぞれの向く方向と考えると分かりやすいです。

基本的には"壁"にぶつかるまでDP方向に直進し、左上から時計回りにぐるぐる回っていくと考えてよいので、詳しい説明は省きます。


また、記憶領域としてスタックを1つ持ち、後入れ後出しで値を取り扱うことができます。


使うことのできる色は、基本色のRYGCBMとその明暗で18色、そして白黒の2色で計20色です。


既に、移動前後の色相と明度の変化により実行される命令が決まる、と書きましたが、それは以下のように決まっています。

命令はそれぞれ、

  • push: スタックに値を1つ入れる

  • pop: スタックから値を1つ出す

  • addsubtractmultiplydivide: スタックから値を2つ出して加減乗除を計算し、結果をスタックに入れる
    (除は「2つ目の値 / 1つ目の値」)

  • mod: スタックから値を2つ出し、2つ目の値を1つ目の値で割った余りをスタックに入れる

  • not: スタックから値を1つ出し、その値が0なら1を、0以外なら0をスタックに入れる

  • greater: スタックから値を2つ出し、2つ目の値が1つ目の値よりも大きければ1を、そうでなければ0をスタックにいれる

  • pointer: スタックから値を1つ出し、その値の回数だけDPを時計回りに回転させる
    (負の数なら反時計回り)

  • switch: スタックから値を1つ出し、その値の回数だけCCを切り替える
    (負の数でも同じ)

  • duplicate: スタックから値を1つ出し、同じ値を2つスタックに入れる

  • roll: スタックから値を2つ出し、1つ目の値の回数だけ、2つ目の値だけの深さへの"roll"操作を行う
    "roll"操作は、一番新しい値を"2つ目の値"番目に入れ、それよりも上にある値を1段ずつ引き上げるという操作
    (回数が負の数なら逆の操作になる)

  • in: 標準入力から値を1つ受け取り、スタックに入れる
    (numberでは数字を、charでは文字を10進数の unicode で受け取る)

  • out: スタックから値を1つ出し、標準出力に出力する
    (numberでは数字を、charでは数字を文字に変換 (10進数の unicode) して出力する)

を意味しています。

また、白では何もせず、黒は"壁"としての役割を持ちます。


基本的な部分は説明できたかと思いますが、不十分な部分もあると思いますので適宜他のサイトを見てください。

[Piet入門編]Piet入門1 - 趣味プログラミングの部屋


コーディング

さて、だらだらと説明を書きましたが、百聞は一見に如かずとも言いますし、やってみた方がわかりやすいかと思います。


実際にコードを"描いて"いくわけですが、ペイントソフトや Excel でやるとなると、1ドット単位での操作や命令の難解さで厳しいです。

エディタが提供されていますので、積極的に利用しましょう。

piet-editor.github.io

github.com


今回は Pietron を使います。

実行画面はこんな感じです。


命令は色相と明度の変化で決まるので、初めの色は何でも大丈夫です。
なので、以下に挙げた12種類のコードは、全て明度変化1 (1 Darker) で、1を push する命令となっています。


自己紹介

それでは、この記事の初めに挙げた「自己紹介」のコードを作ってみます。

自己紹介というタイトルは、その名の通り、名前を出力するコードであるということでつけました。
すなわち、実行することで「taq」と出力されます。


Piet で文字を出力するためには out(c) を使うわけですが、これはスタックにある数字を unicode とみなして文字に変換して出力しています。
なので、「t」「a」「q」それぞれを表す unicode を調べ、それを四則演算等を用いて表現する必要があります。

調べると「t」「a」「q」はそれぞれ 116、97、113 で表されることがわかりました。


続いて、これらをどのように計算するかを考えます。
「t」の 116 を例にしてみます。

$116 = 4 \times (4 \times 7 + 1)$ と表されるので、

4 (push 4)
4 4 (dup)
4 4 4 (dup)
4 4 4 4 (dup)
4 4 4 4 1 (push 1)
4 4 4 3 (sub)
4 4 7 (add)
4 28 (multi)
4 28 1 (push 1)
4 29 (add)
116 (multi)

という風にスタックを操作していって 116 を作りました。
最後に out(c) を実行して「t」を出力します。

10 × 10 のコーデルで「t」を出力するコードは以下のようになります。
念のため"壁"である黒を隅に置きましたが、無くして詰めても問題ないです。

同じ要領で 97 と 113 も作っていくと、

ぴったり収まりました。偶然です。


Pietron にはデバッグモードがあり、F7 を押すことで1つずつ進む様子を見ることができます。

スタックが変化していって、最終的に「taq」が出力される様子が見えて面白いです。


これで完成でもよかったのですが、一本道が明らかに見えていたのでなんとなく周りに関係ない色を置いてカモフラージュしてみました。


書いていませんでしたが、Piet の処理は同じ白道を何回かたどると終了するので、右端はあえて白にしています。


本家に寄せたバージョン

この言語はピエト・モンドリアンに着想を得ているということで、『赤・青・黄のコンポジション』に似せた感じで作ってみました。

経路を矢印で示してみると、

数値に関しては、四角形にぴったり収まらなかったので似た色で調整しています。


まとめ

以上、難解プログラミング言語「Piet」を紹介しました。

私自身 processing を扱っていながら芸術方面にはめっぽう疎いですが、ただの絵だと思っていたら実はソースコードで実行すると隠されたメッセージが読める、みたいなものは現代アートという感じがしてとても好きです。

また、ネタ言語とはいえチューリング完全で、今回は roll やDP・CCの操作を使っていないので、できることはもっとあろうかと思います。

みなさんもぜひ Piet を使って、processing とはまた違った「プログラミング×芸術」を楽しんでみてください。


最後まで読んでいただいてありがとうございました。また次回。

「Stone Story RPG」で遊んでみる

こんにちは。
今回は「Stone Story RPG」というゲームを遊んでみます。


ちょっと長めのイントロが入ります。

友人と Steam の無料ゲームを発掘することがままあるのですが、そこで「Super Auto Pets」というゲームを見つけました。
これがなかなかおもしろく、友人が見つけてくれたんですが、以前にも「Muck」という良ゲーを見つけてくれたので見る目あるなと思っていました。

「Super Auto Pets」ですが、アプリ版もあるようだということでせっかくなので入れることにしました。
アプリを入れるときは、おすすめのゲームみたいなものがいくつか表示されると思うのですが、その中に自分の好きな雰囲気のゲームを見つけ、それが「Stone Story RPG」でした。
実際に遊んでみても自分の好きなタイプのゲームで、これまたいいもの見つけちゃったなと思っていたところ、自分でコードを書いて周回させるという要素が解放されました。
思わぬところでプログラミングに出会ってしまったのでこれは記事にしておこうということで、この記事を書くに至りました。


思い返してみると、確か「INDIE Live Expo」という配信で紹介されていて、そこで周回のためのコードを書くというシステムがおもしろいなと気になっていたものかもしれないと、謎の運命を感じてしまいました。
(全然関係ないですが、この配信に出演されているわいわいさんは私の推しYouTuberなのでぜひ見てみてください。)


Stone Story RPG 自体は Steam にもあり、3400円で販売されています。アプリ版は無料です。


一応「プログラミングの備忘録」という名前なので、コードが関わる部分だけ紹介するに留めます。
加えてネタバレを含むと思われますので、問題ない方のみ読み進めていただければと思います。


目次


Stone Story RPG におけるプログラミング

大前提として、「神殿」のボス「ナーガラージャ」を倒して「明瞭の石」を手に入れる必要があります。
実際のコーディング画面はこんな感じです。

言語は「Stonescript」という、おそらく独自の言語です。
リファレンス等はこちらのサイトにあります。

stonestoryrpg.com


デフォルトでは以下のように書かれています。

?hp < 7
  activate potion
?loc = caves
  equipL sword
  equipR shield
  ?foe = boss
    equip crossbow


「HPが7より少なくなったらポーションを使う」と「ステージが恐怖の洞窟のとき、左手に剣、右手に盾を装備し、ボス戦ではクロスボウを装備する」という内容が書かれています。

この Stonescript という言語は、if文を「?」で書くことに特徴があると思います。
他の言語でも三項演算子で「?」を使いますが、それとはまた違った文法です。


実際にやってみる

レファレンスは先程示したサイトにあるので、軽くですが書きながら勉強してみました。


例えば、普段は氷の大剣と丈夫な盾で戦うけど、体力が減ったら活力の大剣で回復する、というようにすると、

equipL big ice sword // 氷の大剣を装備
equipR compound // 丈夫な盾を装備

? hp < maxhp // 体力が1でも減ったら
  equipL big vigor sword // 活力の大剣に切り替える


ボス戦ではヘビーハンマーを使い、適宜アビリティも発動させる、なら、

? foe = boss // ボス戦のとき
  equip hammer // ヘビーハンマーを装備
  ? item.CanActivate() & // アビリティが発動できて
  ^foe.distance < 5 // 敵との距離が5より小さいとき
    activate right // 発動


装備品の名前については日本語で書くこともできます。

また、1行に書ける文字数は限られているようなので、長くなりそうなときは ^ を使って前の行とつながっていることを示した上で改行して書くと良いです。

他、素材の近くになったら「星の石」で集める、ただの移動の時は「トリスケルの石」で加速する、攻撃されたときは「明瞭の石」で避けるなど、いろいろとできることがあると思います。


APDE や pydroid を使ったときも思いましたが、やはりスマホでコードを書くのは大変ですね。
なぜかアローキーで入力位置の移動もできなかったので、難易度が上がりました。


まとめ

以上、あっさりとですが Stone Story RPG における Stonescript を使った周回の自動化を試してみました。

ただの自動周回では装備が同じであったりポーションや特殊攻撃を使ってくれないなどの、自分で周回すればクリアできるのに…という悩みを解決してくれる素晴らしいシステムだと思います。

「Stone Story RPG」みなさんもぜひ遊んでみてください。


最後まで読んでいただいてありがとうございました。また次回。

フーリエ変換を利用して絵を描く

こんにちは。
今回は「フーリエ変換を利用して絵を描く」ということをやってみます。


フーリエ変換の利用法としてたまに挙げられるのが、「フーリエ変換の結果から周転円をつくって絵を描かせる」というものだと思います。

例えば、


これまでにフーリエ変換シリーズとして2記事書いてきました

taq.hatenadiary.jp

taq.hatenadiary.jp

せっかくフーリエ変換に触れたので、今回はこれらのまとめ的な感じでやってみます。


ちなみに、正式名称はあるのでしょうか。
フーリエ変換の結果から周転円をつくって絵を描かせる」というと長いです。

上記のサイトには「fourier epicycle drawing」という英語があったので、とりあえずそうよんでおくこととします。


目次


方法

そもそもフーリエ変換の根底には「フーリエ級数」という、複雑な周期信号は単純な波の和で表すことができる、という考え方があります。
そして、この "単純な波" がどういうものなのかを探すのがフーリエ変換になります。

単純な波は円を使って表現することができるので、周転円、すなわち波の足し合わせを使うことで元の周期信号を描くことができる、という訳です。


手順としては、

  1. 描きたい絵の座標を用意する

  2. フーリエ変換する (同時に周波数、振幅、位相も計算しておく)

  3. 結果を周転円に変換する

  4. 実行

です。


実装

基本的なつくりは我らが Shiffman 先生の動画を参考にさせて頂きました。

Coding Challenge 125: Fourier Series - YouTube


絵の座標を用意する

大前提として何かしら絵が必要です。
一筆書きかつ始点と終点が同じものだと周期的になるのでよいかと思います。

今回は計算する前にキャンバスに手で描いて入力することにしました。


int N;
float[] signalX = {};
float[] signalY = {};
Complex[] signal;
int points = 0;
boolean pause = true; // ポーズ中

void mouseDragged(){ // マウスをドラッグ中は
  if(pause == true){
    // マウスの軌跡を配列に入れていく
    signalX = splice(signalX, mouseX-width/2, points);
    signalY = splice(signalY, mouseY-height/2, points);
  
    points ++;
  }
}

void mouseReleased(){ // ドラッグを止めたら
  if(pause == true){
    N = points - 1;
    println(N);
    
    // FFTのアルゴリズム上、要素数は2の冪にしたい
    N = int(pow(2, ceil(log(N)/log(2))));
    println(N);
    
    signal = new Complex[N];
    for(int i = 0; i < points-1; i ++){
      signal[i] = new Complex(signalX[i], signalY[i]);
    }
    // 余りの部分は終点の座標を入れておく
    for(int i = points-1; i < N; i ++){
      signal[i] = new Complex(signalX[points-2], signalY[points-2]);
    }
    
    pause = false; // ポーズを止める
  }
}


フーリエ変換する

続いて、フーリエ変換を行います。

離散フーリエ変換 (DFT) でも問題ないとは思いますが、せっかくなので高速フーリエ変換 (FFT) を使って実装します。 FFT() に関しては、本ブログの高速フーリエ変換の記事で作成したものを使っています。

先程は「同時に周波数、振幅、位相も計算しておく」と書きました。
DFTの場合は DFT() 関数内に組み込むことができますが、FFTの場合はアルゴリズムの関係上 (多分) 組み込むことができないので、使う際に計算することにしています。

if(pause == false){
  fourier = new Complex[N];
    
  // フーリエ変換
  fourier = FFT(signal);
  
  for(int i = 0; i < N; i ++){
    float n = i; // 周波数
    float r = fourier[i].abs()/N; // 振幅
    float phi = atan2(fourier[i].im, fourier[i].re); // 位相
  }


結果を周転円に変換

フーリエ変換とその後の計算によって周転円の描画に必要な値が求まったので、これを元に描画します。

if(pause == false){
  // 周転円の円 (N = 0) の中心座標
  float x = width/2;
  float y = height/2;

  for(int i = 0; i < N; i ++){
    float prex = x;
    float prey = y;
    
    // 次の円の中心座標
    x += r * cos(n * theta + phi);
    y += r * sin(n * theta + phi);
    
    // 周転円の描画
    stroke(150);
    noFill();
    ellipse(prex, prey, 2*r, 2*r);
    line(prex, prey, x, y);
    fill(150);
    ellipse(x, y, 5, 5);
  }
}


実行

必要なものは揃ったので、後は変換前後の絵を表示するようにしたり、フーリエ変換の関数や複素数のクラスを組み込んで完成です。

float theta;
int N;
float[] signalX = {};
float[] signalY = {};
Complex[] signal;
Complex[] fourier;
ArrayList<PVector> trace = new ArrayList<PVector>();
int points = 0;
boolean pause = true;

void setup(){
  size(500, 500);
}

void mouseDragged(){
  if(pause == true){
    signalX = splice(signalX, mouseX-width/2, points);
    signalY = splice(signalY, mouseY-height/2, points);
  
    points ++;
  }
}

void mouseReleased(){
  if(pause == true){
    N = points - 1;

    N = int(pow(2, ceil(log(N)/log(2))));
    
    signal = new Complex[N];
    for(int i = 0; i < points-1; i ++){
      signal[i] = new Complex(signalX[i], signalY[i]);
    }
    for(int i = points-1; i < N; i ++){
      signal[i] = new Complex(signalX[points-2], signalY[points-2]);
    }
    
    pause = false;
  }
}

void draw(){
  background(0);
  
  for(int i = 0; i < points-1; i ++){
    stroke(255, 0, 0);
    line(signalX[i]+width/2, signalY[i]+height/2, signalX[i+1]+width/2, signalY[i+1]+height/2);
  }
  
  if(pause == false){
    fourier = new Complex[N];
    
    float x = width/2;
    float y = height/2;
    
    fourier = FFT(signal);
    
    for(int i = 0; i < N; i ++){
      float n = i;
      float r = fourier[i].abs()/N;
      float phi = atan2(fourier[i].im, fourier[i].re);
      
      float prex = x;
      float prey = y;
      
      x += r * cos(n * theta + phi);
      y += r * sin(n * theta + phi);
      
      stroke(150);
      noFill();
      ellipse(prex, prey, 2*r, 2*r);
      line(prex, prey, x, y);
      fill(150);
      ellipse(x, y, 5, 5);
    }

    trace.add(new PVector(x, y));

    translate(0, 0);
    stroke(255);
    noFill();
    beginShape();
    for(PVector t: trace){
        vertex(t.x, t.y);
    }
    endShape();
    
    theta += 2*PI/N;
  }
}

Complex[] FFT(Complex[] input){
  int N = input.length;
  Complex[] output = new Complex[N];
  
  if(N == 1){
    return input;
  }
  else{
    Complex[] even = new Complex[N/2];
    Complex[] odd = new Complex[N/2];
    
    for(int i = 0; i < N/2; i ++){
      even[i] = input[i].add(input[N/2+i]);
      
      Complex w = new Complex(cos(2*PI*i/N), -sin(2*PI*i/N));
      odd[i] = (input[i].sub(input[N/2+i])).mult(w);
    }
    
    Complex[] array_even = new Complex[N/2];
    Complex[] array_odd = new Complex[N/2];
    
    array_even = FFT(even);
    array_odd = FFT(odd);
    
    for(int i = 0; i < N/2; i ++){
      output[2*i] = array_even[i];
      output[2*i+1] = array_odd[i];
    }
   
    return output;
  }
}

class Complex{
  float re, im;
  
  Complex(float _re, float _im){
    re = _re;
    im = _im;
  }
  
  Complex add(Complex a){
    return new Complex(re + a.re, 
                       im + a.im);
  }
  
  Complex sub(Complex a){
    return new Complex(re - a.re, 
                       im - a.im);
  }
  
  Complex mult(Complex a){
    return new Complex(re * a.re - im * a.im, 
                       re * a.im + im * a.re);
  }
  
  float abs(){
    return sqrt(re*re + im*im);
  }
}


実行してみると、

初めの手描きの絵 (赤線) がフーリエ変換され、周転円による描画 (白線) が行われます。
そして、描かれた絵は元の絵と見事に重なっていると思います。


余談

絵について

今回は手書きの絵を使いましたが、座標から指定する方法もあります。
ただしその場合は、少なくとも一筆書きであるという条件は満たしておかないと描画がおかしくなってしまいます。

点の集まりから一筆書きのルートを探す、雑誌やなんかによくあった順番通りになぞると絵が完成するやつみたいなことをしなければなりません。

探してみるとあるもので、いわゆる「巡回セールスマン問題」や「中国人郵便配達問題」とよばれるようなものを考えることで (アニメキャラである必要はありませんが) 点や辺の集まりから一筆書きの最短経路を求めることができます。

任意のアニメキャラを一筆書きで数式にするやつ (Pythonで) - メモ帳の裏の切れ端


ちなみに、アイキャッチの「FOURIER」は自分で座標を計算して一点一点配列に打ち込んで作っています。つまりゴリ押し。


計算・描画方法について

今回は $x$、$y$ を1つの複素数の実部、虚部それぞれにあてはめてフーリエ変換することで、周転円1つで fourier epicycle drawing を行いました。

この $x$、$y$ をそれぞれフーリエ変換して周転円2つに割り当てれば、はじめに紹介した動画のようなものができます。

そして、このような計算を2つの絵に対して同時に行う、例えば絵1の $x$、$y$ を周転円1の実部、周転円2の虚部それぞれに、絵2の $x$、$y$ を周転円1の虚部、周転円2の実部それぞれにあてはめる、というように計算をすると、少し前に話題になったツイート (エックセズ?ポスト?) のようなものができあがります。

(もし本当にあの青い鳥と $\mathbb{X}$ マークの間に隠れた関係があるのなら激アツ展開ですが、実際はどんな絵でも同じようなことができます。)


まとめ

以上、フーリエ変換を利用して任意の絵を描くような動きをさせるという「fourier epicycle drawing」を processing に実装することができました。

つくりたてほやほやのような感じで記事を書いているので、終点に達したときに不自然に間が空いたり、元の絵の始点と終点が同じ所にないと変な線が入ったりと改善点はあります。
その辺りは今後ブラッシュアップしていければと思います。


最後まで読んでいただいてありがとうございました。また次回。