home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / graphic / dwb / ray.c < prev    next >
Text File  |  1989-05-13  |  21KB  |  594 lines
  1. /************************************************************************
  2.  *                                                                      *
  3.  *                  Copyright (c) 1987, David B. Wecker                 *
  4.  *                          All Rights Reserved                         *
  5.  *                                                                      *
  6.  * This file is part of DBW_Render                                      *
  7.  *                                                                      *
  8.  * DBW_Render is distributed in the hope that it will be useful, but    *
  9.  * WITHOUT ANY WARRANTY. No author or distributor accepts               *
  10.  * responsibility to anyone for the consequences of using it or for     *
  11.  * whether it serves any particular purpose or works at all, unless     *
  12.  * he says so in writing. Refer to the DBW_Render General Public        *
  13.  * License for full details.                                            *
  14.  *                                                                      *
  15.  * Everyone is granted permission to copy, modify and redistribute      *
  16.  * DBW_Render, but only under the conditions described in the           *
  17.  * DBW_Render General Public License. A copy of this license is         *
  18.  * supposed to have been given to you along with DBW_Render so you      *
  19.  * can know your rights and responsibilities. It should be in a file    *
  20.  * named COPYING. Among other things, the copyright notice and this     *
  21.  * notice must be preserved on all copies.                              *
  22.  ************************************************************************
  23.  *                                                                      *
  24.  * Authors:                                                             *
  25.  *      DBW - David B. Wecker                                           *
  26.  *      JHL - John H. Lowery, IBM-PC/Microsoft C port                   *
  27.  *                                                                      *
  28.  * Versions:                                                            *
  29.  *      V1.0 870125 DBW - First released version                        *
  30.  *      V1.1 (IBM) 890101 JHL - version 1.0 for the IBM-PC/VGA/MCGA     *
  31.  *                                                                      *
  32.  ************************************************************************/
  33.  
  34. #define MODULE_RAY
  35. #include "ray.h"
  36.  
  37. vector    base,
  38.           vx,
  39.           vy,
  40.           d,
  41.           p,
  42.           sumval,
  43.           sampleval,
  44.           varsum,
  45.           samplep,
  46.           rnd_vr,
  47.           rnd_vu,
  48.           subvr,
  49.           subvu,
  50.           old_ave,
  51.           ave_val,
  52.           eye2,
  53.           focalpoint,
  54.           lens1,
  55.           lens2,
  56.           temp1,
  57.           temp2;
  58.  
  59. float     rnd_wid,
  60.           rnd_hi,
  61.           lensx,
  62.           lensy,
  63.           focallength,
  64.           totalcount,
  65.           count,
  66.           maxcount,
  67.           d_alias,
  68.           d_bestf1,
  69.           d_bestf2;
  70.  
  71. unsigned  modulo;
  72. int       row,
  73.           col,
  74.           i,j,k,
  75.           compthresh,
  76.           compmin,
  77.           bunch,
  78.           fac1,
  79.           fac2,
  80.           bestf1,
  81.           bestf2,
  82.           spatialx,
  83.           spatialy,
  84.           pixelcompute,
  85.           pixelguess,
  86.           didguess;
  87. long      curtime,
  88.           prvtime,
  89.           nxttime,
  90.           xtrtime,
  91.           buftime,
  92.           totaltime,
  93.           totalguess,
  94.           totalcompute;
  95.  
  96. unsigned char cache[9][MAXX][3];
  97. int           computes[MAXX];
  98.  
  99.  
  100. guess2(col,sum,val,dis,xdir,ydir)
  101. int     col,*sum,dis,xdir,ydir;
  102. int     val[3];
  103. {
  104.      int i,x1,y1,x2,y2,retval,shft,isum;
  105.  
  106.      shft    = 6 - dis;
  107.      isum    = 1 << shft;
  108.      y1      = 4 - (dis * ydir);
  109.      x1      = col - (dis * xdir);
  110.      if (y1 < 0)
  111.           y1 = 0;
  112.      if (y1 > 8)
  113.           y1 = 8;
  114.      if (x1 < 0)
  115.           x1 = 0;
  116.      if (x1 >= MAXCOL)
  117.           x1 = MAXCOL-1;
  118.  
  119.      retval = 0;
  120.      y2 = y1;
  121.      while (1) 
  122.      {
  123.           if (cache[y2][x1][0] != 0xFF) 
  124.           {
  125.                for (i = 0; i < 3; i++)
  126.                     val[i] += ((int)cache[y2][x1][i]) << shft;
  127.                *sum    += isum;
  128.                retval   = 1;
  129.           }
  130.           if (y2 == 4)
  131.                break;
  132.           y2 += ydir;
  133.      }
  134.      x2 = x1;
  135.      while (1) 
  136.      {
  137.           if (cache[y1][x2][0] != 0xFF) 
  138.           {
  139.                for (i = 0; i < 3; i++)
  140.                     val[i] += ((int)cache[y1][x2][i]) << shft;
  141.                *sum    += isum;
  142.                retval   = 1;
  143.           }
  144.           if (x2 == col)
  145.                break;
  146.           x2 += xdir;
  147.      }
  148.      return(retval);
  149. }
  150.  
  151. void guess(col)
  152. int col;
  153. {
  154.      int     i,sum,dir,dis;
  155.      int     val[3];
  156.  
  157.      sum = 0;
  158.      dir = 0;
  159.      for (i = 0; i < 3; i++)
  160.           val[i] = 0;
  161.      for (dis = 1; dir != 15 && dis < 7; dis++) 
  162.      {
  163.           if ((dir & 1) != 1)
  164.                if (guess2(col,&sum,val,dis,-1,1))
  165.                     dir |= 1;
  166.           if ((dir & 2) != 2)
  167.                if (guess2(col,&sum,val,dis,1,-1))
  168.                     dir |= 2;
  169.           if ((dir & 4) != 4)
  170.                if (guess2(col,&sum,val,dis,1,1))
  171.                     dir |= 4;
  172.           if ((dir & 8) != 8)
  173.                if (guess2(col,&sum,val,dis,-1,-1))
  174.                     dir |= 8;
  175.           if (dis > 2 && ((dir & 3) == 3 || (dir & 12) == 12))
  176.                break;
  177.      }
  178.      if (sum == 0)
  179.           sum = 1;
  180.      for (i = 0; i < 3; i++) 
  181.      {
  182.           cache[4][col][i] = (unsigned char)(val[i] / sum);
  183.      }
  184. }
  185.  
  186. void do_raytrace()
  187. {
  188.  
  189.      /* set up pixel coordinate */
  190.      vecscale((float)col,vr,vx);
  191.      vecsub(base,vy,p);
  192.      vecsum(p,vx,p);
  193.  
  194.      /* Compute this next pixel.  Distribute 'n' rays to do this */
  195.      veczero(sumval);  /* start with a black pixel */
  196.      veczero(varsum);
  197.      veccopy(backgroundval,ave_val);  /* have to start someplace */
  198.      count = 0.0;
  199.      bunch = antialias;
  200.      spatialx = 0;
  201.      spatialy = 0;
  202.      curr_runs = 0;
  203.  
  204.      while (bunch > 0) 
  205.      {
  206.  
  207.           /* Cast another sample ray for this pixel into the next subregion */
  208.           /* Pick a new random perturbing vector for this sample for spatial
  209.           antialiasing */
  210.           if (count < d_alias) 
  211.           {
  212.                rnd_wid = 0.5;  /* force first bunch through region centers */
  213.                rnd_hi  = 0.5;
  214.           }
  215.           else
  216.           {
  217.                rnd_wid = rnd();
  218.                rnd_hi  = rnd();
  219.           }
  220.  
  221.           vecscale(rnd_wid / d_bestf2,vr,rnd_vr);
  222.           vecscale(rnd_hi  / d_bestf1,vu,rnd_vu);
  223.           vecscale((float) spatialx / d_bestf2,vr,subvr);
  224.           vecscale((float) spatialy / d_bestf1,vu,subvu);
  225.           vecsum(rnd_vr,subvr,samplep);
  226.           vecsum(samplep,p,samplep);
  227.           vecsub(samplep,rnd_vu,samplep);
  228.           vecsub(samplep,subvu,samplep);
  229.           veccopy(eye,eye2);
  230.           direction(eye2,samplep,d);  /* direction from lens center to pixel */
  231.  
  232.           if (aperture > 0.0) 
  233.           {
  234.                vecscale(focus,d,focalpoint);  /* relative point on focus plane */
  235.                lensx = rnd();  /* pick random point on lens */
  236.                lensy = rnd();
  237.                if (rnd() < 0.5)
  238.                     lensx = -lensx;
  239.                if (rnd() < 0.5)
  240.                     lensy = -lensy;
  241.                vecscale(lensx,lens1,temp1);  /* make into vectors in lens plane */
  242.                vecscale(lensy,lens2,temp2);
  243.                vecsum(temp1,temp2,eye2);  /* point in the lens */
  244.                direction(eye2,focalpoint,d);  /* from spot on lens to focal point */
  245.                vecsum(eye,eye2,eye2);  /* calc new 3D eye point in lens */
  246.           }
  247.  
  248.           /* allow aborting if curr_runs > max_runs */
  249.           if (setjmp(env)) 
  250.           {
  251.                computes[col] = max_runs;    /* too many computes */
  252.                cache[8][col][1] = 0;   /* make sure we don't try again */
  253.                return;
  254.           }
  255.  
  256.           /* perform the ray trace */
  257.           dodirection(sampleval,eye2,d,1.0,ambientlight);
  258.  
  259.           count      += 1.0;  /* increment the ray count */
  260.           totalcount += 1.0;
  261.           bunch--;  /* Count down in this bunch */
  262.           spatialx++;  /* move across to next subregion */
  263.           if (spatialx == bestf2) 
  264.           {
  265.                spatialx = 0;  /* move back & down to start of next subregion scan */
  266.                spatialy++;
  267.                if (spatialy == bestf1)
  268.                     spatialy = 0;  /* all done,start at beginning again */
  269.           }
  270.  
  271.           vecsum(sampleval,sumval,sumval);    /* accumulate the samples */
  272.           vecscale(1.0/count,sumval,ave_val); /* compute new average */
  273.           vecsub(sampleval,ave_val,sampleval);/* compute nth variance */
  274.           vecmul(sampleval,sampleval,sampleval); /* square it */
  275.           vecsum(sampleval,varsum,varsum);    /* accumulate variances */
  276.  
  277.           if (bunch == 0) 
  278.           {
  279.                /* When bunch=0,perform the sufficiency test */
  280.                if ((varsum[0] / count) > variance ||
  281.                  (varsum[1] / count) > variance ||
  282.                  (varsum[2] / count) > variance)
  283.                     if (count < 128)  /* maximum of 127+bunch samples */
  284.                          bunch = antialias;  /* do another bunch */
  285.           }
  286.      }  /* while distributing */
  287.  
  288.      if (count > maxcount)
  289.           maxcount = count;  /* new champion pixel */
  290.  
  291.      if (histogram) 
  292.      {
  293.           cv(count/d_alias,count/d_alias,count/d_alias,ave_val);
  294.      }
  295.      else
  296.      {
  297.           vecscale(d_maxgray,ave_val,ave_val);
  298.      }
  299.  
  300.      /* store the pixel away */
  301.      for (i = 0; i < 3; i++) 
  302.      {
  303.           if (ave_val[i] >= d_maxgray)
  304.                ave_val[i] = d_maxgray-1.0;
  305.           cache[8][col][i] = (char) ave_val[i];
  306.      }
  307.  
  308.      /* save full statistics */
  309.      pixelcompute++;
  310.      total_runs += curr_runs;
  311.      computes[col] = curr_runs;
  312. }
  313.  
  314. main(argc,argv)
  315. int     argc;
  316. char    *argv;
  317. {
  318.      stacktop = curstack();
  319.      stackbot = stacktop;
  320.  
  321.      printf(VERSION);
  322.  
  323.      getinput(argc,argv);
  324.      getoutput(argc,argv);
  325.  
  326.      fprintf(fpout, VERSION);
  327.  
  328.      brickmortarwidth  = brickmortar / brickwidth;
  329.      brickmortarheight = brickmortar / brickheight;
  330.      brickmortardepth  = brickmortar / brickdepth;
  331.      d_maxgray = (float) MAXGRAY;
  332.      sqrt3 = (float)sqrt(3.0);
  333.      sqrt3times2 = 2.0 * sqrt3;
  334.      d_alias = (float) antialias;
  335.  
  336.      /* for simulating depth-of-field,we distribute rays scattered from an
  337.      imaginary lens disk,out through a focal point,and into the scene.
  338.      In order to scatter the rays over the lens,we need two vectors
  339.      that are orthogonal to themselves and the central direction of view.
  340.      In essence,these two vectors define the film plane.  Just such two
  341.      vectors were calculated in FIL.C for the EYE point.  */
  342.  
  343.      veccopy(vu,lens1);
  344.      veccopy(vr,lens2);
  345.      vecscale(aperture,lens1,lens1);  /* Scale them accordingly */
  346.      vecscale(aperture,lens2,lens2);
  347.  
  348.  
  349.      vecscale((float) (MAXCOL / DISPWID) / (MAXROW / DISPHI),vu,vu);
  350.      vecscale((float) (MAXROW / 2),vu,vy);
  351.      vecscale((float) (-MAXCOL /2),vr,vx);
  352.      vecsum(vrp,vy,base);
  353.      vecsum(base,vx,base);
  354.  
  355.      bestf1 = 1;  /* Find the best integral factors of antialias */
  356.      bestf2 = antialias;
  357.      for (fac1 = 2; fac1 <= antialias / 2; fac1++) 
  358.      {
  359.           fac2 = antialias / fac1;
  360.           if (fac1 * fac2 == antialias)
  361.                if (abs(fac1 - fac2) < abs(bestf1 - bestf2)) 
  362.                {
  363.                     bestf1 = fac1;
  364.                     bestf2 = fac2;
  365.                }
  366.      }
  367.      d_bestf1        = (float) bestf1;
  368.      d_bestf2        = (float) bestf2;
  369.      totalcount      = 0.0;
  370.      maxcount        = 0.0;
  371.      sline           = startrow;
  372.      buftime         = (long)(maxhours * 3600.0 / 408.0);
  373.      totaltime       = 0L;
  374.      curtime         = time(0L);
  375.      xtrtime         = 0L;
  376.      totalguess      = 0L;
  377.      totalcompute    = 0L;
  378.      pixelcompute    = MAXCOL;
  379.  
  380.      printf("\nApproximately %ld seconds per scan line\n",buftime);
  381.      fprintf(fpout,"\nApproximately %ld seconds per scan line\n",buftime);
  382.  
  383.      /* start out with no guesses as to computing difficulty */
  384.      for (i = 0; i < MAXCOL; i++)
  385.           computes[i] = 0;
  386.  
  387.      /* define the cache scan lines as "unknown" */
  388.      for (i = 0; i < 9; i++)
  389.           for (j = 0; j < MAXCOL; j++)
  390.                cache[i][j][0] = cache[i][j][1] = 0xFF;
  391.  
  392.      /* do each scan line */
  393.      for (row = startrow-4; row < endrow+4; row++) 
  394.      {
  395.  
  396.           /* set up stats for this line */
  397.           compthresh      = (pixelcompute * 2) / 3;
  398.           compmin         = (pixelcompute * 9) / 10;
  399.           pixelguess      = 0;
  400.           pixelcompute    = 0;
  401.  
  402.           /* Decide how much time we have */
  403.           prvtime    = curtime;
  404.           nxttime    = prvtime + buftime + (xtrtime >> 3);
  405.           xtrtime   -= xtrtime >> 3;
  406.  
  407.           if ((row+4-startrow) % 60 == 0)
  408.                printf("\nRow %3d: ",row);
  409.           printf(".");
  410.           fflush(stdout);
  411.  
  412.           /* premultiply out the desired row (for do_raytrace) */
  413.           vecscale((float)row,vu,vy);
  414.  
  415.           /* get random pixels until out of time */
  416.           curtime  = time(0L);
  417.           max_runs = 10;
  418.           while (curtime <= nxttime || pixelcompute <= compmin) 
  419.           {
  420.  
  421.                /* find remaining pixel with smallest estimated computes */
  422.                i = j = (int)(rnd() * (float)MAXCOL);
  423.                k = 32767;
  424.                col = -1;
  425.                do 
  426.                {
  427.                     /* see if we have a possible best new pixel on the line*/
  428.                     if (cache[8][i][0] == 0xFF) 
  429.                     {
  430.                          if (pixelcompute > compthresh) 
  431.                          {
  432.                               col      = i;
  433.                               max_runs = 32767;
  434.                               break;
  435.                          }
  436.                          if (cache[8][i][1] == 0xFF) 
  437.                          {
  438.                               if (col == -1 || k > computes[i]) 
  439.                               {
  440.                                    col = i;
  441.                                    k   = computes[i];
  442.                               }
  443.                          }
  444.                     }
  445.                     if (++i >= MAXCOL)
  446.                          i = 0;
  447.                }
  448.                while (i != j)
  449.                     ;
  450.  
  451.                /* Couldn't find a pixel.... */
  452.                if (col == -1) 
  453.                {
  454.  
  455.                     /* See if there are any that need more computing */
  456.                     j = 1;
  457.                     for (i = 0; i < MAXCOL; i++)
  458.                          if (cache[8][i][0] == 0xFF) 
  459.                          {
  460.                               cache[8][i][1] = 0xFF;
  461.                               j = 0;
  462.                          }
  463.                     if (j)
  464.                          goto DONE_COMPUTE;
  465.  
  466.                     /* see if we can TRY more computing */
  467.                     if (max_runs == 10)
  468.                          max_runs = 32767;
  469.                     else
  470.                          goto DONE_COMPUTE;
  471.                }
  472.                else 
  473.                {
  474.                     if (kbhit())
  475.                          if (getch() == 0x1B)
  476.                          {
  477.                               printf("\nESCAPE Received");
  478.                               fprintf(fpout,"\nESCAPE Received");
  479.                               goto OPERATOR_ABORT;
  480.                          }
  481.                     do_raytrace();
  482.                }
  483.                curtime = time(0L);
  484.           }
  485.  
  486.           /* all done computing */
  487. DONE_COMPUTE:
  488.           totalcompute += (long)pixelcompute;
  489.  
  490.           /* ignore until we get 4 lines into the cache */
  491.           if (row >= startrow) 
  492.           {
  493.  
  494.                modulo = 0;
  495.                for (col = 0; col < MAXCOL; col++) 
  496.                {
  497.                     if (cache[4][col][0] == 0xFF) 
  498.                     {
  499.                          ++pixelguess;
  500.                          guess(col);
  501.                          didguess = 1;
  502.                     }
  503.                     else
  504.                          didguess = 0;
  505.                     i = col / PPW;
  506.                     if (modulo) 
  507.                     {
  508.                          blueline[i]  |= ((int)cache[4][col][0]) << modulo;
  509.                          greenline[i] |= ((int)cache[4][col][1]) << modulo;
  510.                          redline[i]   |= ((int)cache[4][col][2]) << modulo;
  511.                     }
  512.                     else
  513.                     {
  514.                          blueline[i]  = (int) cache[4][col][0];
  515.                          greenline[i] = (int) cache[4][col][1];
  516.                          redline[i]   = (int) cache[4][col][2];
  517.                     }
  518.                     modulo = (((modulo >> 2) + 1) % PPW) << 2;
  519.                     if (didguess)
  520.                          cache[4][col][0] = 0xFF;
  521.                }
  522.  
  523.                /* dump the scan line to the file */
  524.                write_scanline();
  525.           }
  526.  
  527.           curtime     = time(0L);
  528.           xtrtime    += nxttime - curtime;
  529.           nxttime     = curtime - prvtime;
  530.           totaltime  += nxttime;
  531.           totalguess += (long)pixelguess;
  532.  
  533.           /* move the cache entries down */
  534.           for (i = 1; i < 9; i++)
  535.                for (j = 0; j < MAXCOL; j++)
  536.                     if (cache[i][j][0] == 0xFF)
  537.                          cache[i-1][j][0] = cache[i-1][j][1] = 0xFF;
  538.                     else
  539.                          for (k = 0; k < 3; k++)
  540.                               cache[i-1][j][k] = cache[i][j][k];
  541.  
  542.           /* define the last cache scan line as "unknown" */
  543.           for (j = 0; j < MAXCOL; j++)
  544.                cache[8][j][0] = cache[8][j][1] = 0xFF;
  545.      } 
  546.  
  547. OPERATOR_ABORT:
  548.  
  549.      printf("\n");
  550.      fprintf(fpout,"\n");
  551.      if (fclose(fp))
  552.           ERROR( "Error closing output file." );
  553.  
  554.      printf("First Line                 = %d\n", startrow);
  555.      printf("Last Line                  = %d\n", sline);
  556.      printf("Average distribution count = %4.2lf, max = %1.0lf\n",
  557.        ((float)totalcount)/((float)(endrow-startrow)*(float)MAXCOL),
  558.        maxcount);
  559.      printf("Total pixels computed      = %ld\n",totalcompute);
  560.      printf("Total pixels guessed       = %ld\n",totalguess);
  561.      printf("Total time used            = %ld seconds\n",totaltime);
  562.      printf("Ray intersect runs         = %ld\n",total_runs);
  563.      printf("Max intersects for 1 ray   = %ld\n",max_intersects);
  564.      printf("Avg intersects per pixel   = %4.2lf\n",
  565.        ((float)total_runs)/((float)(endrow-startrow)*(float)MAXCOL));
  566.      printf("Total sorts                = %ld\n",sorts);
  567.      if (sorts > 0)
  568.           printf("Average sort size          = %4.2lf\n",
  569.             ((float)sort_size) / (float)sorts);
  570.      stacktop -= stackbot;
  571.      stacktop += 2000L;
  572.      printf("Maximum stack size         = %ld\n",stacktop);
  573.  
  574.      /* statistics file */
  575.      fprintf(fpout,"First Line                 = %d\n", startrow);
  576.      fprintf(fpout,"Last Line                  = %d\n", sline);
  577.      fprintf(fpout,"Average distribution count = %4.2lf, max = %1.0lf\n",
  578.        ((float)totalcount)/((float)(endrow-startrow)*(float)MAXCOL),
  579.        maxcount);
  580.      fprintf(fpout,"Total pixels computed      = %ld\n",totalcompute);
  581.      fprintf(fpout,"Total pixels guessed       = %ld\n",totalguess);
  582.      fprintf(fpout,"Total time used            = %ld seconds\n",totaltime);
  583.      fprintf(fpout,"Ray intersect runs         = %ld\n",total_runs);
  584.      fprintf(fpout,"Max intersects for 1 ray   = %ld\n",max_intersects);
  585.      fprintf(fpout,"Avg intersects per pixel   = %4.2lf\n",
  586.        ((float)total_runs)/((float)(endrow-startrow)*(float)MAXCOL));
  587.      fprintf(fpout,"Total sorts                = %ld\n",sorts);
  588.      if (sorts > 0)
  589.           fprintf(fpout,"Average sort size          = %4.2lf\n",
  590.             ((float)sort_size) / (float)sorts);
  591.      fprintf(fpout,"Maximum stack size         = %ld\n",stacktop);
  592.  
  593. }
  594.