2014年2月20日木曜日

フォトンマップのための八分木の実装


フォトンマップとは?

今までの記事では、パストレーシングを利用してレンダリングしてきましたが、一番の問題は、なんといってもレンダリングに時間がかかりすぎることです。フォトンマップを使ったレンダリングは、この問題を解決してくれます。フォトンマップとは、光源から放射されたフォトン(レイ、光線)が拡散面に衝突した位置(と入射方向や反射方向)とその衝突面の材質を格納するデータ構造です。このフォトンマップを使えば、拡散面における輝度推定が容易になります。つまり、フォトンマップを使ったレンダリング過程は、大きくわけると2段階にわけられます。第1段階は、フォトンマップの構築です。第2段階は、レイトレーシングまたはパストレーシングを使って、拡散面の輝度推定値を構築されたフォトンマップから計算してレンダリングします。

フォトンマップのデータ構造

拡散面に衝突したフォトンの情報が格納されてれば、どんなデータ構造を使おうが問題ありません。ただ、輝度推定の過程で、格納された全てのフォトンの中から目的のフォトンを効率よく見つけ出す必要があるので、この目的に合ったデータ構造を使えば、より短い時間でレンダリングできます。フォトンマップによる輝度推定では、拡散面上の対象の点に近い位置にあるフォトンを集めて輝度を計算します。このような計算に対して、フォトンを効率よく見つけるデータ構造をいろいろ考えだされているようですが、ここでは、パストレーシングに比べてかなり早くなっていれば、他のデータ構造に比べてちょっと遅いくらいは問題ないので、比較的単純な構造である八分木を実装することにします。二分木の3次元版みたいなものだから、イメージしやすいし。

八分木の仕様

八分木は空間を再帰的に8分割して、その子空間にデータを格納する木構造です。この木構造をフォトンマップとして使えるように、次のような仕様にしました。
  • 空間として直方体を使う
  • 直方体の内側に位置するフォトンは、その直方体に格納される
  • 1つの直方体はフォトンに対して容量を持っていて、フォトンを格納できる数が決まっている
  • 直方体のフォトンが容量に達したら、直方体を8分割し、該当する子の直方体にフォトンを格納する

八分木の実装

八分木クラスの主要なメソッドを概観します。

コンストラクタ

コンストラクタでは、直方体を定義します。
  abstract class Octree[B, C <: Octree[B, C]](
    val xmin: Float, val xmax: Float,
    val ymin: Float, val ymax: Float,
    val zmin: Float, val zmax: Float,
    val capacity: Int,
    val parent: Option[C]
  )
型パラメータBは、格納するデータの型です。型パラメータCは、サブクラスの型です。子の直方体を生成する際には、C型のインスタンスを生成するようにします。これは、サブクラスで定義したメンバーを使用できるようにするためです。

子を生成するメソッド

直方体を8分割して子の直方体を生成するときに使用されるメソッドです。コンストラクタと同じ引数を取ります。サブクラスで実装します。
    protected def create(
      xmin: Float, xmax: Float,
      ymin: Float, ymax: Float,
      zmin: Float, zmax: Float,
      capacity: Int,
      parent: Option[C]
    ): C
戻り値の型は、上述したようにサブクラスの型です。

フォトンの追加

フォトンを追加するメソッドです。既に最大容量に達していれば、子を作成し、子に格納する。
10 
11 
    def add( a: Vector3, b: B )(implicit tag: ClassTag[C]): Unit = if( ndata < capacity ) {
      _data = (a, b) :: _data
      ndata += 1
      onUpdate( b )
    } else {
      createChildren // 既に子を持っている場合は子を作成しない
      // このOctreeが持っているデータを子に移す
      _data.foreach { case (p, q) => findChild( p ).foreach( _.add( p, q ) ) }
      findChild( a ).foreach( _.add( a, b ) )
      _data = Nil
    }

フォトンの探索

指定した球に含まれる、または接している最下層の直方体を見つけ、クロージャに渡す。
10 
11 
    def filterChildren( b: BoundingSphere, f: (Int, C) => Unit ): Unit = boundingSphere.contains( b ) match {
      // 離れている場合
      case 3 => ()
      // bがこのOctreeを含んでいる場合
      case 2 => f( 2, this )
      // このOctreeがbを含んでいるか接している場合
      case n => _children match {
        case None => f( n, this )
        case Some( cs ) => cs.foreach { _.filterChildren( b, f ) }
      }
    }

2014年2月13日木曜日

魚眼レンズの実装

大域照明で3次元モデルをレンダリングする場合、モデル中のある点における輝度を計算するために、図1のようにその点を中心とする単位半球面を通過する光の量を求める必要があります。つまり、その点からモデル全体を見渡した時に、その点から見て光源はどの方向にありどのくらい入射してくるのか、また光源じゃない部分からもその点に間接的に光は届くので、その部分からの光の量はどのくらいなのかということを計算しなければいけません。前回のパストレーシングは、ある点に光源から直接届く光と光源から1度別の点で反射して届く光のみを考慮して輝度を計算しました。したがって、完全な大域照明と言えないでしょう。

図1 ある点に入射する光。球面に沿って面積分して光の総和を求める

モデル中のある点に入射してくる光の量というのは、言い換えれば、その点を中心とする単位半球面に沿った光の量の面積分のことだと言えるでしょう。パストレーシングは、その面積分を行っているわけですね。つまり、直接光に関して言えば、半球面上に投影されるモデル中の光源の面積が直接光の量に比例します。であるからには、ある点からモデル全体を見渡した時、その半球面上に描かれるモデルを視覚化してみたくなりました。そこで、魚眼レンズのようなものを実装し、魚眼レンズに移った像を平面に描いてみました。球面を平面で表現する方法は、地球を地図で表現する方法と同じようにすれば良さそうですね。学生のころ地理の授業でいろいろならった記憶があります。メルカトル図法とか。ここでは、ランベルト正積方位図法を使って2次元に映し出してみることにしました。ランベルト正積方位図法っていうのは、面積比が実際と比べてだいたい同じくらいになる図法です。正確に面積比が同じなるわけではないようです。ですが、簡単に球面上の面積を視覚化するには、ちょうどいい方法だと思い、この図法を使うことにしました。

平面(スクリーン)上のピクセル値から球面座標を求める

球面上の点をスクリーンに投影するには、球面上の点とスクリーン上のピクセル値を結び付けなければなりません。図2に示すように、ランベルト正積方位図法では、球面上の点P'がスクリーン上に投影する点P''は、天頂Aと点P'が作る線分AP'を半径とする円とスクリーン2と交わる点のうち近い方の点とします。なので、スクリーン上のピクセル値を決めれば、その点に対応する球面上の点を決めることができる。


図2 ランベルト方位図法

球面上の点は、極座標で表現するのが簡単なので、天頂角AOP'を求めるのがよいでしょう。この天頂角AOP'は、余弦定理で計算できます。そんな定理ありましたね。懐かしいですね。ほとんど忘れてるのでWikipediaで調べてみました。
cos(x) = (b^2 + c^2 - a^2) / (2bc)
三角形の三辺の長さが分かれば、角度が分かるよっていう定理ですね。
ここで三角形AOP'を見ると、辺AP'は分かってます。でも、辺AOと辺OP'は分かってないですね。ただ、辺AOと辺OP'は円の半径OBなので同じなので、求めるべきは円の半径OBということになります。

ところで、点Bは、魚眼レンズの一番端の部分です。この点をスクリーン上で考えると、やはりスクリーンの一番端になります。スクリーンが正方形だと仮定して、その幅をwとすると線分AB'の長さは、w/2となります。線分AB'と線分ABの長さは等しいので、線分ABの長さもw/2になります。三角形AOBは、直角二等辺三角形なので、各辺の長さの比AO:OB:ABは、1:1:sqrt(2)となります。つまり、辺AOの長さは、辺ABの長さを1/sqrt(2)倍したw/(2 x sqrt(2))です。先ほど述べたように、線分AOと線分P'Oの長さは等しいので、線分P'Oの長さもw/(2 x sqrt(2))です。

これで、三角形AOP'の全ての辺の長さが分かりました。余弦定理に当てはめると天頂角AOP'を求めることができます。次のソースコードは、Scalaでスクリーン上のピクセル値から半球上の点の座標を求めるコードです。

10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
    /** スクリーンの幅 */
    private lazy val s = screen.ymax - screen.ymin

    /** スクリーンの半径 */
    private lazy val sr = s/( 2f * sqrt( 2f ) )

    def position( x: Int, y: Int ) = {
      // ピクセル値からシステム座標系に変換
      val sx = screen.screenX( x )
      val sy = screen.screenY( y )
      // OBの長さの2乗
      val ss = s * s
      // AP''の長さ
      val rr = sx * sx + sy * sy
      val r = sqrt( rr )
      // 天頂角AOP'
      val theta = acos( 1f - 4 * rr / ss )
      // OP'のxy成分
      val xy = sr * sin( theta )
      // 点P'の座標
      if( x == screen.width/&& y == screen.width/) Vector3( 0, 0, sr * cos( theta ) )
      else Vector3( sx * xy / r, sy * xy / r, sr * cos( theta ) )
    }

半球面の点P'に投影されるモデル中の点を求める

スクリーン上のピクセル値から球面上の点P'の座標は求めることができました。では、この点P'にはモデル中のどの点が投影されるのでしょうか? この計算はさほど難しくありません。なぜなら、ベクトルOP'とシーンの中の全てのポリゴンとの交点のうち点P'と最も距離の近いものが、求める点になります。これは、レイトレーシングの時に登場する計算と同じです。

魚眼レンズを使ってレンダリング

前回のシーンを実装した魚眼レンズを使ってレイトレーシングでレンダリングしました。比較の為に普通のカメラを使ってレンダリングした画像を載せておきます。左が普通のカメラ。右が魚眼レンズ。



拡散面から見たモデル全体

下左図の赤丸の中の白い点から魚眼レンズでシーン全体を見渡してみました。下右図は、レイトレーシングのレンダリング結果です。


以下はその時のコード。
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
package examples {

  import net.ruffy.marble.graphics2d.{ Color }
  import net.ruffy.marble.graphics3d.{ Screen, DefaultCamera, LambertCamera, RayTracer }
  import net.ruffy.marble.math.{ Math, Vector3 }
  import Math._

  object LambertMap {

    def main( args: Array[String] ) {

      val renderer = new RayTracer {
        val width = 300
        val height = 300
        val screen = Screen( width, height )
        val root = cornelBox
        val camera = {
          val cam = DefaultCamera( screen, 17f, 88f.toRadians, 90f.toRadians, 5 )
          root.toCamera( cam ).nearestIntersection( cam.createPhoton(width/8, height/2) ) match {
            case Some( i ) =>
              val wp = cam.toWorld( i.point )
              val wq = cam.toWorld( i.point + i.normal )
              LambertCamera( screen, wp, wp -> wq, 0f )
            case _ => LambertCamera( screen, 5f, 88f.toRadians, 90f.toRadians )
          }
        }
        val numberOfTimesToTrace = 5
      }

      renderer.render
      renderer.screen.write( "png""images/raytrace_lambert.png" )

    }

  }

}

2014年2月3日月曜日

パストレーシングの収束2

前回は、数字でパストレーシングの収束を考えました。今回は実際にレンダリングして画像のノイズの減り具合を目で比べてみました。光源から拡散面に1度反射してカメラに到達したレイ(直接光)と光源から拡散面に2度反射してカメラに到達したレイ(環境光)を使って輝度を計算しました。環境光として2度拡散面に反射してカメラに到達したレイしか考慮してないためか、全体的に暗めな印象になってます。各ピクセルの輝度計算に使用したレイは、30000 paths/pixelくらいまでです。かなり時間かかりました。下の画像のように、10000 paths/pixelくらいでノイズが目立たなくなって来ています。