メトロポリス輸送の事始め

2014-12-21
#processing #math

メトロポリス輸送の事始め

前から双方向パストレーシングを理解したいと思ってたんだけど、視点からのパスとライトからの パスの結びつけるときのウェイトの与え方がわからなくて、なんかいい資料がないかとさまよっていた ところ、双方向パストレーシングとは違うんだけどメトロポリス光輸送の論文が引っかかった。

A Practical Introduction to Metropolis Light Transport

手を動かしながら読み進めていきたい。

１章は導入だけなので、２章から。

2 Metropolis Transport

２章のタイトルは「メトロポリス輸送」とのことで、まず何をやりたいのかを語る：

  • 関数 f を近似したいとして、
  • f に比例するサンプリング分布を生成し、その分布から 得られるサンプルのヒストグラムを作る
  • 結果のヒストグラムは f に比例する（当然）
  • f を近似するようにヒストグラムをスケール

f に比例するサンプリング分布が作れるという前提からして、ちょっと何言ってるのかわからない ですなんだけど、メトロポリス輸送 はこの方法を使って関数を近似する、ということらしい。

2.1 Creating a Sampling Distribution

  • メトロポリスアルゴリズムは詳細釣り合い というアイディアで、f に比例する分布を作り出す

で、その後詳細釣り合いについてあれこれ書かれているんだけど、なにを書いているのか さっぱりわからず…。

  • すでに f に比例するヒストグラムが存在するとして
  • あるビンから他のビンに流れ出す量が戻る量と同じなら 定常分布 が保たれる （これを 詳細釣り合い という）

なんちゃらかんちゃら…

遷移関数を定義する

  • 関数 K は 暫定遷移関数 T （または 変異戦略 と呼ばれる）を使って定義される
  • は、現在のサンプル位置 から次のサンプル位置の提案として が選ばれる 確率
  • 確率 で から へ遷移し、 で に留まる
  • 図2のmakeHistogram関数はメトロポリス輸送を使って画像をコピーする
  • すごく単純な変異戦略で、一様確率で画像のランダム位置を選ぶ
  • しかし幅広い遷移関数が使える、後でMLTの力はいい変異戦略を選ぶことに依ることを見る

でCのプログラムが載っている。こいつをProcessingで動かせるようにするために、 Javaに書き直してみる：

void makeHistogram(int w, int h, float[] F, float[] histogram, int mutations) {
  // Create an initial sample point
  int x0 = randomInteger(0, w - 1);
  int x1 = randomInteger(0, h - 1);
  float Fx = F[x1 * w + x0];

  // In this example, the tentative trasition function T simply chooses
  // a random pixel location, so Txy and Tyx are always equal.
  float Txy = 1.0 / (w * h);
  float Tyx = 1.0 / (w * h);

  // Create a histogram of values using Metropolis sampling.
  for (int i = 0; i < mutations; ++i) {
    // choose a tentative next sample according to T.
    int y0 = randomInteger(0, w - 1);
    int y1 = randomInteger(0, h - 1);
    float Fy = F[y1 * w + y0];
    float Axy = min(1.0, (Fy * Txy) / (Fx * Tyx));  // equation 2.
    if (randomReal(0.0, 1.0) < Axy) {
      x0 = y0;
      x1 = y1;
      Fx = Fy;
    }
    histogram[x1 * w + x0] += 1;
  }
}

これに画像を読み込んで各ピクセルの明るさから目的の関数Fの配列を生成して、 関数makeHistogramを呼び出して得られたhistogramを描画するプログラムを書いてやる：

int randomInteger(int min, int max) {
  return floor(random(max - min + 1)) + min;
}

float randomReal(float min, float max) {
  return random(max - min) + min;
}

/* @pjs preload="mandrill.jpg"; */
PImage img;
float[] F;
float fAve;
float[] histogram;

void setup() {
  img = loadImage("mandrill.jpg");
  //size(img.width * 2, img.height);  // Javascriptでエクスポートした時にうまく動かない
  size(512, 256);

  F = new float[img.width * img.height];
  float fTotal = 0;
  for (int i = 0; i < img.height; ++i) {
    for (int j = 0; j < img.width; ++j) {
      color c = img.pixels[i * img.width + j];
      float b = brightness(c);
      F[i * img.width + j] = b;
      fTotal += b;
    }
  }
  fAve = fTotal / (img.width * img.height);
  drawHistogram(img.width, 0, img.width, img.height, F, 1);

  histogram = new float[img.width * img.height];
}

int count;
void draw() {
  makeHistogram(img.width, img.height, F, histogram, img.width * img.height);
  ++count;

  float hAve = count;  // because one makeHistogram call emits samples same as image size.
  float s = fAve / hAve;
  drawHistogram(0, 0, img.width, img.height, histogram, s);

  textSize(12);
  textAlign(LEFT, TOP);
  fill(255, 0, 0);
  text("" + count, 0, 0);
}

void drawHistogram(int x0, int y0, int w, int h, float[] histogram, float s) {
  for (int i = 0; i < h; ++i) {
    for (int j = 0; j < w; ++j) {
      float v = histogram[i * w + j] * s;
      set(j + x0, i + y0, color((int) v));
    }
  }
}

デモ

僕の頭では説明を読んだだけではわからなかったけど、プログラムを動かしてみるとなんとなく動作が 理解できた。これによると、

  • 暗い点から同じ明るさまたは明るい点には100%遷移する（Axyが1になるので）
  • 逆に明るい点から暗い点には、その明るさの比によって確率で遷移するかしないか決まる

この規則で、なんで f に比例する分布が得られるのかさっぱりわからないんだけど、 動かした結果は確かにそれっぽい。不思議。

関連記事

新着記事