API tools faq
paste
Guest User

Eric Graham's Juggler RT1.C

a guest
Apr 9th, 2015
766
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 11.13 KB | None | 0 0
raw download clone embed print report
  1. /*      RT1.C    Ray tracing program
  2.         Copyright 1987 Eric Graham
  3.         Permission is granted to copy and modify this file, provided
  4.         that this notice is retained.
  5. */
  6.  
  7.  
  8.  
  9.  
  10. #define BIG 1.0e10
  11. #define SMALL 1.0e-3
  12. #define DULL    0
  13. #define BRIGHT  1
  14. #define MIRROR  2
  15.  
  16. double dot();   /* Vector dot product */
  17. struct lamp {
  18.     double pos[3];      /* position of lamp */
  19.     double color[3];    /* color of lamp */
  20.     double radius;      /* size of lamp */
  21. };
  22.  
  23. struct sphere {
  24.     double pos[3];      /* position of sphere */
  25.     double color[3];    /* color of sphere */
  26.     double radius;      /* size of sphere */
  27.     int type;   /* type of surface, DULL, BRIGHT or MIRROR */
  28. };
  29.  
  30. struct patch {          /* a small bit of something visible */
  31.     double pos[3];      /* position */
  32.     double normal[3];   /* direction 90 degrees to surface */
  33.     double color[3];    /* color of patch */
  34. };
  35.  
  36. struct world {  /* everything in the universe, except observer */
  37.     int numsp;          /* number of spheres */
  38.     struct sphere *sp;  /* array of spheres */
  39.     int numlmp;         /* number of lamps */
  40.     struct lamp *lmp;   /* array of lamps */
  41.     struct patch horizon[2]; /* alternate squares on the ground */
  42.     double illum[3];    /* background diffuse illumination */
  43.     double skyhor[3];   /* sky color at horizon */
  44.     double skyzen[3];   /* sky color overhead */
  45. };
  46.  
  47. struct observer {       /* now the observer */
  48.     double obspos[3];   /* his position */
  49.     double viewdir[3];  /* direction he is looking */
  50.     double uhat[3];     /* left to right in view plane */
  51.     double vhat[3];     /* down to up in view plane */
  52.     double fl,px,py;    /* focal length and pixel sizes */
  53.     int nx,ny;          /* number of pixels */
  54. };
  55.  
  56. main()
  57. {
  58.     double line[6],brite[3];
  59.     struct observer o;  struct world w;
  60.     int i,j,ii,jj,skip; short int si,sj;
  61.  
  62.     setup(&o,&w,&skip); /*  Provide this function to set up the
  63.                             observer and the world */
  64.     si=1+(o.nx-1)/skip;
  65.     sj=1+(o.ny-1)/skip;
  66.     initsc(si,sj);      /*  Set up the screen for Hold and Modify
  67.                             mode. See the ROM Kernel manual */
  68.     for (jj=j=0; j<o.ny; j+=skip,jj++) {
  69.         for (ii=i=0; i<o.nx; i+=skip,ii++){
  70.             pixline(line,&o,i,j); raytrace(brite,line,&w);
  71.             ham(ii,jj,brite);     /* Provide this function to */
  72.         }                         /* a pixel */
  73.     }
  74.     cleanup(0);         /* Free up resources allocated */
  75. }                       /* by initsc() */
  76.  
  77. raytrace(brite,line,w)  /* Do the raytracing */
  78. double brite[3],*line;  struct world *w;
  79. {
  80.     double t,tmin,pos[3];  int k;
  81.     struct patch ptch;  struct sphere *spnear;
  82.     struct lamp *lmpnear;
  83.  
  84.     tmin=BIG;  spnear=0;        /* can we see some spheres */
  85.     for (k=0; k<w->numsp; ++k)
  86.         if (intsplin(&t,line,w->sp+k)) {
  87.             if (t<tmin) {tmin=t; spnear=w->sp+k;}
  88.         }
  89.     lmpnear=0;                  /* are we looking at a lamp */
  90.     for (k=0; k<w->numlmp; ++k)
  91.         if (intsplin(&t,line,w->lmp+k)) {
  92.             if (t < tmin) {tmin=t; lmpnear=w->lmp+k;}
  93.         }
  94.     if (lmpnear) {              /* we see a lamp! */
  95.         for (k=0; k<3; ++k)
  96.             brite[k]=lmpnear->color[k]/(lmpnear->radius*
  97.                      lmpnear->radius);
  98.          return 0;
  99.         }
  100.     if (inthor(&t,line))        /* do we see the ground? */
  101.         if (t<tmin) {
  102.             point(pos,t,line);  k=gingham(pos); /* cheap vinyl */
  103.             veccopy(w->horizon[k].pos,pos);
  104.             pixbrite(brite,&(w->horizon[k]),w,0);
  105.             return 0;
  106.         }
  107.     if (spnear) {               /* we see a sphere */
  108.         point(ptch.pos,tmin,line);  setnorm(&ptch,spnear);
  109.         colorcpy(ptch.color,spnear->color);
  110.         switch(spnear->type) {  /* treat the surface type */
  111.             case BRIGHT:        /* is it a highlight? */
  112.                 if (glint(brite,&ptch,w,spnear,line)) return 0;
  113.             case DULL:
  114.                 pixbrite(brite,&ptch,w,spnear); return 0;
  115.             case MIRROR:
  116.                 mirror(brite,&ptch,w,line); return 0;
  117.         }
  118.          return 0;
  119.     }
  120.     skybrite(brite,line,w);     /* nothing else, must be sky */
  121. }
  122.  
  123. skybrite(brite,line,w)          /* calculate sky color */
  124. double brite[3],*line;
  125. struct world *w;
  126. {   /* Blend a sky color from the zenith to the horizon */
  127.     double sin2,cos2;  int k;
  128.     sin2=line[5]*line[5];
  129.     sin2/=(line[1]*line[1]+line[3]*line[3]+sin2);
  130.     cos2=1.0-sin2;
  131.     for (k=0; k<3; ++k)
  132.         brite[k]=cos2*w->skyhor[k]+sin2*w->skyzen[k];
  133. }
  134.  
  135. pixline(line,o,i,j)             /* calculate ray for pixel i,j */
  136. double *line;  struct observer *o;  int i,j;
  137. {
  138.     double x,y,tp[3];  int k;
  139.     y=(0.5*o->ny-j)*o->py;
  140.     x=(i-0.5*o->nx)*o->px;
  141.     for (k=0; k<3; ++k)
  142.         tp[k]=o->viewdir[k]*o->fl+y*o->vhat[k]+
  143.               x*o->uhat[k]+o->obspos[k];
  144.     genline(line,o->obspos,tp); /* generate equation of line */
  145. }
  146.  
  147. vecsub(a,b,c)                   /* a=b-c for vectors */
  148. double *a,*b,*c;
  149. {
  150.     int k;
  151.     for (k=0; k<3; ++k) a[k]=b[k]-c[k];
  152. }
  153.  
  154. intsplin(t,line,sp)     /* intersection of sphere and line */
  155. double *t,*line;  struct sphere *sp;
  156. {/* t returns the parameter for where on the line */
  157.     double a,b,c,d,p,q,tt,sqrt();  int k; /* the sphere is hit */
  158.     a=b=0.0;  c=sp->radius; c=-c*c;
  159.     for (k=0; k<3; ++k) {
  160.         p=(*line++)-sp->pos[k];  q=*line++;
  161.         a=q*q+a;  tt=q*p;  b=tt+tt+b;  c=p*p+c;
  162.     } /* a,b,c are coefficients of quadratic equation for t */
  163.     d=b*b-4.0*a*c;
  164.     if (d <= 0) return 0;       /* line misses sphere */
  165.     d=sqrt(d);  *t=-(b+d)/(a+a);
  166.     if (*t<SMALL) *t=(d-b)/(a+a);
  167.     return *t>SMALL;    /* is sphere is in front of us? */
  168. }
  169.  
  170. qintsplin(line,sp)      /* as above, but we don't need t */
  171. double *line;
  172. struct sphere *sp;
  173. {
  174.     double a,b,c,d,p,q;  int k;
  175.     a=b=0.0;  c=sp->radius; c=-c*c;
  176.     for (k=0; k<3; ++k) {
  177.         p=(*line++)-sp->pos[k];  q=*line++;
  178.         a+=q*q;  b+=2.0*p*q;  c+=p*p;
  179.     }
  180.     d=b*b-4.0*a*c;  return d > 0.0;
  181. }
  182.  
  183. inthor(t,line)  /* intersection of line with ground */
  184. double *t,*line;
  185. {
  186.     if (line[5] == 0.0) return 0;
  187.     *t=-line[4]/line[5];  return *t > SMALL;
  188. }
  189.  
  190. genline(l,a,b)  /* generate the equation of a line through the */
  191. double *l,*a,*b;/* two points a and b */
  192. {
  193.     int k;
  194.     for (k=0; k<3; ++k) {*l++=a[k]; *l++=b[k]-a[k];}
  195. }
  196.  
  197. double dot(a,b) /* dot product of 2 vectors */
  198. double *a,*b;
  199. {
  200.  return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
  201. }
  202.  
  203. point(pos,t,line)  /* calculate position of a point on the */
  204. double *pos,t,*line; /* line with parameter t */
  205. {
  206.     int k;  double a;
  207.     for (k=0; k<3; ++k) {
  208.         a=*line++;  pos[k]=a+(*line++)*t;
  209.     }
  210. }
  211.  
  212. glint(brite,p,w,spc,incident)   /* are we looking at a */
  213. double brite[3];                /* highlight? */
  214. struct patch *p;  struct world *w;  struct sphere *spc;
  215. double *incident;
  216. {
  217.     int k,l,firstlite;  static double minglint=0.95;
  218.     double line[6],t,r,lp[3],*pp,*ll,cosi;
  219.     double incvec[3],refvec[3],ref2;
  220.     firstlite=1;
  221.     for (l=0; l<w->numlmp; ++l) {
  222.         ll=(w->lmp+l)->pos;  pp=p->pos;
  223.         vecsub(lp,ll,pp);  cosi=dot(lp,p->normal);
  224.         if (cosi <= 0.0) continue; /* not with this lamp! */
  225.         genline(line,pp,ll);
  226.         for (k=0; k<w->numsp; ++k) {
  227.             if (w->sp+k == spc) continue;
  228.             if (intsplin(&t,line,w->sp+k)) goto cont;
  229.         }
  230.         if (firstlite) {
  231.             incvec[0]=incident[1];  incvec[1]=incident[3];
  232.             incvec[2]=incident[5];
  233.             reflect(refvec,p->normal,incvec);
  234.             ref2=dot(refvec,refvec);  firstlite=0;
  235.         }
  236.         r=dot(lp,lp);  t=dot(lp,refvec);
  237.         t*=t/(dot(lp,lp)*ref2);
  238.         if (t > minglint) { /* it's a highlight */
  239.             for (k=0; k<3; ++k) brite[k]=1.0;
  240.             return 1;
  241.         }
  242.     cont:
  243.     }
  244.     return 0;
  245. }
  246.  
  247. mirror(brite,p,w,incident) /* bounce ray off mirror */
  248. double brite[3];  struct patch *p;
  249. struct world *w;  double *incident;
  250. {
  251.     int k;  double line[6],incvec[3],refvec[3],t;
  252.     incvec[0]=incident[1];  incvec[1]=incident[3];
  253.     incvec[2]=incident[5];  t=dot(p->normal,incvec);
  254.     if (t >= 0) { /* we're inside a sphere, it's dark */
  255.         for (k=0; k<3; ++k) brite[k]=0.0;
  256.         return 0;
  257.     }
  258.     reflect(refvec,p->normal,incvec);  line[0]=p->pos[0];
  259.     line[2]=p->pos[1];  line[4]=p->pos[2];  line[1]=refvec[0];
  260.     line[3]=refvec[1];  line[5]=refvec[2];
  261.     raytrace(brite,line,w);  /* recursion saves the day */
  262.     for (k=0; k<3; ++k) brite[k]=brite[k]*p->color[k];
  263.     return 1;
  264. }
  265.  
  266. pixbrite(brite,p,w,spc)  /* how bright is the patch? */
  267. double brite[3];  struct patch *p;
  268. struct world *w;  struct sphere *spc;
  269. {
  270.     int k,l;  double line[6],t,r,lp[3],*pp,*ll,cosi,diffuse;
  271.     double sqrt();
  272.     static double zenith[3]={0.0,0.0,1.0},f1=1.5,f2=0.4;
  273.     diffuse=(dot(zenith,p->normal)+f1)*f2;
  274.     for (k=0; k<3; ++k) brite[k]=diffuse*w->illum[k]*p->color[k];
  275.     if (p && w) {
  276.         for (l=0; l<w->numlmp; ++l) {
  277.             ll=(w->lmp+l)->pos;  pp=p->pos;  vecsub(lp,ll,pp);
  278.             cosi=dot(lp,p->normal);  if (cosi <= 0.0) goto cont;
  279.             genline(line,pp,ll);
  280.             for (k=0; k<w->numsp; ++k) {
  281.                 if (w->sp+k == spc) continue;
  282.                 if (intsplin(&t,line,w->sp+k)) goto cont;
  283.             }
  284.             r=sqrt(dot(lp,lp));  cosi=cosi/(r*r*r);
  285.             for (k=0; k<3; ++k)
  286.                 brite[k]=brite[k]+cosi*p->color[k]
  287.                          *w->lmp[l].color[k];
  288.             cont:
  289.         }
  290.     }
  291. }
  292.  
  293. setnorm(p,s)    /* normal (radial) direction of sphere */
  294. struct patch *p;  struct sphere *s;
  295. {
  296.     double *t,a;  int k;
  297.     vecsub(t=p->normal,p->pos,s->pos);  a=1.0/s->radius;
  298.     for (k=0; k<3; ++k) {*t=(*t)*a; ++t;}
  299. }
  300.  
  301. colorcpy(a,b)  /* a=b for colors */
  302. double *a,*b;
  303. {
  304.  int k;
  305.  for (k=0; k<3; ++k) a[k]=b[k];
  306. }
  307.  
  308. veccopy(a,b)  /* a=b for vectors */
  309. double *a,*b;
  310. {int k;
  311.  for (k=0; k<3; ++k) a[k]=b[k];
  312. }
  313.  
  314. gingham(pos) /* are we on 'black' or 'white' tile? */
  315. double *pos;
  316. {       /* tiles are 3 units wide */
  317.     double x,y;  int kx,ky;
  318.     kx=ky=0;  x=pos[0]; y=pos[1];
  319.     if (x < 0) {x=-x; ++kx;}
  320.     if (y < 0.0) {y=-y; ++ky;}
  321.     return ((((int)x)+kx)/3+(((int)y)+ky)/3)%2;
  322. }
  323.  
  324. reflect(y,n,x)  /* law of reflection, n is unit normal, */
  325. double *y,*n,*x; /* x is incoming ray, y is outgoing ray */
  326. {
  327.     double u[3],v[3],vv,xn,xv;  int k;
  328.     vecprod(u,x,n);      /* normal to the plane of n and y */
  329.     if (veczero(u)) {    /* bounce right back */
  330.         y[0]=-x[0];  y[1]=-x[1];  y[2]=-x[2];  return 0;
  331.     }
  332.     vecprod(v,u,n);          /* u,v and n are orthogonal */
  333.     vv=dot(v,v);  xv=dot(x,v)/vv;  xn=dot(x,n);
  334.     for (k=0; k<3; ++k) y[k]=xv*v[k]/(xn*n[k]);
  335. }
  336.  
  337. vecprod(a,b,c)          /* vector product a=b^c */
  338. double *a,*b,*c;
  339. {
  340.     a[0]=b[1]*c[2]-b[2]*c[1];
  341.     a[1]=b[2]*c[0]-b[0]*c[2];
  342.     a[2]=b[0]*c[1]-b[1]*c[0];
  343. }
  344.  
  345. veczero(v)              /* is vector null? */
  346. double *v;
  347. {
  348.     if (v[0] != 0.0) return 0;  if (v[1] != 0.0) return 0;
  349.     if (v[2] != 0.0) return 0; return 1; }
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!