PBRマテリアルでパストレーシングしてみた！(Disney Principled BRDF)

float3 ImportanceSampleGGX( float2 Xi, float Roughness, float3 N )
{
    float a = Roughness * Roughness;
    float Phi = 2 * PI * Xi.x;
    float CosTheta = sqrt( (1 - Xi.y) / ( 1 + (a*a - 1) * Xi.y ) );
    float SinTheta = sqrt( 1 - CosTheta * CosTheta );

    float3 H;
    H.x = SinTheta * cos( Phi );
    H.y = SinTheta * sin( Phi );
    H.z = CosTheta;

    float3 UpVector = abs(N.z) < 0.999 ? float3(0,0,1) : float3(1,0,0);
    float3 TangentX = normalize( cross( UpVector, N ) );
    float3 TangentY = cross( N, TangentX );
    // Tangent to world space
    return TangentX * H.x + TangentY * H.y + N * H.z;
}

#include <math.h>   // Original: smallpt, a Path Tracer by Kevin Beason, 2008
#include <stdlib.h> // Make : g++ -O3 -fopenmp smallpt.cpp -o smallpt
#include <stdio.h>  //   Modified by tyfkda, added PBR material.
struct Vec {
  double x, y, z; Vec(double xyz_=0) { x=y=z=xyz_; }
  Vec(double x_, double y_, double z_){ 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() const { return *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) const {return Vec(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);} };
inline Vec reflect(const Vec& v, const Vec& n) { return v-n*2*n.dot(v); }
inline Vec refract(const Vec& v, const Vec& n, double nnt) {
  double ddn=v.dot(n), cos2t=1-nnt*nnt*(1-ddn*ddn);
  return (v*nnt - n*(ddn*nnt+sqrt(cos2t))).norm(); }
struct Ray { Vec o, d; Ray(const Vec& o_, const Vec& d_) : o(o_), d(d_) {} };
struct Sphere {
  double rad; Vec p, e, c; double roughness, metallic, opacity;
  Sphere(double rad_, Vec p_, Vec e_, Vec c_, double rgh=.2, double metl=0, double opc=1):
    rad(rad_), p(p_), e(e_), c(c_), roughness(rgh), metallic(metl), opacity(opc) {}
  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
    double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad;
    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)),//Left
  Sphere(1e5, Vec( 1e5+99,40.8,81.6),Vec(),Vec(.25,.25,.75), .2, 1),//Rght
  Sphere(1e5, Vec(50,40.8,-1e5),     Vec(),Vec(.75,.75,.75)),//Back
  Sphere(1e5, Vec(50,40.8, 1e5+296), Vec(),Vec()           ),//Frnt
  Sphere(1e5, Vec(50,-1e5, 81.6),    Vec(),Vec(.75,.75,.75)),//Botm
  Sphere(1e5, Vec(50, 1e5+81.6,81.6),Vec(),Vec(.75,.75,.75)),//Top
  Sphere(16.5,Vec(27,16.5,47),       Vec(),Vec(1, .70, .29), .35, 1),//Mirr
  Sphere(16.5,Vec(73,16.5,78),       Vec(),Vec(.5,1,1), .2, 0, .5),//Glas
  Sphere(8,Vec(47,8,110),            Vec(),Vec(1,0,0), .1),
  Sphere(600, Vec(50,681.6-.27,81.6),Vec(12,12,12),  Vec() ,0),//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){
  int n=sizeof(spheres)/sizeof(Sphere); double d, inf=t=1e20;
  for(int i=n;i--;) if((d=spheres[i].intersect(r))&&d<t){t=d;id=i;} return t<inf; }
Vec F(const Vec& F0, const Vec& V, const Vec& H) { // Fresnel reflection
  return F0 + (Vec(1.0) - F0) * pow(1.0 - fmax(V.dot(H), 0.0), 5.0); }  // Schlick approximation
Vec importanceSampleGGX(unsigned short *Xi, double Roughness, const Vec& N) {
  double a = Roughness * Roughness, Phi = 2.0 * M_PI * erand48(Xi), xiy = erand48(Xi);
  double CosTheta = sqrt((1.0 - xiy) / (1.0 + (a * a - 1.0) * xiy));
  double SinTheta = sqrt(1.0 - CosTheta * CosTheta);
  Vec H(SinTheta * cos(Phi), SinTheta * sin(Phi), CosTheta);
  Vec UpVector = abs(N.z) < 0.999 ? Vec(0.0, 0.0, 1.0) : Vec(1.0, 0.0, 0.0);
  Vec TangentX = (UpVector % N).norm(), TangentY = N % TangentX;
  return TangentX * H.x + TangentY * H.y + N * H.z; }
Vec radiance(const Ray &r, int depth, unsigned short *Xi){
  double t; int id=0; if (!intersect(r, t, id)) return Vec(); // if miss, return black
  const Sphere &obj = spheres[id]; ++depth;    // the hit object
  Vec x=r.o+r.d*t, n=(x-obj.p).norm(), c=obj.c;
  if (c.x == .75 && c.y == .75 && c.z == .75 && n.y > -.1) { // Checkered pattern.
    double s=0.1; int tx=fmod((x.x+s*1000)*s,2), ty=fmod(x.y*s,2), tz=fmod(x.z*s,2);
    if (((tx+ty+tz) & 1) == 0) c=c*(1./16); }
  bool into=n.dot(r.d)<0; Vec nl=into?n:n*-1, V=r.d*-1, N, L;
  do { N = importanceSampleGGX(Xi, obj.roughness, nl); L = reflect(r.d, N); } while (L.dot(nl)<=0);
  const Vec &H = N /*==(V + L).norm()*/; double m = obj.metallic, opacity = obj.opacity;
  Vec F0 = Vec(0.04) * (1-m) + c * m, Ks = F(F0, V, H), tdir=L;
  if (opacity < 1.0) {
    double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(N);
    if (1-nnt*nnt*(1-ddn*ddn)>=0) {  // Refraction
      tdir=refract(r.d,N,nnt); double a=nt-nc, b=nt+nc, R0=a*a/(b*b), cc=1-(into?-ddn:tdir.dot(n));
      double Re=R0+(1-R0)*pow(cc,5); Ks = Vec(Re) + (c - Vec(Re)) * (m * opacity);
    } else Ks = Vec(1); /* Total internal reflection */ }
  Vec Kd = Vec(1.0) - Ks; double q = fmax(Ks.x,fmax(Ks.y,Ks.z)), p=1;
  if (erand48(Xi) < q){                   // SPECULAR reflection
    return obj.e + Ks.mult(radiance(Ray(x,L),depth,Xi)) * (1/q);
  } else if (erand48(Xi) < opacity * opacity){  // DIFFUSE reflection
    if (depth>5){p=fmax(c.x,fmax(c.y,c.z))*.99; if (erand48(Xi)>=p) return obj.e;} //R.R.
    Vec d = importanceSampleGGX(Xi, 1.0, n);  // Random sampling for outward hemisphere
    return obj.e + c.mult(radiance(Ray(x,d),depth,Xi)).mult(Kd) * ((1-m)/(p*(1-q)));
  } else {                                // Dielectric REFRACTION
    Vec c2 = Vec(sqrt(c.x), sqrt(c.y), sqrt(c.z));
    if (depth>5){p=fmax(c2.x,fmax(c2.y,c2.z))*.99; if (erand48(Xi)>=p) return obj.e;} //R.R.
    return obj.e + c2.mult(radiance(Ray(x,tdir),depth,Xi)).mult(Kd) * (1/(p*(1-q))); }
}
int main(int argc, char *argv[]){
  int w=1024/2, h=768/2; long samps = argc==2 ? atoi(argv[1]) : 100; // # samples
  Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm()); // cam pos, dir
  Vec cx=Vec(w*.5135/h,0,0), cy=(cx%cam.d).norm()*.5135, *c=new Vec[w*h];
  for (long s=0, ns=8; ++s <= samps; ){
    fprintf(stderr, "\rRendering (%ld/%ld spp) %5.2f%%", s, samps, 100.*s/samps);
#pragma omp parallel for schedule(dynamic, 1)                  // OpenMP
    for (int y=0; y<h; y++){                     // Loop over image rows
      unsigned short Xi[3]; Xi[0]=(-3*y)*((s*131)>>7); Xi[1]=s*s; Xi[2]=y*y*y;
      for (int x=0, i; i=(h-y-1)*w+x, x<w; x++) {           // Loop cols
        Vec d = (cx*((erand48(Xi)+x)/w-.5) + cy*((erand48(Xi)+y)/h-.5) + cam.d).norm();
        c[i] = c[i] + radiance(Ray(cam.o,d),0,Xi); } }
    if (s >= ns || s >= samps) {
      FILE *f = fopen("image.ppm", "w"); 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/s), toInt(c[i].y/s), toInt(c[i].z/s));
      fclose(f); ns <<= 1; } } }