Awesome

minpt: A path tracer in 300 lines of C++

minpt is a path tracer. The entire code is written with 300 lines of code in C++, without any dependencies except for standard libraries and OpenMP. This project is inspired and motivated by great smallpt path tracer by Kevin Beason, which can generate beautiful image of Cornell Box only with 100 lines of code. The purpose of this project is to answer a question - "what I can do if we are allowed to add a few times more codes than smallpt, with the power of modern standard libraries of C++". This is what I did. I tried to introduce various features which we can often find in various rendering systems, within three times more code than smallpt - 300 lines of code.

Features

• Global illumination renderer with path tracing
• 300 lines, 80 columns of C++ code
• No library dependency except for standard libraries and OpenMP
• Triangle mesh support
• OBJ model support, including MTL files (partial)
• Acceleration structure with SAHBVH
• Parallel rendering support via OpenMP
• Parallel BVH construction
• Diffuse, specular, and glossy BSDFs
• Cook-Torrance model with anisotropic GGX normal distribution
• Texture mapping with alpha textures
• Area light, environment light support
• Next event estimation (NEE)
• Importance sampling for BSDFs and environment light
• Multiple importance sampling among NEE and BSDF sampling
• Pinhole camera / realistic camera

Code

#define _CRT_SECURE_NO_WARNINGS
#include <random>
#include <algorithm>
#include <numeric>
#include <optional>
#include <fstream>
#include <unordered_map>
#include <filesystem>
#include <queue>
#include <thread>
#include <atomic>
#include <condition_variable>
#include <cstring>
#include <omp.h>
#ifdef _MSC_VER
using namespace std; namespace fs=experimental::filesystem; using I=int;
I bswap(I x){return _byteswap_ulong(x);} string pp(string p){return p;}
#elif defined(__GNUC__)
using namespace std; namespace fs=filesystem; using I=int;
I bswap(I x){return __builtin_bswap32(x);}
string pp(string p){replace(p.begin(),p.end(),'\\','/'); return p;}
#endif
template<class T> using op=optional<T>; using F=double;
using C=char; constexpr F Inf=1e+10,Eps=1e-4,Pi=3.14159265358979323846;
struct V{F x,y,z; V(F v=0):V(v,v,v){} V(F x,F y,F z):x(x),y(y),z(z){}
 F operator[](I i)const{return(&x)[i];} F m()const{return std::max({x,y,z});}};
V operator+(V a,V b){return{a.x+b.x,a.y+b.y,a.z+b.z};}
V operator-(V a,V b){return{a.x-b.x,a.y-b.y,a.z-b.z};}
V operator*(V a,V b){return{a.x*b.x,a.y*b.y,a.z*b.z};}
V operator/(V a,V b){return{a.x/b.x,a.y/b.y,a.z/b.z};}
V operator-(V v){return{-v.x,-v.y,-v.z};} F sq(F v){return v*v;}
V vmin(V a,V b){return{min(a.x,b.x),min(a.y,b.y),min(a.z,b.z)};}
V vmax(V a,V b){return{max(a.x,b.x),max(a.y,b.y),max(a.z,b.z)};}
F dot(V a,V b){return a.x*b.x+a.y*b.y+a.z*b.z;}
V norm(V v){return v/sqrt(dot(v,v));} V refl(V w,V n){return 2*dot(w,n)*n-w;}
V intp(V a,V b,V c,F u,F v) {return a*(1-u-v)+b*u+c*v;}
V cross(V a,V b){return{a.y*b.z-a.z*b.y,a.z*b.x-a.x*b.z,a.x*b.y-a.y*b.x};}
tuple<V,V> odr(V n){F s=copysign(1,n.z),a=-1/(s+n.z),b=n.x*n.y*a;
 V u(1+s*n.x*n.x*a,s*b,-s*n.x),v(b,s+n.y*n.y*a,-n.y); return{u,v};}
op<V> refr(V wi,V n,F et){F t=dot(wi,n),t2=1-et*et*(1-t*t);
 return t2>0?et*(n*t-wi)-n*sqrt(t2):op<V>{};}
struct Rng{mt19937 eng; uniform_real_distribution<F> dist;
 Rng(){}; Rng(I seed){eng.seed(seed); dist.reset();}
 F u(){return dist(eng);} V uD(){F r=sqrt(u()),t=2*Pi*u();
  F x=r*cos(t),y=r*sin(t); return V(x,y,sqrt(max(.0,1-x*x-y*y)));}};
struct Dist{vector<F> c{0}; void add(F v){c.push_back(c.back()+v);}
 void norm(){F sum=c.back(); for(F& v:c){v/=sum;}}
 F p(I i)const{return (i<0||i+1>=I(c.size()))?0:c[i+1]-c[i];}
 I samp(Rng& rn)const{auto it=upper_bound(c.begin(),c.end(),rn.u());
  return clamp(I(distance(c.begin(),it))-1,0,I(c.size())-2);}};
struct Dist2{vector<Dist> ds; Dist m; I w,h;
 void init(const vector<F>& v,I a,I b){w=a; h=b; ds.assign(h,{});
  for(I i=0;i<h;i++){auto& d=ds[i]; for(I j=0;j<w;j++){d.add(v[i*w+j]);}
   m.add(d.c.back()); d.norm();} m.norm();}
 F p(F u,F v)const{I y=min(I(v*h),h-1); return m.p(y)*ds[y].p(I(u*w))*w*h;}
 tuple<F,F> samp(Rng& rn)const{I y=m.samp(rn),x=ds[y].samp(rn);
  return{(x+rn.u())/w,(y+rn.u())/h};}}; struct Ind{I p=-1,t=-1,n=-1;};
namespace T{enum{None=0,AreaL=1<<0,EnvL=1<<1,E=1<<2,D=1<<3,G=1<<4,M=1<<5,
 FresRefl=1<<6,FresTran=1<<7,PRefl=1<<8,Refl=D|G|PRefl|FresRefl,L=AreaL|EnvL,
 Tran=FresTran|M,Fres=FresRefl|FresTran,S=PRefl|Fres|M,NonS=D|G};}
struct Geo{vector<V> ps,ns,ts;}; struct R{V o,d;}; struct Surf{V p,n,t,u,v;
 Surf(){} Surf(V p,V n,V t):p(p),n(n),t(t){tie(u,v)=odr(n);}
 bool op(V wi,V wo)const{return dot(wi,n)*dot(wo,n)<=0;}
 tuple<V,V,V> obn(V wi)const{I i=dot(wi,n)>0;return{i?n:-n,u,i?v:-v};}};
struct H{Surf sp; struct Obj* o;}; struct Tex{I w,h; vector<F> cs,as;
 I fl(I i){I j=i/3,x=j%w,y=j/w; return 3*((h-y-1)*w+x)+i%3;}
 F pf(I i,F e,vector<uint8_t>& ct){return pow(F(ct[fl(i)])/e,2.2);}
 F pf(I i,F e,vector<float>& ct){if(e<0){return ct[fl(i)];}
  auto m=bswap(*(int32_t*)&ct[fl(i)]); return *(float*)&m;}
 template <class T> void loadpxm(vector<F>& c,string p){puts(p.c_str());
  static vector<T> ct; FILE* f=fopen(p.c_str(),"rb"); if(!f){return;}
  F e; fscanf(f,"%*s %d %d %lf%*c",&w,&h,&e); I sz=w*h*3,s=sizeof(T);
  ct.assign(sz,0); c.assign(sz,0); fread(ct.data(),s,sz,f);
  for(I i=0;i<sz;i++){c[i]=pf(i,e,ct);} fclose(f);}
 void loadpfm(string p){loadpxm<float>(cs,p);} void load(string p){
  auto b=fs::path(p); auto pc=b.replace_extension(".ppm").string(); 
  auto pa=(b.parent_path()/fs::path(b.stem().string()+"_alpha.ppm")).string();
  loadpxm<uint8_t>(cs,pc); loadpxm<uint8_t>(as,pa);}
 V ev(V t,bool alpha=0)const{F u=t.x-floor(t.x),v=t.y-floor(t.y);
  I x=clamp(I(u*w),0,w-1),y=clamp(I(v*h),0,h-1),i=w*y+x;
  return alpha?V(as[3*i]):V(cs[3*i],cs[3*i+1],cs[3*i+2]);}};
struct B{V mi=V(Inf),ma=V(-Inf); V center()const{return (mi+ma)*.5;}
 F sa()const{V d=ma-mi; return 2*(d.x*d.y+d.y*d.z+d.z*d.x);}
 V operator[](I i)const{return(&mi)[i];}
 bool isect(const R& r,F tl,F th)const{for(I i=0;i<3;i++){F vd=1/r.d[i];
  F t1=(mi[i]-r.o[i])*vd,t2=(ma[i]-r.o[i])*vd;if(vd<0){swap(t1,t2);}
  tl=max(t1,tl); th=min(t2,th); if(th<tl){return 0;}} return 1;};};
B merge(B b,V p){return{vmin(b.mi,p),vmax(b.ma,p)};}
B merge(B a,B b){return{vmin(a.mi,b.mi),vmax(a.ma,b.ma)};}
F gt(const Surf& s1, const Surf& s2){V d=s2.p-s1.p; F L2=dot(d,d);
 d=d/sqrt(L2); return abs(dot(s1.n,d))*abs(dot(s2.n,-d))/L2;};
struct S{I t; R r; V w; F pcs;}; struct SL{V wo; F d; V fs; F p;};
struct Mat{I t=T::None; V Kd,Ks,Ke; Tex* mapKd=0; F Ni,Ns,an,ax,ay;};
struct Lens{F cr,t,e,ar;}; struct Obj{I t; Mat* M=0; vector<Ind> fs;
 struct{V p,u,v,w; F tf,a,fs,id,ss; vector<Lens> ls; vector<op<B>> pbs;} E;
 struct{Dist dist; F invA;} L; struct{Tex map; F rot; Dist2 dist;} EL;
 op<R> trl(R r)const{F z=0; for(I i=I(E.ls.size())-1;i>=0;i--){
  auto& l=E.ls[i]; z-=l.t; struct LH{F t; V n;}; auto h=[&]()->op<LH>{
   if(l.cr==0){F t=(z-r.o.z)/r.d.z; return t<0?op<LH>{}:LH{t,{}};}
   V c(0,0,z+l.cr),oc=c-r.o; F b=dot(oc,r.d),dt=b*b-dot(oc,oc)+l.cr*l.cr;
   if(dt<0){return{};} F t0=b-sqrt(dt),t1=b+sqrt(dt);
   F t=(r.d.z>0)^(l.cr<0)?min(t0,t1):max(t0,t1); if(t<0){return{};}
   V n=(r.o+t*r.d-c)/l.cr; n=dot(n,-r.d)<0?-n:n; return LH{t,n}; }();
  if(!h){return{};} V p=r.o+r.d*h->t; if(p.x*p.x+p.y*p.y>l.ar*l.ar){return{};}
  r.o=p; if(l.cr==0){continue;} F et=l.e/(i>0&&E.ls[i-1].e!=0?E.ls[i-1].e:1);
  auto wt=refr(-r.d,h->n,et); if(!wt){return{};} r.d=*wt; } return r; }
 F ffd(F id)const{op<R> r; auto& lb=E.ls.back(); for(I i=9;i>0;i--){
  if(r=trl({V(0,0,-lb.t+id),norm(V(lb.ar*i/10,0,-id))})){break;}}
  if(!r){return Inf;} F t=-r->o.x/r->d.x,z=(r->o+t*r->d).z;
  F sz=0; for(auto& l:E.ls){sz+=l.t;} return z<sz?-z:Inf;};
 void initE(string p,V e,V c,V u,F fv,F fd,F fs,F ss,F a){puts(p.c_str());
  E.p=e; E.tf=tan(fv*Pi/180*.5); E.a=a; E.ss = ss; E.fs=fs*.001; 
  E.w=norm(e-c); E.u=norm(cross(u,E.w)); E.v=cross(E.w,E.u); C l[4096];
  ifstream f(p); while(f.getline(l,4096)){if(l[0]=='#'||l[0]=='\0'){continue;}
   F cr,t,eta,ar; sscanf(l,"%lf %lf %lf %lf",&cr,&t,&eta,&ar);
   E.ls.push_back({cr*.001,t*.001,eta,ar*.001*.5});} if(E.ls.empty()){return;}
  F lo=Eps,hi=Inf; for(I i=0;i<99;i++){F mi=(lo+hi)*.5; (ffd(mi)<fd?hi:lo)=mi;}
  E.id=hi; Rng rn(42); I n=64; E.pbs.assign(n,{});
  F cfv=-1,sy=sqrt(E.fs*E.fs/(1+E.a*E.a)),sx=E.a*sy; I m=1<<12;
  for(I i=0;i<n;i++){B b; bool f=0; auto& lb=E.ls.back(); V p(0,0,-lb.t+E.id);
   for(I j=0;j<m;j++){p.x=(i+rn.u())/n*E.fs*.5; F r=sqrt(rn.u()),t=2*Pi*rn.u();
    V pl(r*cos(t)*lb.ar,r*sin(t)*lb.ar,-lb.t); auto rt=trl({p,norm(pl-p)});
    if(rt){f=1; b=merge(b,pl); if(p.x<sx*.5){cfv=max(cfv,rt->d.z);}}}
   if(f){E.pbs[i]=b;}} printf("%lf\n",atan(tan(Pi-acos(cfv))/E.a)*180/Pi*2);}
 void initEL(string p,F rot){EL.map.loadpfm(p); EL.rot=rot*Pi/180;
  auto& cs=EL.map.cs; I w=EL.map.w,h=EL.map.h; vector<F> ls(w*h);
  for(I i=0;i<w*h;i++){V v(cs[3*i],cs[3*i+1],cs[3*i+2]);
   ls[i]=v.m()*sin(Pi*(i/w+.5)/h);} EL.dist.init(ls,w,h);}
 void initAreaL(const Geo& G){for(size_t fi=0;fi<fs.size();fi+=3){
   V a=G.ps[fs[fi].p],b=G.ps[fs[fi+1].p],c=G.ps[fs[fi+2].p],cr=cross(b-a,c-a);
   L.dist.add(sqrt(dot(cr,cr))*.5);} L.invA=1/L.dist.c.back(); L.dist.norm();}
 F D(V wh,V u,V v,V n)const{F x=M->ax,y=M->ay;
  return 1/(Pi*x*y*sq(sq(dot(wh,u)/x)+sq(dot(wh,v)/y)+sq(dot(wh,n))));}
 F G(V wi,V wo,V u,V v,V n)const{auto G1=[&](V w){F c=dot(w,n),s=sqrt(1-c*c);
  F cp=dot(w,u)/s,cs=dot(w,v)/s,a2=sq(cp*M->ax)+sq(cs*M->ay);
  return c==0?0:2/(1+sqrt(1+a2*sq(s/c)));}; return G1(wi)*G1(wo);}
 op<S> samp(Rng& rn,I ty,const Surf& sp,const V& wi)const{
  if(ty&T::E){V rp=2*wi-1; I n=I(E.pbs.size()); auto& lb=E.ls.back();
   if(E.ls.empty()){V d=-norm(V(E.a*E.tf*rp.x,E.tf*rp.y,1));
     return S{T::E,E.p,E.u*d.x+E.v*d.y+E.w*d.z,V(1),1};}
   F sy=sqrt(E.fs*E.fs/(1+E.a*E.a)); V o=rp*V(E.a*sy*.5,sy*.5,0); o.z=E.id;
   F l=sqrt(o.x*o.x+o.y*o.y); I i=clamp(I(l/E.fs*2*n),0,n-1); auto& b=E.pbs[i];
   if(!b){return {};} V bl=b->ma-b->mi; V p=b->mi+bl*V(rn.u(),rn.u(),0);
   F s=l!=0?o.y/l:0,c=l!=0?o.x/l:1; V d=norm(V(c*p.x-s*p.y,s*p.x+c*p.y,p.z)-o);
   auto r=trl({o,d}); if(!r){return{};} F A=bl.x*bl.y,Z=lb.t+E.id;
   F w=d.z*d.z*d.z*d.z*A/(Z*Z); o=r->o; d=r->d; o=E.u*o.x+E.v*o.y+E.w*o.z+E.p;
   d=E.u*d.x+E.v*d.y+E.w*d.z; return S{T::E,o,d,V(w*E.ss),1};}
  else if(ty&T::NonS){auto* mp=M->mapKd; V Kd=(mp?mp->ev(sp.t):M->Kd);
   F wd=Kd.m(),ws=M->Ks.m(); if(wd==0&&ws==0){wd=1;ws=0;}
   F s=wd+ws; wd/=s; ws/=s; ty=rn.u()<wd?T::D:T::G;
   if(ty&T::D){auto[n,u,v]=sp.obn(wi); bool usea=mp&&!mp->as.empty();
    F a=usea?mp->ev(sp.t,1).x:1; V d=rn.uD(); return rn.u()>a ?
     S{T::M,sp.p,-wi,V(1),wd}:S{T::D,sp.p,u*d.x+v*d.y+n*d.z,Kd,wd};}
   else if(ty&T::G){auto[n,u,v]=sp.obn(wi); F u1=rn.u()*2*Pi,u2=rn.u();
    V wh=norm(sqrt(u2/(1-u2))*(M->ax*cos(u1)*u+M->ay*sin(u1)*v)+n);
    V wo=refl(wi,wh); if(sp.op(wi,wo)){return{};}
    return S{T::G,sp.p,wo,ev(T::G,sp,wi,wo)/pdf(T::G,sp,wi,wo),ws};}}
  else if(ty&T::PRefl){return S{T::PRefl,sp.p,refl(wi,sp.n),V(1),1};}
  else if(ty&T::Fres){I i=dot(wi,sp.n)>0;V n=i?sp.n:-sp.n;F et=i?1/M->Ni:M->Ni;
   auto wt=refr(wi,n,et); F Fr=!wt?1:[&](){F cos=i?dot(wi,sp.n):dot(*wt,sp.n);
    F r=(1-M->Ni)/(1+M->Ni); r=r*r; return r+(1-r)*pow(1-cos,5);}();
   return rn.u()<Fr?S{T::FresRefl,sp.p,refl(wi,sp.n),V(1),1}
    : S{T::FresTran,sp.p,*wt,V(et*et),1};} return{};}
 F pdf(I type,const Surf& sp,const V& wi,const V& wo)const{
  if(sp.op(wi,wo)){return 0;} if(type&T::D){return 1/Pi;}
  else if(type&T::G){ V wh=norm(wi+wo); auto[n,u,v]=sp.obn(wi);
   return D(wh,u,v,n)*dot(wh,n)/(4*dot(wo,wh)*dot(wo,n));} return 0;}
 op<SL> sampL(Rng& rn,const Geo& G,const Surf& sp)const{
  if(t&T::AreaL){I i=L.dist.samp(rn); F s=sqrt(max(.0,rn.u()));
   V a=G.ps[fs[3*i].p],b=G.ps[fs[3*i+1].p],c=G.ps[fs[3*i+2].p];
   Surf spL(intp(a,b,c,1-s,rn.u()*s),norm(cross(b-a,c-a)),{});
   V pp=spL.p-sp.p,wo=norm(pp); F p=pdfL(sp,spL,-wo);
   return p==0?op<SL>{}:SL{wo,sqrt(dot(pp,pp)),ev(T::AreaL,spL,{},-wo),p};}
  else if(t&&T::EnvL){auto[u,v]=EL.dist.samp(rn); F t=Pi*v,st=sin(t);
   F p=2*Pi*u+EL.rot; V wo(st*sin(p),cos(t),st*cos(p)); F pL=pdfL(sp,{},-wo);
   return pL==0?op<SL>{}:SL{wo,Inf,ev(T::EnvL,{},{},-wo),pL};} return{};}
 F pdfL(const Surf& sp, const Surf& spL, const V& wo)const{
  if(t&T::AreaL){F G=gt(sp,spL); return G==0?0:L.invA/G;}
  else if(t&&T::EnvL){V d=-wo; F at=atan2(d.x,d.z); at=at<0?at+2*Pi:at;
   F t=(at-EL.rot)*.5/Pi; F u=t-floor(t),v=acos(d.y)/Pi; F st=sqrt(1-d.y*d.y);
   return st==0?0:EL.dist.p(u,v)/(2*Pi*Pi*st*abs(dot(-wo,sp.n)));} return 0;}
 V ev(I ty,const Surf& sp,const V& wi,const V& wo)const{
  if(ty&T::AreaL){return dot(wo,sp.n)<=0 ? V() : M->Ke;}
  else if(ty&T::EnvL){V d=-wo; F at=atan2(d.x,d.z); at=at<0?at+2*Pi:at;
   F t=(at-EL.rot)*.5/Pi; return EL.map.ev({t-floor(t),acos(d.y)/Pi,0});}
  else if(ty&T::G){if(sp.op(wi,wo)){return{};} V wh=norm(wi+wo);
   auto[n,u,v]=sp.obn(wi); V Fr=M->Ks+(1-M->Ks)*pow(1-dot(wo,wh),5);
   return M->Ks*Fr*(D(wh,u,v,n)*G(wi,wo,u,v,n)/(4*dot(wi,n)*dot(wo,n)));}
  else if(ty&T::D){if(sp.op(wi,wo)){return{};}
   auto* mp=M->mapKd; F a=(mp&&!mp->as.empty())?mp->ev(sp.t,1).x:1;
   return (mp?mp->ev(sp.t):M->Kd)*(a/Pi);} return{};}};
struct Scene{Geo G; Obj E; vector<Obj> os; vector<unique_ptr<Tex>> ts;
 vector<Mat> ms; vector<I> Ls; unordered_map<string,I> msmap,tsmap;
 bool useEL=0; Obj* envL(){return useEL?&os.front():nullptr;}
 tuple<const Obj&,F> sampL(Rng& rn)const{I n=I(Ls.size());
  I i=clamp(I(rn.u()*n),0,n-1); return{os.at(Ls[i]),1./n};}
 bool ws(C c){return c==' '||c=='\t';}; F pdfL(){return 1./Ls.size();}
 bool cm(C*& t,const C* c,I n){return !strncmp(t,c,n)&&ws(t[n]);}
 void ss(C*& t){t+=strspn(t," \t");} void sc(C*& t){t+=strcspn(t,"/ \t");}
 F nf(C*& t){ss(t); F v=atof(t); sc(t); return v;}
 V nv(C*& t){V v; v.x=nf(t); v.y=nf(t); v.z=nf(t); return v;}
 I pi(I i,I vn){return i<0?vn+i:i>0?i-1:-1;} Ind pind(C*& t){Ind i; ss(t);
  i.p=pi(atoi(t),I(G.ps.size())); sc(t); if(t++[0]!='/'){return i;}
  i.t=pi(atoi(t),I(G.ts.size())); sc(t); if(t++[0]!='/'){return i;}
  i.n=pi(atoi(t),I(G.ns.size())); sc(t); return i;}
 void nn(C*& t,C name[]){sscanf(t,"%s",name);};
 void loadmtl(string p){ifstream f(p); C l[4096],name[256]; puts(p.c_str());
  while(f.getline(l,4096)){auto* t=l; ss(t); if(cm(t,"newmtl",6)){
    nn(t+=7,name); msmap[name]=I(ms.size()); ms.emplace_back(); continue;}
   if(ms.empty()){continue;} Mat& m=ms.back(); if(cm(t,"Kd",2)){m.Kd=nv(t+=3);}
   else if(cm(t,"Ks",2)) m.Ks=nv(t+=3); else if(cm(t,"Ni",2)) m.Ni=nf(t+=3);
   else if(cm(t,"Ns",2)) m.Ns=nf(t+=3); else if(cm(t,"aniso",5)) m.an=nf(t+=5);
   else if(cm(t,"Ke",2)){m.Ke=nv(t+=3); m.t|=m.Ke.m()>0?T::L:0;}
   else if(cm(t,"illum",5)){ss(t+=6); I v=atoi(t);
    m.t|=v==7?T::Fres:v==5?T::PRefl:T::NonS;}
   else if(cm(t,"map_Kd",6)){nn(t+=7,name); auto it=tsmap.find(name);
    if(it!=tsmap.end()){m.mapKd=ts[it->second].get(); continue;}
    tsmap[name]=I(ts.size()); ts.emplace_back(new Tex);
    ts.back()->load(pp((fs::path(p).remove_filename()/name).string()));
    m.mapKd=ts.back().get();} else{continue;}}}
 void load(string p,string env,F rot){C l[4096],name[256]; ifstream f(p); 
  if(!env.empty()){useEL=1; os.push_back({T::EnvL}); os.back().initEL(env,rot);
  Ls.push_back(0);} while(f.getline(l,4096)){C* t=l; ss(t); I on=I(os.size());
   if     (cm(t,"v" ,1)){G.ps.emplace_back(nv(t+=2));}
   else if(cm(t,"vn",2)){G.ns.emplace_back(nv(t+=3));}
   else if(cm(t,"vt",2)){G.ts.emplace_back(nv(t+=3));}
   else if(cm(t,"f" ,1)){t+=2; if(ms.empty()){ms.push_back({T::D,1});
    os.push_back({T::D,&ms.back()});} Ind is[4]; for(auto& i:is) i=pind(t);
    auto& fs=os.back().fs; fs.insert(fs.end(),{is[0],is[1],is[2]});  
    if(is[3].p!=-1){fs.insert(fs.end(),{is[0],is[2],is[3]});}}
   else if(cm(t,"usemtl",6)){t+=7; nn(t,name); auto& M=ms[msmap.at(name)];
    os.push_back({!M.t?(M.t=T::NonS):M.t,&M}); if(M.t&T::L) Ls.push_back(on);}
   else if(cm(t,"mtllib",6)){nn(t+=7,name);
    loadmtl(pp((fs::path(p).remove_filename()/name).string()));}else continue;}
  for(auto& o:os) if(o.M&&o.M->t&T::AreaL){o.initAreaL(G);}
  for(auto& m:ms){if(!(m.t&T::NonS)){continue;} F r=2/(2+m.Ns);
   F as=sqrt(1-m.an*.9); m.ax=max(1e-3,r/as); m.ay=max(1e-3,r*as);}}
 struct TH{F t,u,v;}; struct Tri{V p1,e1,e2,n; B b; V c; I oi,fi;
  Tri(const V& p1,const V& p2,const V& p3,I oi,I fi)
   : p1(p1),oi(oi),fi(fi){e1=p2-p1; e2=p3-p1; n=norm(cross(p2-p1,p3-p1));
   b=merge(b,p1); b=merge(b,p2); b=merge(b,p3); c=b.center();}
  op<TH> isect(const R& r,F tl,F th)const{V p=cross(r.d,e2),tv=r.o-p1;
   V q=cross(tv,e1); F d=dot(e1,p),ad=abs(d),s=copysign(1,d);
   F u=dot(tv,p)*s,v=dot(r.d,q)*s; if(ad<1e-8||u<0||v<0||u+v>ad) return{};
   F t=dot(e2,q)/d; return t<tl||th<t ? op<TH>{} : TH{t,u/ad,v/ad};}};
 struct N{B b; bool leaf=0; I s,e,c1,c2;}; vector<N> ns;
 vector<Tri> trs; vector<I> ti; void build(){queue<tuple<I,I,I>> q;
  for(size_t oi=0;oi<os.size();oi++){auto& fs=os[oi].fs;
   for(size_t fi=0;fi<fs.size();fi+=3){V a=G.ps[fs[fi].p],b=G.ps[fs[fi+1].p];
    V c=G.ps[fs[fi+2].p]; trs.emplace_back(a,b,c,I(oi),I(fi));}}
  I nt=I(trs.size()); q.push({0,0,nt}); ns.assign(2*nt-1,{}); ti.assign(nt,0);
  iota(ti.begin(),ti.end(),0); mutex mu; condition_variable cv;
  atomic<I> pr=0; I nn=1; bool done=0; using ul=unique_lock<mutex>;
  auto process=[&](){while(!done){auto [ni,s,e]=[&]()->tuple<I,I,I>{ul lk(mu);
    if(!done&&q.empty()){cv.wait(lk,[&](){return done||!q.empty();});}
    if(done) return{}; auto v=q.front(); q.pop(); return v;}();
   if(done) break; N& n=ns[ni]; for(I i=s;i<e;i++) n.b=merge(n.b,trs[ti[i]].b);
   auto st=[&,s=s,e=e](I ax){sort(&ti[s],&ti[e-1]+1,[&](I i1, I i2){
    return trs[i1].c[ax]<trs[i2].c[ax];});}; F b=Inf; I bi,ba;
   auto lf=[&,s=s,e=e](){n.leaf=1; n.s=s; n.e=e; pr+=e-s;if(pr==I(trs.size())){
    ul lk(mu);done=1;cv.notify_all();}}; if(e-s<2){lf();continue;}
   for(I a=0;a<3;a++){thread_local vector<F> l(nt+1),r(nt+1); st(a); B bl,br;
    for(I i=0;i<=e-s;i++){I j=e-s-i; l[i]=bl.sa()*i; r[j]=br.sa()*i;
    bl=i<e-s?merge(bl,trs[ti[s+i]].b):bl;br=j>0?merge(br,trs[ti[s+j-1]].b):br;}
    for(I i=1;i<e-s;i++){F c=1+(l[i]+r[i])/n.b.sa();if(c<b){b=c;bi=i;ba=a;}}}
   if(b>e-s){lf(); continue;} st(ba); I m=s+bi; unique_lock<mutex> lk(mu);
   q.push({n.c1=nn++,s,m}); q.push({n.c2=nn++,m,e}); cv.notify_one();}};
  vector<thread> ths(omp_get_max_threads());
  for(auto& th:ths){th=thread(process);} for(auto& th:ths){th.join();}}
 op<H> isect(const R& r,F tl=Eps,F th=Inf){op<TH> mh,h; I mi,s[99]{},si=0;
  while(si>=0){auto& n=ns.at(s[si--]); if(!n.b.isect(r,tl,th)){continue;}
   if(!n.leaf){s[++si]=n.c1; s[++si]=n.c2; continue;}
   for(I i=n.s;i<n.e;i++) if(h=trs[ti[i]].isect(r,tl,th)){mh=h;th=h->t;mi=i;}}
  if(!mh){return{};} Tri& tr=trs[ti[mi]]; Obj& o=os[tr.oi]; V p=r.o+th*r.d;
  F u=mh->u,v=mh->v; Ind& i=o.fs[tr.fi],j=o.fs[tr.fi+1],k=o.fs[tr.fi+2];
  V n=i.n<0?tr.n:norm(intp(G.ns[i.n],G.ns[j.n],G.ns[k.n],u,v));
  V t=i.t<0?0:intp(G.ts[i.t],G.ts[j.t],G.ts[k.t],u,v); return H{{p,n,t},&o};}};
I main(I,C**v){Scene sc; puts(v[1]); sc.load(v[1],v[2],atof(v[7])); sc.build();
 I ns=atoi(v[5]),mb=atoi(v[6]),w=atoi(v[8]),h=atoi(v[9]),sz=w*h*3;
 vector<V> D(w*h,V()); sc.E.initE(v[3],V(atof(v[10]),atof(v[11]),atof(v[12])),
  V(atof(v[13]),atof(v[14]),atof(v[15])),V(0,1,0),F(atof(v[16])),
  F(atof(v[17])),F(atof(v[18])),F(atof(v[19])),F(w)/h);
 auto L=[&](Rng& rn,V&& uv){V L,T(1),wi=uv; I ty=T::E; Obj* o=&sc.E; Surf sp;
  for(I b=0;b<mb;b++){bool nee=!sc.Ls.empty()&&b>0&&ty&T::NonS;
   auto s=o->samp(rn,ty,sp,wi); if(!s){break;} auto[t,r,w,pcs]=*s; T=T/pcs;
   if(nee){auto[oL,pLs]=sc.sampL(rn); auto sL=oL.sampL(rn,sc.G,sp);
    if(sL){auto[wo,d,fL,pL]=*sL; if(!sc.isect({sp.p,wo},Eps,d*(1-Eps))){
     L=L+T*o->ev(t,sp,wi,wo)*fL/(o->pdf(t,sp,wi,wo)+pL*pLs);}}}
   T=T*w; auto hi=sc.isect(r); if(!hi){hi=H{{},sc.envL()}; if(!hi->o){break;}}
   auto[spH,oH]=*hi; if(oH->t&T::L){ L=L+T*oH->ev(oH->t&T::L,spH,{},-r.d)/
    (b==0||t&T::S?1:oH->pdfL(sp,spH,-r.d)*sc.pdfL()/o->pdf(t,sp,wi,r.d)+1);}
   if(b>3){F q=max(.2,1-T.m()); T=T/(1-q); if(rn.u()<q) break;}
   wi=-r.d; o=oH; sp=spH; ty=o->t&~T::L;} return L;};
 #pragma omp parallel for schedule(dynamic,1)
 for(I i=0;i<w*h;i++){thread_local Rng rn(42+omp_get_thread_num());
  for(I j=0;j<ns;j++) D[i]=D[i]+L(rn,V((i%w+rn.u())/w,(i/w+rn.u())/h,0))/ns;}
 FILE* f=fopen(v[4],"wb"); fprintf(f, "PF\n%d %d\n-1\n",w,h);
 vector<float> d(sz); for(I y=0;y<h;y++) for(I x=0;x<w;x++) for(I i=0;i<3;i++){
  d[3*(y*w+x)+i]=float(D[(h-1-y)*w+(w-1-x)][i]);}
 fwrite(d.data(),4,d.size(),f); fclose(f); return 0;}

Build

Windows

Tested with Visual Studio 2017 Community (15.6).

> "C:\Program Files (x86)\Microsoft Visual Studio\2017\Community\VC\Auxiliary\Build\vcvarsall.bat" x64
> cl /EHsc /std:c++latest /openmp minpt.cpp

Linux

Tested with GCC 8.1.

$ g++ -std=c++17 -fopenmp -O3 -w minpt.cpp -lstdc++fs -o minpt

Usage

All command line arguments are positional and required.

$ ./minpt \
    scene_path envmap_path lensfile_path output_path \
    spp max_path_length envmap_rotation image_width image_height \
    camera_pos_x camera_pos_y camera_pos_z \
    camera_lookat_pos_x camera_lookat_pos_y camera_lookat_pos_z \
    vertical_fov focus_distance sensor_diagonal_length sensor_sensitivity

• General arguments
  • scene_path: Input path to .obj file
  • envmap_path: Input path to .pfm file for environment map ("": black background)
  • lensfile_path: Input path to lens file ("": use pinhole camera)
  • output_path: Output path to rendered image
  • spp: Number of samples per pixel
  • max_path_length: Maximum path length (=number of bounces+1)
  • envmap_rotation: Rotation angle of the environment map around up vector
  • image_width image_height: Output image resolution
• Camera related arguments
  • camera_pos_{x,y,z}: Camera position
  • camera_lookat_pos_{x,y,z}: Look-at position of the camera
    • Up vector is fixed to (0,1,0)
  • Pinhole camera specific arguments
    • vertical_fov: Vertical field of view
  • Realistic camera specific arguments
    • focus_distance: Focus distance on the optical axis [m]
    • sensor_diagonal_length: Diagonal length of the sensor [mm]
    • sensor_sensitivity: Sensitivity multiplier of the sensor

Input/Output Description

Image file

Minpt only supports binary .ppm format (P6) as input textures, and .pfm format as environment maps. The rendered images are also generated with .pfm format. To check the rendered images, you want to use tools capable of viewing .pfm files, such as PCG HDRITools or HDRView.

OBJ file (.obj)

Triangle and quad meshes are supported.

OBJ material file (.mtl)

Minpt uses illum parameter to select the material type:

• illum=7: Fresnel reflection, refraction BSDF
• illum=5: Perfect mirror BRDF
• otherwise: Mixtured Lambertian and Cook-Torrance BRDF

Anisotropy of GGX normal distribution can be controlled via aniso parameter.

Lens description file

As an input to the realistic camera, minpt requires a lens description file, describing the configuration of the lens system to be simulated. We used the same format referred in Fig.1 of [Kolb et al. 1995]. The file is compatible with lens description file of pbrt-v3. For instance, the following data is wide.22mm.dat in pbrt-v3 scene repository.

35.98738  1.21638 1.54  23.716
11.69718  9.9957  1     17.996
13.08714  5.12622 1.772 12.364
-22.63294 1.76924 1.617 9.812
71.05802  0.8184  1     9.152
0         2.27766 0     8.756
-9.58584  2.43254 1.617 8.184
-11.28864 0.11506 1     9.152
-166.7765 3.09606 1.713 10.648
-7.5911   1.32682 1.805 11.44
-16.7662  3.98068 1     12.276
-7.70286  1.21638 1.617 13.42
-11.97328 0       1     17.996

Examples

Example scenes can be obtained from Morgan McGuire's Computer Graphics Archive. Slight modification of .mtl files and conversion of image files to .ppm are needed to obtain the rendered images.

CornellBox-Sphere

$ ./minpt \
    ./CornellBox/CornellBox-Sphere.obj "" "" result.pfm \
    100 20 0 1920 1080 \
    0 1 5 0 1 0 30 0 0 0

fireplace_room (Teaser 1)

$ ./minpt \
    ./fireplace_room/fireplace_room.obj "" ./wide.22mm.dat \
    result.pfm \
    1000 20 0 1920 1080 \
    5.101118 1.083746 -2.756308 4.167568 1.078925 -2.397892 \
    43.001194 3 35 10

Modification to the .mtl file:

$ diff fireplace_room.mtl fireplace_room_orig.mtl
19c19
< Ns 10
---
> Ns 5
32c32
< Ks 0.2 0.2 0.2
---
> Ks 0.1 0.1 0.1
41,42c41,42
< Ks 0.4 0.4 0.4
< Ns 50
---
> Ks 0.02 0.02 0.02
> Ns 10.00
51c51
< Ni 1.5
---
> Ni 1.00
70c70
< illum 5
---
> illum 3
123,124c123,124
< Ks 0.2 0.2 0.2
< Ns 30
---
> Ks 0.04 0.04 0.04
> Ns 0.06

rungholt (Teaser 2)

The environment map is converted from flower_road_4k.hdr in HDRI Haven. This scene requires approximately 4GB of free memory.

$ ./minpt \
    ./rungholt/rungholt.obj ./flower_road_4k.pfm ./wide.22mm.dat \
    result.pfm \
    1000 20 0 1920 1080 \
    243.523193 155.362595 172.186203 242.806686 154.837311 171.727173 \
    43.273157 100 35 10

Tips

Configuring camera parameters

The following snippet can be used to extract (a part of) camera parameters using Blender. When the realistic camera is enabled, minpt prints effective vertical FoV which can be used to match the image with the image using pinhole camera.

import bpy
import mathutils
from math import pi
import os 
scene = bpy.context.scene
camera = bpy.data.objects["Camera"]
m = camera.matrix_world.copy()
m = mathutils.Matrix.Rotation(-pi/2,4,'X')*m
m.transpose()
p = mathutils.Vector(m[3][0:3])
e = p-mathutils.Vector(m[2][0:3])
fov = bpy.data.cameras["Camera"].angle_y*180/pi
print("%d %d %f %f %f %f %f %f %f" % (
  scene.render.resolution_x,scene.render.resolution_y,
  p[0],p[1],p[2],e[0],e[1],e[2],fov))

Details

• Ray-triangle intersection
  • Möller–Trumbore algorithm [Möller & Trumbore 1997]
• Ray-AABB intersection
  • http://psgraphics.blogspot.de/2016/02/new-simple-ray-box-test-from-andrew.html
• Orthonormal basis computation
  • Method by Duff et al. [2017]
  • Branchless version using copysign
• Parallel BVH build
  • Simple producer-consumer design
  • Blocking queue
• Environment light
  • Image-based importance sampling [Colbert et al. 2010]
• Realistic camera [Kolb et al. 1995] [Steinert et al. 2011]
  • Importance sampling of exit pupils
  • Autofocus using binary search

References

• [Möller & Trumbore 1997] Fast, Minimum Storage Ray-Triangle Intersection. JGT. 1997.
• [Duff et al. 2017] Building An Orthonormal Basis, Revisited. JCGT. 2017.
• [Colbert et al. 2010] Importance Sampling for Production Rendering. SIGGRAPH Course. 2010.
• [Kolb et al. 1995] A Realistic Camera Model for Computer Graphics. SIGGRAPH. 1995.
• [Steinert et al. 2011] General Spectral Camera Lens Simulation. CGF. 2011.
• [Wald 2007] On fast Construction of SAH based Bounding Volume Hierarchies. Proc. Eurographics/IEEE Symposium on Interactive Ray Tracing. 2007.

Author

Hisanari Otsu (Personal site, Twitter)

License

This software is licensed under MIT license. See LICENSE file for detail.