穴日記

どうだ明るくなったろう

AO BenchをGo言語に移植しました

AO Bench(http://code.google.com/p/aobench/)というAmbient Occlusionによるレンダリングを行うベンチマークがあり、結構いろんな言語に移植されてるのですが、Go言語版がなかったので移植しました。
Core2の2.66GHzで7.8 secといったところでした。前のOCaml版よりちょっと遅いですね。
ただ、非常にべたな移植であんまりGo言語の癖を把握してないので意図せず遅い感じになってる部分もあるかもしれません。
あと、本当なら並列実行とか並行実行とかの仕組みを取り入れるべきでしょう。オブジェクト指向とかも。ただあんまり改変しすぎると俺Benchになってしまいますね。
また、ファイル出力がようわからなかったので普通に標準出力にppmの結果を出してたりしてちょっとアレです。

package main

import "math"
import "fmt"
import "math/rand"

const (
      WIDTH  = 256
      HEIGHT = 256
      NSUBSAMPLES = 2
      NAO_SAMPLES = 8
)

type vec struct {
     x, y, z float64     
}

type Isect struct {
     t float64
     p, n vec
     hit int
}

type Sphere struct {
     center vec
     radius float64
}

type Plane struct {
     p, n vec
}

type Ray struct {
     org, dir vec
}

var spheres [3]Sphere
var plane Plane

func vdot(v0, v1 vec) float64 {
     return v0.x * v1.x + v0.y * v1.y + v0.z * v1.z
}

func vcross(c *vec, v0, v1 vec) {
      c.x = v0.y * v1.z - v0.z * v1.y;
      c.y = v0.z * v1.x - v0.x * v1.z;
      c.z = v0.x * v1.y - v0.y * v1.x;
}

func vnormalize(c *vec) {
     length := math.Sqrt(vdot(*c, *c)) 
     if math.Abs(length) > 1.0e-17 {
     	c.x /= length
	c.y /= length
	c.z /= length
     }
}

func ray_sphere_intersect(isect *Isect, ray *Ray, sphere *Sphere) {
     rs := vec{ray.org.x - sphere.center.x,
     	       ray.org.y - sphere.center.y,  
	       ray.org.z - sphere.center.z}
     B := vdot(rs, ray.dir) 
     C := vdot(rs, rs) - sphere.radius * sphere.radius
     D := B * B - C
     if D > 0.0 {
     	t := -B - math.Sqrt(D)
	if (t > 0.0) && (t < (*isect).t) {
	   isect.t = t
	   isect.hit = 1
	   
	   isect.p.x = ray.org.x + ray.dir.x * t
	   isect.p.y = ray.org.y + ray.dir.y * t
	   isect.p.z = ray.org.z + ray.dir.z * t

	   isect.n.x = isect.p.x - sphere.center.x
	   isect.n.y = isect.p.y - sphere.center.y
	   isect.n.z = isect.p.z - sphere.center.z

	   vnormalize(&(isect.n))
	}
     }
}

func ray_plane_intersect(isect *Isect, ray *Ray, plane *Plane) {
     d := -vdot(plane.p, plane.n)
     v := vdot(ray.dir, plane.n)
     if math.Abs(v) < 1.0e-17 {
     	return
     }
     t := -(vdot(ray.org, plane.n) + d) / v
     if (t > 0.0) && (t < isect.t) {
     	isect.t = t
	isect.hit = 1
	
	isect.p.x = ray.org.x + ray.dir.x * t
	isect.p.y = ray.org.y + ray.dir.y * t
	isect.p.z = ray.org.z + ray.dir.z * t
	
	isect.n = plane.n
     }
}

func orthoBasis(basis *[3]vec, n vec) {
     basis[2] = n
     basis[1] = vec{0.0, 0.0, 0.0}
     if (n.x < 0.6) && (n.x > -0.6) {
     	basis[1].x = 1.0
     } else if (n.y < 0.6) && (n.y > -0.6) {
        basis[1].y = 1.0 
     } else if (n.z < 0.6) && (n.z > -0.6) {
        basis[1].z = 1.0
     } else {
        basis[1].x = 1.0
     }
     vcross(&basis[0], basis[1], basis[2])
     vnormalize(&basis[0])
     vcross(&basis[1], basis[2], basis[0])
     vnormalize(&basis[1])
}

func ambient_occlusion(col *vec, isect *Isect) {
     ntheta := NAO_SAMPLES
     nphi   := NAO_SAMPLES
     const eps float64 = 0.0001
     p := vec{isect.p.x + eps * isect.n.x,
              isect.p.y + eps * isect.n.y,
              isect.p.z + eps * isect.n.z}
     var basis [3]vec
     orthoBasis(&basis, isect.n)     

     var occlusion float64 = 0.0

     for j := 0; j < ntheta; j ++ {
     	 for i := 0; i < nphi; i ++ {
	     theta := math.Sqrt(rand.Float64())
	     phi   := 2.0 * math.Pi * rand.Float64()
	     x := math.Cos(phi) * theta
	     y := math.Sin(phi) * theta
	     z := math.Sqrt(1.0 - theta * theta)
	     
	     rx := x * basis[0].x + y * basis[1].x + z * basis[2].x
	     ry := x * basis[0].y + y * basis[1].y + z * basis[2].y
	     rz := x * basis[0].z + y * basis[1].z + z * basis[2].z

	     var ray Ray
	     ray.org = p
	     ray.dir = vec{rx, ry, rz}
	     
	     var occIsect Isect
	     occIsect.t = 1.0e+17
	     occIsect.hit = 0
	     ray_sphere_intersect(&occIsect, &ray, &spheres[0])	
	     ray_sphere_intersect(&occIsect, &ray, &spheres[1])
	     ray_sphere_intersect(&occIsect, &ray, &spheres[2])
	     ray_plane_intersect(&occIsect, &ray, &plane)
	     if occIsect.hit > 0 {
	     	occlusion += 1.0
	     }
        }
     }
     var tmp float64 = float64(ntheta * nphi)
     occlusion = (tmp - occlusion) / tmp
     *col = vec{occlusion, occlusion, occlusion}
}

func clamp(f float64) uint8 {
     i := int32(f * 255.5)
     if i < 0 {
     	i = 0
     }
     if i > 255 {
        i = 255
     }
     return uint8(i)
}

func render(img []uint8, w int32, h int32, nsubsamples int32) {
     var fimg []float64 = make([]float64, w * h * 3)
     var x, y, u, v int32
     for y = 0; y < h; y ++ {
     	 for x = 0; x < w; x ++ {
	     for v = 0; v < nsubsamples; v ++ {
	     	 for u = 0; u < nsubsamples; u ++ {
		     px := (float64(x) + (float64(u) / float64(nsubsamples)) - (float64(w) / 2.0)) / (float64(w) / 2.0)
		     py := -(float64(y) + (float64(v) / float64(nsubsamples)) - (float64(h) / 2.0)) / (float64(h) / 2.0)
		     
		     var ray Ray
		     ray.org = vec{0.0, 0.0, 0.0}
		     ray.dir = vec{px, py, -1.0}
		     vnormalize(&ray.dir)
		     
		     var isect Isect
		     isect.t = 1.0e+17
		     isect.hit = 0
		     
		     ray_sphere_intersect(&isect, &ray, &spheres[0])
		     ray_sphere_intersect(&isect, &ray, &spheres[1])
		     ray_sphere_intersect(&isect, &ray, &spheres[2])
		     ray_plane_intersect(&isect, &ray, &plane)

		     if isect.hit > 0 {
		     	var col vec
			ambient_occlusion(&col, &isect)
			fimg[3 * (y * w + x) + 0] += col.x
			fimg[3 * (y * w + x) + 1] += col.y
			fimg[3 * (y * w + x) + 2] += col.z
		     }
                  }
              }
	      fimg[3 * (y * w + x) + 0] /= float64(nsubsamples * nsubsamples)
	      fimg[3 * (y * w + x) + 1] /= float64(nsubsamples * nsubsamples)
	      fimg[3 * (y * w + x) + 2] /= float64(nsubsamples * nsubsamples)

	      img[3 * (y * w + x) + 0] = clamp(fimg[3 * (y * w + x) + 0])
	      img[3 * (y * w + x) + 1] = clamp(fimg[3 * (y * w + x) + 1])
	      img[3 * (y * w + x) + 2] = clamp(fimg[3 * (y * w + x) + 2])
        }
     }
}

func init_scene() {
     spheres[0].center = vec{-2.0, 0.0, -3.5}
     spheres[0].radius = 0.5
     spheres[1].center = vec{-0.5, 0.0, -3.0}
     spheres[1].radius = 0.5
     spheres[2].center = vec{ 1.0, 0.0, -2.2}
     spheres[2].radius = 0.5
     plane.p = vec{0.0, -0.5, 0.0}
     plane.n = vec{0.0,  1.0, 0.0}
}

func saveppm(fname string, w int32, h int32, img []uint8) {
     fmt.Printf("P3\n")
     fmt.Printf("%d %d\n", w, h)
     fmt.Printf("255\n")
     for i := 0; i < WIDTH*HEIGHT; i ++ {
       fmt.Printf("%d %d %d\n", img[i*3], img[i*3+1], img[i*3+2])
     }     
}

func main() {
     var img []uint8 = make([]uint8, WIDTH * HEIGHT * 3)
     init_scene()
     render(img[:], WIDTH, HEIGHT, NSUBSAMPLES)

     saveppm("ao.ppm", WIDTH, HEIGHT, img[:])
}