smallpt on C language

smallpt.c

// Original: https://www.kevinbeason.com/smallpt/
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
typedef struct { double x, y, z; } Vec;    // position, also color (r,g,b)
Vec vadd(Vec a, Vec b) { return (Vec){a.x+b.x, a.y+b.y, a.z+b.z}; }
Vec vsub(Vec a, Vec b) { return (Vec){a.x-b.x, a.y-b.y, a.z-b.z}; }
Vec vscale(Vec a, double b) { return (Vec){a.x*b, a.y*b, a.z*b}; }
Vec vmult(Vec a, Vec b) { return (Vec){a.x*b.x, a.y*b.y, a.z*b.z}; }
double vdot(Vec a, Vec b) { return a.x*b.x+a.y*b.y+a.z*b.z; }
double vlen(Vec a) { return sqrt(vdot(a, a)); }
Vec vnorm(Vec a) { double s=1/vlen(a); return (Vec){a.x*s, a.y*s, a.z*s}; }
Vec vcross(Vec a, Vec b) { return (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};}
typedef struct { Vec o, d; } Ray;
typedef enum { DIFF, SPEC, REFR } Refl_t;  // material types, used in radiance()
typedef struct { double rad; Vec p, e, c; Refl_t refl; } Sphere;
double sphere_intersect(const Sphere *s, const Ray *r) { // returns distance, 0 if nohit
  Vec op = vsub(s->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=vdot(op, r->d), det=b*b-vdot(op, op)+s->rad*s->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},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){.999,.999,.999},SPEC},//Mirr
  (Sphere){16.5,(Vec){73,16.5,78},    (Vec){},(Vec){.999,.999,.999},REFR},//Glas
  (Sphere){600, (Vec){50,681.6-.27,81.6},(Vec){12,12,12},  (Vec){}, DIFF},//Lite
};
double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; }
int toInt(double x){ return (int)(pow(clamp(x),1/2.2)*255+.5); }
int 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=sphere_intersect(&spheres[i], r))&&d<*t){*t=d;*id=i;}
  return *t<inf;
}
Vec radiance(const Ray *r, int depth, unsigned short *Xi){
  double t; int id=0;             // distance to intersection, 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=vadd(r->o,vscale(r->d,t)), n=vnorm(vsub(x,obj->p)), nl=vdot(n,r->d)<0?n:vscale(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=vscale(f,1/p); else return obj->e;} //R.R.
  Vec rrr;
  if (obj->refl == DIFF){                  // Ideal DIFFUSE reflection
    double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
    Vec w=nl, u=vnorm(vcross(fabs(w.x)>.1?(Vec){0,1}:(Vec){1},w)), v=vcross(w,u);
    Vec d=vnorm(vadd(vadd(vscale(u,cos(r1)*r2s), vscale(v,sin(r1)*r2s)), vscale(w,sqrt(1-r2))));
    rrr = radiance(&(Ray){x,d},depth,Xi);
  } else {  // Mirror or glass
    Ray reflRay={x, vsub(r->d, vscale(n,2*vdot(n,r->d)))};     // Ideal dielectric REFRACTION
    int into = vdot(n,nl)>0;              // Ray from outside going in?
    double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=vdot(r->d,nl), cos2t;
    if (obj->refl == SPEC || (cos2t=1-nnt*nnt*(1-ddn*ddn))<0) { // Ideal SPECULAR reflection, or Total internal reflection
      rrr = radiance(&reflRay,depth,Xi);
    } else {
      Vec tdir = vnorm(vsub(vscale(r->d,nnt), vscale(n, (into?1:-1)*(ddn*nnt+sqrt(cos2t)))));
      double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:vdot(tdir,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);
      if (depth>2) rrr = (erand48(Xi)<P) ? vscale(radiance(&reflRay,depth,Xi),RP) : vscale(radiance(&(Ray){x,tdir},depth,Xi),TP);
      else         rrr = vadd(vscale(radiance(&reflRay,depth,Xi),Re), vscale(radiance(&(Ray){x,tdir},depth,Xi),Tr));
    }
  }
  return vadd(obj->e, vmult(f, rrr));
}
int main(int argc, char *argv[]){
  int w=1024, h=768, samps = argc==2 ? atoi(argv[1])/4 : 1; // # samples
  Ray cam={(Vec){50,52,295.6}, vnorm((Vec){0,-0.042612,-1})}; // cam pos, dir
  Vec cx={w*.5135/h}, cy=vscale(vnorm(vcross(cx,cam.d)), .5135);
  Vec *c=calloc(w*h, sizeof(Vec));
#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));
    for (unsigned short x=0, Xi[3]={0,0,y*y*y}; 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++){ Vec 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 = vadd(vadd(vscale(cx, ((sx+.5 + dx)/2 + x)/w - .5),
                              vscale(cy, ((sy+.5 + dy)/2 + y)/h - .5)), cam.d);
            d = vnorm(d);
            r = vadd(r, vscale(radiance(&(Ray){vadd(cam.o, vscale(d, 140)),d},0,Xi), 1./samps));
          } // Camera rays are pushed ^^^^^ forward to start in interior
          c[i] = vadd(c[i], vscale((Vec){clamp(r.x),clamp(r.y),clamp(r.z)}, .25));
        }
  }
  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));
}