AOBenchをF#で実装した。
http://lucille.atso-net.jp/aobench/ で開催されているAOBenchをF#にて実装しました。
実装としてはC#版からの多くを拝借して移植。
あんまり関数型言語っぽくしてない。画像処理のような大きな配列を書き換える（リストから別のリストを作るのはやりやすいんだけどね）作業は、手続き型的に書いたほうがいいと思う。
ヴァリアントとマッチをうまくつかうと素直な流れで書くことができますね。
速度はCore2Duo T8100 2.1GHzで約9.4秒でした。
以下ソース
#light open System open System.Drawing open System.Drawing.Imaging open System.IO open System.Runtime.InteropServices open System.Diagnostics let IMAGE_WIDTH = 256; let IMAGE_HEIGHT = 256; let NSUBSAMPLES = 2; let NAO_SAMPLES = 8; // # of sample rays for ambient occlusion ///-------------------------------------------------------------------------------------- type Vec(x:double, y:double, z:double) = member v.x = x member v.y = y member v.z = z member v.length = sqrt(x*x+y*y+z*z) let vdot(v0:Vec, v1:Vec) = v0.x * v1.x + v0.y * v1.y + v0.z * v1.z let vadd(v0:Vec, v1:Vec) = Vec(v0.x+v1.x,v0.y+v1.y,v0.z+v1.z) let vsub(v0:Vec, v1:Vec) = Vec(v0.x-v1.x,v0.y-v1.y,v0.z-v1.z) let vcross(v0:Vec, v1:Vec) = Vec(v0.y * v1.z - v0.z * v1.y,v0.z * v1.x - v0.x * v1.z,v0.x * v1.y - v0.y * v1.x) let vnormalize(v:Vec) = let l = 1.0/v.length in Vec(v.x*l,v.y*l,v.z*l) ///-------------------------------------------------------------------------------------- type HitInfo = | Hit of double * Vec * Vec (* distance, position, normal *) | NoHit ///-------------------------------------------------------------------------------------- type Ray(o:Vec,d:Vec) = member r.org = o member r.dir = d member r.pos t = vadd(r.org, Vec(r.dir.x*t,r.dir.y*t,r.dir.z*t)) ///-------------------------------------------------------------------------------------- type IGeometry = abstract intersect: Ray * double -> HitInfo ///-------------------------------------------------------------------------------------- type Sphere(c:Vec,r:double) = member s.c = c member s.r2 = r * r interface IGeometry with member s.intersect(r:Ray, dist:double)= let rs = vsub(r.org,s.c) in let B = vdot(rs,r.dir) in let C = vdot(rs,rs) - s.r2 in let D = B*B-C in if D > 0.0 then let t = -B - sqrt(D) in if (t > 0.0) && (t < dist) then let p = r.pos(t) in Hit(t,p,vsub(p,s.c)) else NoHit else NoHit ///-------------------------------------------------------------------------------------- type Plane(p:Vec,n:Vec) = member s.p = p member s.n = n interface IGeometry with member s.intersect(r:Ray, dist:double)= let d = -vdot(s.p ,n) in let v = vdot(r.dir,n) in if (abs(v) < 1.0e-6) then NoHit else let t = -(vdot(r.org,s.n) + d) / v in if ( (t > 0.) && (t < dist) ) then let p = r.pos(t) in Hit(t,p,n) else NoHit ///-------------------------------------------------------------------------------------- let rnd = System.Random();; let spheres = [Sphere(Vec(-2.0,0.0,-3.5),0.5) ; Sphere(Vec(-0.5, 0.0, -3.0), 0.5) ; Sphere(Vec(1.0, 0.0, -2.2), 0.5)];; let planes = [Plane(Vec(0.0, -0.5, 0.0), Vec(0.0, 1.0, 0.0))];; ///-------------------------------------------------------------------------------------- let intersection_geometry(g:IGeometry, r:Ray, t:double) = g.intersect(r,t); let rec intersection_list(lst:'a list, info:HitInfo, r:Ray, t:double) = match lst with | [] -> info | y::ys -> let h = intersection_geometry(y,r,t) in match h with | Hit(d,p,n) -> intersection_list(ys,h,r,d) | NoHit -> intersection_list(ys,info,r,t) let intersection(r:Ray, t:double) = let h = intersection_list(spheres,NoHit,r,t) in match h with |Hit(d,p,n) -> intersection_list(planes,h,r,d) |NoHit -> intersection_list(planes,h,r,t) ///-------------------------------------------------------------------------------------- let ortho_basis(normal:Vec) = let v = Vec(0., 0., 0.) let orthoBasis = Array.create 3 v orthoBasis.[2] <- normal let tmp = if (normal.x < 0.6) && (normal.x > -0.6) then Vec(1.0,0.0,0.0) elif (normal.y < 0.6) && (normal.y > -0.6) then Vec(0.0,1.0,0.0) elif (normal.z < 0.6) && (normal.z > -0.6) then Vec(0.0,0.0,1.0) else Vec(1.0,0.0,0.0) orthoBasis.[0] <- vnormalize(vcross(tmp, orthoBasis.[2])) orthoBasis.[1] <- vnormalize(vcross(orthoBasis.[2], orthoBasis.[0])) orthoBasis let ambient_occlusion(p:Vec, n:Vec) = let ntheta = NAO_SAMPLES let nphi = NAO_SAMPLES let eps = 0.0001 let rp = Vec(p.x+eps*n.x,p.y+eps*n.y,p.z+eps*n.z) let basis = ortho_basis(n) let mutable count = 0 for j=0 to ntheta-1 do for i=0 to nphi-1 do let theta = Math.Sqrt(rnd.NextDouble()) let phi = 2.0 * Math.PI * rnd.NextDouble() let x = Math.Cos(phi) * theta let y = Math.Sin(phi) * theta let z = Math.Sqrt(1.0 - theta * theta) let dx = x*basis.[0].x + y * basis.[1].x + z * basis.[2].x let dy = x*basis.[0].y + y * basis.[1].y + z * basis.[2].y let dz = x*basis.[0].z + y * basis.[1].z + z * basis.[2].z let rd = Vec(dx,dy,dz) let r = Ray(rp,rd) let h = intersection(r,10000.0) let pc = match h with | Hit(t,p,n) -> 1 | NoHit -> 0 count <- count + pc let occlusionRatio = double(ntheta * nphi - count) / (double)(ntheta * nphi) Vec(occlusionRatio,occlusionRatio,occlusionRatio) let shoot(x:double,y:double) = let t = 10000.0 in let r = Ray(Vec(0.,0.,0.),vnormalize(Vec(x,y,-1.))) in match intersection(r,t) with |Hit(t,p,n) -> ambient_occlusion(p,n) |NoHit -> Vec(0.,0.,0.) ///-------------------------------------------------------------------------------------- let save_image(pngFilePath:string, width:int, height:int, byteImage:byte array) = let bmp = new Bitmap(width, height,PixelFormat.Format24bppRgb) let bd = bmp.LockBits(new Rectangle(new Point(), bmp.Size), ImageLockMode.WriteOnly, bmp.PixelFormat); Marshal.Copy( byteImage, 0, bd.Scan0, bd.Stride * bmp.Height); bmp.UnlockBits(bd); bmp.Save(pngFilePath); ///-------------------------------------------------------------------------------------- let clamp(c:double) = if c < 0.0 then byte(0) elif c > 1.0 then byte(255) else byte(c*255.) let render(byteImage:byte array, width:int, height:int, nsubsamples:int) = let w2 = double(width)*0.5 let h2 = double(height)*0.5 let total = width*height*3 let img = Array.create total 0.0 for y=0 to height-1 do for x=0 to width-1 do for v=0 to nsubsamples-1 do for u=0 to nsubsamples-1 do let px = (double(x) + (double(u) / double(nsubsamples)) - w2)/w2 let py = -(double(y) + (double(v) / double(nsubsamples)) - h2)/h2 let ao = shoot(px,py) let idx = 3*(y*width+x) let r = img.[idx+0] let g = img.[idx+1] let b = img.[idx+2] img.[idx+0] <- r+ao.x img.[idx+1] <- g+ao.y img.[idx+2] <- b+ao.z let total = 1.0/double(nsubsamples*nsubsamples) let idx = 3*(y*width+x) let r = img.[idx+0] let g = img.[idx+1] let b = img.[idx+2] img.[idx+0] <- r*total img.[idx+1] <- g*total img.[idx+2] <- b*total byteImage.[idx+0] <- clamp(img.[idx+0]) byteImage.[idx+1] <- clamp(img.[idx+1]) byteImage.[idx+2] <- clamp(img.[idx+2]) let run() = let w = IMAGE_WIDTH let h = IMAGE_HEIGHT let total = w*h*3 let b0 = byte(0) let img = Array.create total b0 let sw = Stopwatch() sw.Start() render(img, w, h, NSUBSAMPLES) sw.Stop() Console.WriteLine("Elapsed: {0}[sec]", double(sw.ElapsedMilliseconds) / 1000.0) save_image("ao.png",w,h,img) Console.WriteLine("Hit any key....") let c = Console.ReadKey() () ///-------------------------------------------------------------------------------------- run();; ///--------------------------------------------------------------------------------------
