## C# vs C++ 全局照明渲染性能比试

2010-06-23 09:48 by Milo Yip, ... 阅读, ... 评论, 收藏, 编辑

512x512像素，每像素1000采样，C#版本渲染时间为40分47秒

## 实验内容

• 使用蒙地卡罗路径追踪(Monte Carlo path-tracing)来产生全局照明效果
• 支持三种双向反射分布函数(bidirectional reflectance distribution function, BRDF): 镜射(specular)、漫射(diffuse)和玻璃(即纯折射的介质)
• 使用2x2超采样(super-sampling)去实现反锯齿
• 使用OpenMP作并行运算，充份利用多核性能

## C++版本

#include <math.h>   // smallpt, a Path Tracer by Kevin Beason, 2008
#include <stdlib.h> // Make : g++ -O3 -fopenmp smallpt.cpp -o smallpt
#include <stdio.h>  //        Remove "-fopenmp" for g++ version < 4.2
#include <time.h>		// MILO
#include "erand48.inc"	// MILO
#define M_PI 3.141592653589793238462643	// MILO
struct Vec {        // Usage: time ./smallpt 5000 && xv image.ppm
double x, y, z;                  // position, also color (r,g,b)
Vec(double x_=0, double y_=0, double z_=0){ x=x_; y=y_; z=z_; }
Vec operator+(const Vec &b) const { return Vec(x+b.x,y+b.y,z+b.z); }
Vec operator-(const Vec &b) const { return Vec(x-b.x,y-b.y,z-b.z); }
Vec operator*(double b) const { return Vec(x*b,y*b,z*b); }
Vec mult(const Vec &b) const { return Vec(x*b.x,y*b.y,z*b.z); }
Vec& norm(){ return *this = *this * (1/sqrt(x*x+y*y+z*z)); }
double dot(const Vec &b) const { return x*b.x+y*b.y+z*b.z; } // cross:
Vec operator%(const Vec &b){return Vec(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);}
};
struct Ray { Vec o, d; Ray(const Vec &o_, const Vec &d_) : o(o_), d(d_) {} };
enum Refl_t { DIFF, SPEC, REFR };  // material types, used in radiance()
struct Sphere {
Vec p, e, c;      // position, emission, color
Refl_t refl;      // reflection type (DIFFuse, SPECular, REFRactive)
Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_):
double intersect(const Ray &r) const { // returns distance, 0 if nohit
Vec op = p-r.o; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
if (det<0) return 0; else det=sqrt(det);
return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0);
}
};
Sphere spheres[] = {//Scene: radius, position, emission, color, material
Sphere(1e5, Vec( 1e5+1,40.8,81.6), Vec(),Vec(.75,.25,.25),DIFF),//Left
Sphere(1e5, Vec(-1e5+99,40.8,81.6),Vec(),Vec(.25,.25,.75),DIFF),//Rght
Sphere(1e5, Vec(50,40.8, 1e5),     Vec(),Vec(.75,.75,.75),DIFF),//Back
Sphere(1e5, Vec(50,40.8,-1e5+170), Vec(),Vec(),           DIFF),//Frnt
Sphere(1e5, Vec(50, 1e5, 81.6),    Vec(),Vec(.75,.75,.75),DIFF),//Botm
Sphere(1e5, Vec(50,-1e5+81.6,81.6),Vec(),Vec(.75,.75,.75),DIFF),//Top
Sphere(16.5,Vec(27,16.5,47),       Vec(),Vec(1,1,1)*.999, SPEC),//Mirr
Sphere(16.5,Vec(73,16.5,78),       Vec(),Vec(1,1,1)*.999, REFR),//Glas
Sphere(600, Vec(50,681.6-.27,81.6),Vec(12,12,12),  Vec(), DIFF) //Lite
};
inline double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; }
inline int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); }
inline bool intersect(const Ray &r, double &t, int &id){
double n=sizeof(spheres)/sizeof(Sphere), d, inf=t=1e20;
for(int i=int(n);i--;) if((d=spheres[i].intersect(r))&&d<t){t=d;id=i;}
return t<inf;
}
Vec radiance(const Ray &r, int depth, unsigned short *Xi){
double t;                               // distance to intersection
int id=0;                               // id of intersected object
if (!intersect(r, t, id)) return Vec(); // if miss, return black
const Sphere &obj = spheres[id];        // the hit object
Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c;
double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
if (++depth>5) if (erand48(Xi)<p) f=f*(1/p); else return obj.e; //R.R.
if (depth > 100) return obj.e; // MILO
if (obj.refl == DIFF){                  // Ideal DIFFUSE reflection
double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
} else if (obj.refl == SPEC)            // Ideal SPECULAR reflection
Ray reflRay(x, r.d-n*2*n.dot(r.d));     // Ideal dielectric REFRACTION
bool into = n.dot(nl)>0;                // Ray from outside going in?
double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0)    // Total internal reflection
Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ?   // Russian roulette
}
int main(int argc, char *argv[]){
clock_t start = clock(); // MILO
int w=512, h=512, samps = argc==2 ? atoi(argv[1])/4 : 250; // # samples
Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm()); // cam pos, dir
Vec cx=Vec(w*.5135/h), cy=(cx%cam.d).norm()*.5135, r, *c=new Vec[w*h];
#pragma omp parallel for schedule(dynamic, 1) private(r)       // OpenMP
for (int y=0; y<h; y++){                       // Loop over image rows
fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samps*4,100.*y/(h-1));
unsigned short Xi[3]={0,0,y*y*y}; // MILO
for (unsigned short x=0; x<w; x++)   // Loop cols
for (int sy=0, i=(h-y-1)*w+x; sy<2; sy++)     // 2x2 subpixel rows
for (int sx=0; sx<2; sx++, r=Vec()){        // 2x2 subpixel cols
for (int s=0; s<samps; s++){
double r1=2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1);
double r2=2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2);
Vec d = cx*( ( (sx+.5 + dx)/2 + x)/w - .5) +
cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + cam.d;
} // Camera rays are pushed ^^^^^ forward to start in interior
c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25;
}
}
printf("\n%f sec\n", (float)(clock() - start)/CLOCKS_PER_SEC); // MILO
FILE *f = fopen("image.ppm", "w");         // Write image to PPM file.
fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);
for (int i=0; i<w*h; i++)
fprintf(f,"%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z));
}


## C#版本

using System;
using System.IO;
namespace smallpt_cs {
struct Vec {        // Usage: time ./smallpt 5000 && xv image.ppm
public double x,y,z;                 // position,also color (r,g,b)
public Vec(double x_,double y_,double z_) {x=x_;y=y_;z=z_;}
public static Vec operator +(Vec a,Vec b) {return new Vec(a.x+b.x,a.y+b.y,a.z+b.z);}
public static Vec operator -(Vec a,Vec b) {return new Vec(a.x-b.x,a.y-b.y,a.z-b.z);}
public static Vec operator *(Vec a,double b) {return new Vec(a.x*b,a.y*b,a.z*b);}
public Vec mult(Vec b) { return new Vec(x*b.x,y*b.y,z*b.z);}
public Vec norm() { return this=this*(1/Math.Sqrt(x*x+y*y+z*z));}
public double dot(Vec b) { return x*b.x+y*b.y+z*b.z;}//cross:
public static Vec operator %(Vec a,Vec b) { return new Vec(a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x);}
}
enum Refl_t { DIFF,SPEC,REFR }; // material types,used in radiance()
struct Ray { public Vec o,d;public Ray(Vec o_,Vec d_) { o=o_;d=d_;} }
class Sphere {
public Vec p,e,c;     // position,emission,color
public Refl_t refl;     // reflection type (DIFFuse,SPECular,REFRactive)
public Sphere(double rad_,Vec p_,Vec e_,Vec c_,Refl_t refl_) {
}
public double intersect(Ray r)
{ // returns distance,0 if nohit
Vec op=p-r.o;// Solve t^2*d.d+2*t*(o-p).d+(o-p).(o-p)-R^2=0
if (det<0) return 0;else det=Math.Sqrt(det);
return (t=b-det) > eps?t : ((t=b+det) > eps?t : 0);
}
};
class smallpt {
static Random random=new Random();
static double erand48() { return random.NextDouble();}
new Sphere(1e5,new Vec( 1e5+1,40.8,81.6), new Vec(),new Vec(.75,.25,.25),Refl_t.DIFF),//Left
new Sphere(1e5,new Vec(-1e5+99,40.8,81.6),new Vec(),new Vec(.25,.25,.75),Refl_t.DIFF),//Rght
new Sphere(1e5,new Vec(50,40.8,1e5),      new Vec(),new Vec(.75,.75,.75),Refl_t.DIFF),//Back
new Sphere(1e5,new Vec(50,40.8,-1e5+170), new Vec(),new Vec(),           Refl_t.DIFF),//Frnt
new Sphere(1e5,new Vec(50,1e5,81.6),      new Vec(),new Vec(.75,.75,.75),Refl_t.DIFF),//Botm
new Sphere(1e5,new Vec(50,-1e5+81.6,81.6),new Vec(),new Vec(.75,.75,.75),Refl_t.DIFF),//Top
new Sphere(16.5,new Vec(27,16.5,47),      new Vec(),new Vec(1,1,1)*.999, Refl_t.SPEC),//Mirr
new Sphere(16.5,new Vec(73,16.5,78),      new Vec(),new Vec(1,1,1)*.999, Refl_t.REFR),//Glas
new Sphere(600,new Vec(50,681.6-.27,81.6),new Vec(12,12,12), new Vec(),  Refl_t.DIFF) //Lite
};
static double clamp(double x) { return x<0?0 : x > 1?1 : x;}
static int toInt(double x) { return (int)(Math.Pow(clamp(x),1 / 2.2)*255+.5);}
static bool intersect(Ray r,ref double t,ref int id) {
double d,inf=t=1e20;
for (int i=spheres.Length-1;i >= 0;i--)
if ((d=spheres[i].intersect(r)) != 0 && d<t) { t=d;id=i;}
return t<inf;
}
static Vec radiance(Ray r,int depth) {
double t=0;                              // distance to intersection
int id=0;                              // id of intersected object
if (!intersect(r,ref t,ref id)) return new Vec();// if miss,return black
Sphere obj=spheres[id];       // the hit object
Vec x=r.o+r.d*t,n=(x-obj.p).norm(),nl=n.dot(r.d)<0?n:n*-1,f=obj.c;
double p=f.x>f.y&&f.x>f.z?f.x:f.y>f.z?f.y:f.z;//max refl
if (++depth > 5) if (erand48()<p) f=f*(1 / p);else return obj.e;//R.R.
if (depth > 100) return obj.e;
if (obj.refl == Refl_t.DIFF) {                  // Ideal DIFFUSE reflection
double r1=2*Math.PI*erand48(),r2=erand48(),r2s=Math.Sqrt(r2);
Vec w=nl,u=((Math.Abs(w.x)>.1?new Vec(0,1,0):new Vec(1,0,0))%w).norm(),v=w%u;
Vec d=(u*Math.Cos(r1)*r2s+v*Math.Sin(r1)*r2s+w*Math.Sqrt(1-r2)).norm();
}
else if (obj.refl == Refl_t.SPEC)            // Ideal SPECULAR reflection
Ray reflRay=new Ray(x,r.d-n*2*n.dot(r.d));//IdealdielectricREFRACTION
bool into=n.dot(nl) > 0;               // Ray from outside going in?
double nc=1,nt=1.5,nnt=into?nc / nt : nt / nc,ddn=r.d.dot(nl),cos2t;
if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0)    //Total internal reflection
Vec tdir=(r.d*nnt-n*((into?1:-1)*(ddn*nnt+Math.Sqrt(cos2t)))).norm();
double a=nt-nc,b=nt+nc,R0=a*a/(b*b),c=1-(into?-ddn:tdir.dot(n));
double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
return obj.e+f.mult(depth > 2?(erand48()<P?  // Russian roulette
}
public static void Main(string[] args) {
DateTime start=DateTime.Now;
int w=256,h=256,samps=args.Length==2?int.Parse(args[1])/4:25;// # samples
Ray cam=new Ray(new Vec(50,52,295.6),new Vec(0,-0.042612,-1).norm());//cam pos,dir
Vec cx=new Vec(w*.5135/h,0,0),cy=(cx%cam.d).norm()*.5135,r;Vec[] c=new Vec[w*h];
for (int y=0;y<h;y++) {                        // Loop over image rows
Console.Write("\rRendering ({0}spp) {1:F2}%",samps*4,100.0*y/(h-1));
for (int x=0;x<w;x++)   // Loop cols
for (int sy=0,i=(h-y-1)*w+x;sy<2;sy++)     // 2x2 subpixel rows
for (int sx=0;sx<2;sx++) {               // 2x2 subpixel cols
r=new Vec();
for (int s=0;s<samps;s++) {
double r1=2*erand48(),dx=r1<1?Math.Sqrt(r1)-1:1-Math.Sqrt(2-r1);
double r2=2*erand48(),dy=r2<1?Math.Sqrt(r2)-1:1-Math.Sqrt(2-r2);
Vec d=cx*(((sx+.5+dx)/2+x)/w-.5)+
cy*(((sy+.5+dy)/2+y)/h-.5)+cam.d;
} // Camera rays are pushed ^^^^^ forward to start in interior
c[i]=c[i]+new Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25;
}
}
Console.WriteLine("\n{0} sec",(DateTime.Now-start).TotalSeconds);
using (StreamWriter sw=new StreamWriter("image.ppm")) {
sw.Write("P3\r\n{0} {1}\r\n{2}\r\n",w,h,255);
for (int i=0;i<w*h;i++)
sw.Write("{0} {1} {2}\r\n",toInt(c[i].x),toInt(c[i].y),toInt(c[i].z));
sw.Close();
}
}
}
}


Vec和Ray需要不断在计算中产生实例，所以设它们为struct，struct在C#代表值类型(value type)，ibpp在堆栈上高效分配内存的，不需使用GC。渲染时，Sphere是只读对象，因此用class作为引用类型(reference type)去避免不必要的复制。

## 实验结果和分析

 测试版本 需时(秒) (a) C++ 45.548 (b) C# 61.044 (c) C++ SIMD 20.500 (d) C++(OpenMP) 7.397 (e) C++ SIMD(OpenMP) 3.470 (f)* C++ LCG 17.365 (g)* C# LCG 59.623 (h)* C++ LCG (OpenMP) 3.427

*2010/6/23 加入(f)(g)(h)，見更新1

## 更新

1. 2010/6/23: 园友嗷嗷在本文#78楼回覆中指出，C++版本的很大部分开销在于erand48()函数里调用Runtime库函数。所以，决定用简单的LCG随机数实现，取代原来的库函数(包括C++和C#)，再测试三个版本(f)(g)(h)。结果: C++版本(f)比C#版本(g)快2.43倍。 使用OpenMP(h)是没用OpenMP(f)版本的5.06倍。此LCG版本代码可于此下载。再次感谢这位园友。