fork download
copy 
  1. #include <math.h> 
  2. #include <stdlib.h> 
  3. #include <stdio.h> 
  4. struct Vec{double x, y, z;Vec(double x_=0,double y_=0,double z_=0){x=x_;y=
  5. y_; z=z_; } Vec operator+(const Vec &b) const { return Vec(x+b.x,y+b.y,
  6. z+b.z); } Vec operator-(const Vec &b) const { return Vec(x-b.x,y-b.y,z-
  7. b.z); } Vec operator*(double b) const { return Vec(x*b,y*b,z*b); } Vec 
  8. mult(const Vec &b) const { return Vec(x*b.x,y*b.y,z*b.z); } Vec& norm(){ 
  9. return *this = *this * (1/sqrt(x*x+y*y+z*z)); } double dot(const Vec &
  10. b) const { return x*b.x+y*b.y+z*b.z; }  Vec operator%(Vec&b){return 
  11. 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(
  12. Vec o_, Vec d_) : o(o_), d(d_) {} };enum Refl_t { DIFF, SPEC, REFR }; 
  13. struct Sphere { double rad;  Vec p, e, c;  Refl_t refl;  Sphere(double 
  14. rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_): rad(rad_), p(p_), e(e_), 
  15. c(c_), refl(refl_) {} double intersect(const Ray &r) const {  Vec op = 
  16. p-r.o;  double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad; 
  17. if (det<0) return 0; else det=sqrt(det); return (t=b-det)>eps ? t : ((
  18. t=b+det)>eps ? t : 0); }};Sphere spheres[] = { Sphere(1e5, Vec( 1e5+1,40.8,81.6), 
  19. Vec(),Vec(.75,.25,.25),DIFF), Sphere(1e5, Vec(-1e5+99,40.8,81.6),Vec(),
  20. Vec(.25,.25,.75),DIFF), Sphere(1e5, Vec(50,40.8, 1e5), Vec(),Vec(.75,.75,.75),
  21. DIFF), Sphere(1e5, Vec(50,40.8,-1e5+170), Vec(),Vec(), DIFF), Sphere(1
  22. e5, Vec(50, 1e5, 81.6), Vec(),Vec(.75,.75,.75),DIFF), Sphere(1e5, Vec(50,-1
  23. e5+81.6,81.6),Vec(),Vec(.75,.75,.75),DIFF), Sphere(16.5,Vec(27,16.5,47), 
  24. Vec(),Vec(1,1,1)*.999, SPEC), Sphere(16.5,Vec(73,16.5,78), Vec(),Vec(1,1,1)*.999, 
  25. REFR), Sphere(600, Vec(50,681.6-.27,81.6),Vec(12,12,12), Vec(), DIFF) };
  26. inline double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; }inline 
  27. int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); }inline 
  28. bool intersect(const Ray &r, double &t, int &id){ double n=sizeof(
  29. spheres)/sizeof(Sphere), d, inf=t=1e20; for(int i=int(n);i--;) if((d=
  30. spheres[i].intersect(r))&&d<t){t=d;id=i;} return t<inf;}Vec radiance(
  31. const Ray &r, int depth, unsigned short *Xi){ double t;  int id=0;  if (!
  32. intersect(r, t, id)) return Vec();  const Sphere &obj = spheres[id];  
  33. Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c; 
  34. double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z;  if (++
  35. depth>5) if (erand48(Xi)<p) f=f*(1/p); else return obj.e;  if (obj.
  36. refl == DIFF){  double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(
  37. r2); Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u; Vec 
  38. d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm(); return obj.
  39. e + f.mult(radiance(Ray(x,d),depth,Xi)); } else if (obj.refl == SPEC)  
  40. return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi)); 
  41. Ray reflRay(x, r.d-n*2*n.dot(r.d));  bool into = n.dot(nl)>0;  double 
  42. nc=1, nt=1.5, nnt=into?/nt:/nc, ddn=r.d.dot(nl), cos2t; if ((cos2t=1-
  43. nnt*nnt*(1-ddn*ddn))<0)  return obj.e + f.mult(radiance(reflRay,depth,
  44. Xi)); Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).
  45. norm(); double a=nt-nc, b=nt+nc, R0=a*/(b*b), c = 1-(into?-ddn:tdir.
  46. dot(n)); double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=/P,TP=/(1-
  47. P); return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ?  radiance(reflRay,
  48. depth,Xi)*RP:radiance(Ray(x,tdir),depth,Xi)*TP) : radiance(reflRay,
  49. depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr);}int main(int argc, 
  50. char *argv[]){ int w=1024, h=768, samps = argc==2 ? atoi(argv[1])/4 : 1;  
  51. Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm());  Vec cx=Vec(w*.5135/
  52. h), cy=(cx%cam.d).norm()*.5135, r, *c=new Vec[w*h];
  53. #pragma omp parallel for schedule(dynamic, 1) private(r)  
  54. for (int y=0; y<h; y++){fprintf(stderr,"\rRendering(%d spp) %5.2f%%",
  55. samps*4,100.*y/(h-1));for (unsigned short x=0, Xi[3]={0,0,y*y*y}; x<w;
  56. x++)for (int sy=0,i=(h-y-1)*w+x; sy<2; sy++)  for (int sx=0; sx<2; sx++, 
  57. r=Vec()){for(int s=0; s<samps; s++){double r1=2*erand48(Xi),dx=r1<1?sqrt
  58. (r1)-1:1-sqrt(2-r1); double r2=2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt
  59. (2-r2);Vec d=cx*(((sx+.5+dx)/2+x)/w-.5)+cy*(((sy+.5+dy)/2+y)/
  60. h-.5)+cam.d;r=r+radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./
  61. samps);}c[i]=c[i]+Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25;}} 
  62. FILE *f=fopen("image.ppm","w");fprintf(f,"P3\n%d %d\n%d\n",w,h,255); 
  63. for(int i=0;i<w*h;i++)fprintf(f,"%d %d %d ",toInt(c[i].x),toInt(
  64. c[i].y),toInt(c[i].z));}
Not running #stdin #stdout 0s 0KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/fYDPfp
language:
C++14 (gcc 8.3)
created:
12 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!