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

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

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

例えば、

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

taq.hatenadiary.jp

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

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

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

方法
実装
余談
- 絵について
- 計算・描画方法について
まとめ

方法

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

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

手順としては、

描きたい絵の座標を用意する
フーリエ変換する (同時に周波数、振幅、位相も計算しておく)
結果を周転円に変換する
実行

です。

実装

基本的なつくりは我らが 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の実部それぞれにあてはめる、というように計算をすると、少し前に話題になったツイート (エックセズ？ポスト？) のようなものができあがります。