2008-07-12

回転行列の固有値

数学

回転行列の固有値を考えましょう。二次元の回転では、回転により方向を変えないベクトルは無いので、結果として固有値は複素数です。これはよく知られていますね。ちなみに三次元では、回転軸は回転で方向（と大きさ）を変えないので、固有値のひとつは1であることは自明。他の二つはたぶん複素数ですね。

上記の説明は直感的に理解できますが、本当かどうか気になったので試しに手計算で確認してみました。愚直に特性方程式から求めてますが、もっとスマートな方法があるような気がします。

2次元の場合

回転行列Rは
$R = \begin{pmatrix} cos\theta & -sin\theta \\ sin\theta & cos\theta \end{pmatrix}$
ですので、その特性方程式は
$|R - \lambda E| = \lambda^2 - 2cos\theta\lambda + 1 = 0$
となります。この時、次の判別式Dを使ってこの方程式の実数解が存在するかどうかをチェックすると、
$D = 4(cos^2\theta - 1)$
0<θ<2π の範囲において、-1 < cos^2θ-1 < 0 ですので、D < 0 です。したがって、特性方程式の実数解は存在せず固有ベクトルも存在しません。

3次元の場合

説明を簡単にするため、Ｚ軸周りの回転に限定します。回転行列Rは
$R = \begin{pmatrix} cos\theta & -sin\theta & 0 \\ sin\theta & cos\theta & 0 \\ 0 & 0 & 1 \end{pmatrix}$
ですので、その特性方程式は
$|R - \lambda E| = (1 - \lambda)(\lambda^2 - 2cos\theta\lambda + 1) = 0$
で実数解はλ=1だけです。回転はベクトルの長さを変えないので固有値が1となるのは当然ですね。では、λ=1に属する固有ベクトルを求めてみましょう。とりあえずλ=1を
$(R - \lambda E)\vec{x} = \vec{0}$
に代入すると、
$\begin{pmatrix} cos\theta - 1 & -sin\theta & 0 \\ sin\theta & cos\theta - 1 & 0 \\ 0 & 0 & 0 \end{pmatrix}\begin{pmatrix} x \\ y \\ z \end{pmatrix} = \vec{0}$
となります。この方程式を満たすベクトルを適当に選ぶと、(x, y, z) = (0, 0, k) が得られますが*1、これはＺ軸、すなわち回転軸と平行なベクトルとなります。結果として、３次元の回転行列の固有ベクトルはその回転軸と等しいことが分かりました。

*1:ここの導出がちょっと怪しい。

2008-06-23

シュリックの反射モデル (1)

CG プログラミング

自作レンダラーにシュリックの反射モデルを導入してみました。このモデルの発案者であるシュリック氏の論文「An Inexpensive BRDF Model for Physically-based Rendring」に従って素直に実装しただけですが、案の定、不自然な箇所が多くまだまだこれからです。
§
シュリックの反射モデルは、物体表面の荒さを調整するパラメータによって、拡散反射から鏡面反射までを連続して変化させることが可能なのが特徴です。滑らかな（荒さが小さい）物質ほど鏡面反射率が高くなるという性質は、とても直感的で理解しやすいですね。また、フォンの反射モデルと違い、拡散反射項と鏡面反射項を個別に設定する必要がなく扱いやすいのも利点です。ちなみに、原著論文では、拡散反射項d・光沢反射項g・鏡面反射項sを荒さrの値から自動的に決定する以下の式が提案されています。各項を手動であえて決定しないのは、おそらくパラメータ設定の単純さというこのモデルの利点を損なわないためでしょう。なので、もしも自分でモデルをカスタマイズするならば、独自の補間関数を使う方がベターかと思います。

if ( r < 0.5 ) { 
  g = 4*r*(1 - r);
  d = 0;
  s = 1 - g;
} else {
  g = 4*r*(1 - r);
  d = 1 - g;
  s = 0;
}

以下は、rを変化させて得られた画像で、左からr=0.1 r = 0.5 r = 0.9 となっています。鏡面反射（r = 0.1）のハイライトがもっと目立ってもいいような気がしますが、何か間違っているのでしょうか。今は r < 0.5 とすると、シーン全体が大変暗くなる傾向があります。

2008-06-11

ヒープ（バイナリヒープ）

プログラミング

Good Math, Bad Math : Binary Heaps
ヒープは、スタックやキューと並びシンプルで強力なデータ構造です。ヒープとは、木構造の一種で要素の挿入と最大値（又は最小値）の要素の削除の計算量が共にO(logN)。このヒープを2分木で表現したものを、バイナリヒープと呼びます。バイナリヒープの場合、要素の挿入と削除が行われると、常に根ノードと最大値の要素が常に等しくなるよう、ノードの入れ替えが行われます。
このようなヒープの性質は、優先順位付きキューなど、何らかの基準値によって要素が順位付けされた集合を表現するデータ構造に向いています。また、仕組みが単純な分、実装が容易な点も魅力的です。
恐らく多くの方々にとって、ヒープは一般的なアルゴリズムの教科書で説明されているように、整列プログラムを実装するためのデータ構造としてよく知られているのではないでしょうか。上のヒープの性質を応用すれば、簡単に配列の要素を整列出来ることがすぐに分かるでしょう。すなわち、配列から一つずつ要素を取り出してヒープに挿入すれば、後から自然と要素が整列された順番でヒープから要素が取り出すことが出来ます（なんて簡単！）。
ヒープソートの計算量は、ヒープの挿入と削除の計算量はO(logN)なので、Ｎ個の要素を整列する場合、O(NlogN + NlogN) = O(NlogN) です。代表的な整列アルゴリズムであるクイックソートの場合、平均計算量はO(NlogN)ですが、配列を分割するピボットの選び方によって最悪の場合O(N^2)かかる傾向があります。これに対して、ヒープソートは平均と最悪の場合で変わらず常にO(NlogN)しかかかりません。ただし、空間効率という点においてヒープソートは、ヒープを構築するために最低で元の配列の2倍のメモリサイズが必要なので、他の整列アルゴリズムに劣ります。しかし、PCのメモリ資源が潤沢な現在、計算コスト＞空間コストのような状況においては、ヒープソートの採用は十分に検討する価値があると思います。
以下は、バイナリヒープのサンプルプログラムです。こちらで紹介されているソースコードをRubyで再実装してみました。参考にどうぞ。

class Heap
  
  def initialize
    @list = []
  end

  def insert(value)
    @list << value
    if @list.length > 1
      last = @list.length - 1
      bubble_up(last)      
    end
  end
  
  def remove_largest
    root = @list.shift
    bubble_down(0)
    root
  end
  
  def to_s
    "[#{@list.join(',')}]"
  end  
  
  private
  
  def parent_of(pos)
    if pos == 0
      nil
    elsif pos == 1 || pos == 2
      0
    else
      if pos % 2 == 0
        (pos - 2)/2
      else
        (pos - 1)/2
      end  
    end
    
  end
  
  def children_of(pos)
    n = @list.length
    left = 2*pos + 1
    right = 2*pos + 2
    if n <= left
      []
    elsif n == right
      [left]
    else
      [left, right]      
    end
  end  

  def swap(i, j)
    tmp = @list[i]
    @list[i] = @list[j]
    @list[j] = tmp
  end
  
  def bubble_up(pos)
    unless pos == 0
      parent = parent_of(pos)
      if @list[parent] < @list[pos]
        swap(parent, pos)
        bubble_up(parent)
      end            
    end    
  end
  
  def bubble_down(pos)
    children = children_of(pos)
    unless children.empty?
      if children.length == 1
        i = children[0]
        if @list[pos] < @list[i]
          swap(pos, i)
          bubble_down(i)
        end
      else
        (i, j) = children
        left = @list[i]
        right = @list[j]
        root = @list[pos]
        if (left > root) || (right > root)
          if left > right
            swap(i, pos)
            bubble_down(i)
          else
            swap(j, pos)
            bubble_down(j)
          end
        end
      end
    end
  end
  
end

def heap_use_example
  heap = Heap.new

  heap.insert(1)
  heap.insert(3)
  heap.insert(5)
  heap.insert(7)
  heap.insert(9)
  p heap.to_s

  list = []
  list << heap.remove_largest
  p heap.to_s
  list << heap.remove_largest
  p heap.to_s
  list << heap.remove_largest
  p heap.to_s
  list << heap.remove_largest
  p heap.to_s
  list << heap.remove_largest
  p heap.to_s

  p list.join(',')  
end

heap_use_example

2008-05-26

微小立体角

CG 数学

これまでレンダリング方程式の中で表れる微小立体角（Δω）の導出がいまいち理解できてなかったのですが、たまたま本屋で「ゲームプログラミングのための3Dグラフィックス数学」を立ち読みしたところ、とても分かりやすかったのでメモしておきます。
§
天頂角をθ、方位角をφとした場合、3次元空間上の任意の点は球座標を使って以下のように表現できます。
$\begin{align} x &= r \cdot \sin\theta \, \cos\phi \\ y &= r \cdot \sin\theta \, \sin\phi \\ z &= r \cdot \cos\theta. \end{align}$
円弧の長さは半径と平面角（ラジアン）の積なので、方位角方向（XY軸）の弧の微小量は、 $r \cdot sin\theta\Delta\phi$ 　
同様に、天頂角方向（Z軸）の弧の微小量は、 $r \cdot \Delta\theta$ 　
これら天頂角方向と方位角方向の弧が成す微小平面の面積は各微小量の積と等しく、 $r^{2} \cdot sin\theta\Delta\theta\Delta\phi$ 　
ところで立体角とは、弧が成す球面上の面積を球の半径の2乗で割った値の大きさと定義されています（単位はステラジアン）。したがってその微小量、すなわち微小立体角は、
$\Delta\omega = sin\theta\Delta\theta\Delta\phi$
と求められます。
§
「ゲームプログラミングのための〜」は、スリーディーの加納さんも勧められてて是非手元に置いておきたい一冊ではあります。が、結構値が張るので、懐が気になる無職生活の身にはなかなか手が出しにくいところ。もう少し手頃な価格にならないかなあ。

2008-05-14

ノイズの削減と間接照明

CG プログラミング

約一ヶ月休んでましたが、色々蹴つまづきながらも開発は進んでます。いまは人生初めての就職活動に挑むにあたって、身辺がばたばたして心が落ち着きません。ああ、早く定職に就きたい．．．

ここ１週間は、レンダリング処理に様々な問題が見つかり、その解決に手間取りました。

ノイズ

サンプル数を増やしてもなかなかノイズが減少する気配がないので、疑問に思い詳しく調べたところ、シャドウレイによる光源の可視判定で予期しない失敗をしていることが分かりました。レイとポリゴンとの交差点から光源に向かってシャドウレイを追跡する時、始点を含むオブジェクトを誤って最近傍のオブジェクトとして判定してしまうパターンが多々発生し、これがオブジェクト表面に黒い斑点を生む原因となっていたようです。
そこで、レイの始点から交差点までの距離を互いに比較する際に、境界値を調整する補正値を加えたところ、期待通りにノイズが一気に減少しました。また、新しく層化サンプリングを実装したところ、今までオブジェクト表面に縞模様となって浮き上がっていたノイズも見事に消えて無くなりました。
以下は、修正後のプログラムでコーネルボックスをレンダリングした結果です。左が10パス／ピクセルで、右が100パス／ピクセルです。以前は100パス／ピクセルでも黒いノイズが派手に目立っていたのですが、修正後はそのようなノイズは見当たりません。

今回は「アルゴリズムは間違っていないのに、なぜ？」という思い込みが問題解決の邪魔をしてしまって、余計に時間がかかってしまいました。浮動小数点演算は慎重に行わないと怖いですね。以後の反省にします。

間接照明

これまで間接照明は、通常の光線追跡によって再帰的に求めた放射輝度にBRDFをかけて計算していました。しかし、BRDFの値が大変小さくなる傾向があり、間接照明の効果がほとんど得られていませんでした。実はこの間、試しに直接照明と間接照明を分けてレンダリングしてみたところ、間接照明だけの場合、シーン全体が真っ暗闇になってしまって「駄目だこりゃ」（いかりや長助調）となったしだい。
そこで、手元にある「 Realistic Ray Tracing 」（以下、RRT）の『12. Path Tracing』をじっくり読んでみると、放射輝度の計算式は、確率密度関数p(x)を
$p(\vec{k_i}) = \frac{\rho(\vec{k_i}, \vec{k_o})cos\theta}{R(\vec{k_o})}$
と定義すると、次の式に近似できるから単純に反射率Rをかけたらいいと書いてあります（最初に読んだ時は、意味が分からなくて完全にスルーしてました）。
$L_s(\vec{k_o}) = R(\vec{k_o})L_f(\vec{k_i})$
ここでいう反射率Rとは、あるオブジェクト表面において、全方向に出射される光量と一定方向から入射してくる光量の比であり、完全な拡散反射（ランバート面）だけに限定した場合、
$R(\vec{k_o}) = \pi \rho(\vec{k_i}, \vec{k_o})$
と定義されます。RRTでは、このRを「Directional Hemispherical Reflectance」と読んでいますが、正確な訳が分からないので、とりあえず反射率としておきます。
さて、これにしたがってBRDFに変えて反射率をかけるようプログラムを変更したところ、以前よりずっと柔らかい陰影をもつ画像が得られるようになりました（上記画像）。モンテカルロ積分において、確率密度関数がp(x)である関数F(x)の期待値は、
$g(x) = f(x)p(x)$
とおくと
$E(f(x)) = \frac{1}{N}\sum^{N}_{i = 0}\frac{g(x_i)}{p(x_i)}$
なので、放射輝度をp(x)で割る必要があったのですね。計算式の項を一つ入れ替えただけでここまで結果が変わるなんて、どうしたらいいのか分からず悩んでた時間はなんだったんだ！と叫びたい思いです。
ただ、ここで一つ疑問があります。確率密度関数p(x)は、どのような方法で決定されるのでしょうか。RRTでは、直接照明の計算において、p(x)を光源（面光源）の面積の逆数を指定するのが良いとしていますが、それ以上の突っ込んだ説明はなされていません。つまるところ、直接照明と間接照明ごとに、ある一定条件を満たす関数を適当に自分で見つけろということでしょうか。

2008-04-15

対応フォーマットの件

CG プログラミング

決定を保留していた対応フォーマットの件ですが、正直、ベンダー独自のフォーマットの仕様を一から理解するのが面倒くさいので、とりあえず既に知っていてるPLYに対応することにします。実装に関しては、スタンフォード大学がCのライブラリを公開しているので、そのまま本体に組み込む予定です。ただし、PLYはあくまで即席対応という認識であって、将来的には、既存のシーン記述言語（RIBまたはPOV-Ray）を標準フォーマットととし、これに一本化したいところ*1。バッドノウハウであるファイルフォーマット関連の作業は、さっさと対応を済ませて、肝心のアルゴリズムの方に注力していきたいです。
実は他のレンダラーを参考にしようと、http://nyaxtstep.com/projects/nytrを調べてみたのですが、独自フォーマットのテキストデータを正規表現を駆使してパージングしていたのはちょっと驚きました。XMLのようにベンダーに依存しない業界標準的なフォーマットが存在しないため、皆さん悩んでらっしゃるのでしょうか*2。
シーンのインポート機能と現在後回しにしている部分の実装（層化サンプリングとか乱数生成器とか）が終わったら、ひとまず公開したいと思います。その後どのように展開していくかは追々考えるということで。

*1:DSLの代わりに汎用のスクリプト言語で記述できると、表現方法の自由度が高くなってうれしいかも。ただし、インタープリタを内蔵するか単なる中間フォーマットへのトランスレータに留めるかは迷うところです。

*2:VRMLの後継として国際標準を目指したX3Dはどうなんでしょうか？

2008-04-14

「おもてなしの経営学」

読書

本書を読んで私が一番引っかかったのは、「ビジネスが分かる」とは、いったいどういう能力なのかということ。
一般的に日本のエンジニアは、技術に関しては興味が尽きないがビジネスには全く疎い人達が多いと言われます。本書でも、日本のソフトウェア業界の事情をそのように評して、日米間の技術者の違いを論じながら、今は「ビジネスが分かる技術者」が必要なんだと著者の中島さんは主張します。過去にMS本社でWindows95やIEの開発に直接携わった方ゆえにかなり説得力があります。そして同時に、明言はされないものの、日本のエンジニアのふがいなさに対する中島さんの強い苛立ちと挑発を感じずにはいられませんでした。「おまえら、もっと金儲けに貪欲になれ」と。
中島さんが暗に指摘するように、私自身いつまでもナイーブなままでいるつもりはないですし、一生技術者としてやっていきたいと本気で考えるならば、お金の事から絶対に目を背けることはできないと思っています*1。自分の力で稼ぐことができなければ、組織内で自分の意見を貫くことが困難になり、その結果、将来やりたい事が見つかっても仕事として挑戦しずらい状況に追い込まれるからです。
ところが、そこから一歩進んで「ビジネスが分かるようになるには、どのような能力が必要なのか」を考えようとすると、自分は途端に思考が停まってしまう傾向があります。技術に関する場合、分からない事に直面してもどのように学習していけばいいか、これまでの経験からあまり悩むことはないですが、殊にビジネスとなると、どこから手を付けていっていいか、なかなか見当がつきません。残念ながら、本書からこの疑問に対する回答は得られませんでした。
もしも「ギークのためのビジネス入門」みたいな本があったらぜひ読んでみたいです。最初に何から取り組めばいいか教えてくれる手引き書となる本。ソフトウェア関連の技術書に例えると、「プログラミング作法」とか「達人プログラマー」のような本がそれに相当するのでしょうか。本来、ビジネス感覚というものは読書だけで身に付くものではなく、現実の仕事を通して培われていくものだと思いますが、私のように手探りな状況の人間にとっては心強い一冊になってくれると思います。

おもてなしの経営学アップルがソニーを超えた理由 (アスキー新書)

作者: 中島聡
出版社/メーカー: アスキー
発売日: 2008/03/10
メディア: 新書
購入: 12人クリック: 531回
この商品を含むブログ (270件) を見る

*1:実際のところ、自分の中にまだ素朴な技術者像に惹かれる部分があって、それがビジネス的な発想を阻んでいるところに自分の未熟さを感じています。