home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 3 / Amiga Tools 3.iso / grafik / raytracing / rayshade-4.0.6.3 / libray / libobj / fl.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-09  |  7.7 KB  |  310 lines
  1. /*
  2.  *    fl.c
  3.  *
  4.  *    fracland: a fractal landscape primitive for rayshade.
  5.  *
  6.  *     Copyright (C) 1992, Lawrence K. Coffin.
  7.  *    Parts Copyright (C) 1989, 1991, Craig E. Kolb.
  8.  *    All rights reserved.
  9.  * 
  10.  *    Parts of this program are from "Numerical Recipes in C" by William H.
  11.  *    Press, Brian P. Flannery, Saul A. Teukolsky, and William T. Vetterling.
  12.  *
  13.  *    This software may be freely copied, modified, and redistributed
  14.  *    provided that this copyright notice is preserved on all copies.
  15.  *
  16.  *    You may not distribute this software, in whole or in part, as part of
  17.  *    any commercial product without the express consent of the authors.
  18.  * 
  19.  *    There is no warranty or other guarantee of fitness of this software
  20.  *    for any purpose.  It is provided solely "as is".
  21.  */
  22.  
  23. #include "geom.h"
  24. #include "hf.h"
  25. #include "fl.h"
  26. #include <math.h>
  27. #if defined(AMIGA) && defined(__GNUC__)
  28. #define pow mypow
  29. extern double mypow(const double x, const double y);
  30. #endif
  31. #include <stdlib.h>
  32.  
  33. static Methods *iHfMethods = NULL;
  34. static char hfName[] = "heighfield";
  35.  
  36. static void integrate_grid(), QueueTri();
  37. static int DDA2D(), CheckCell();
  38. static Float intHftri();
  39. static float minalt(), maxalt();
  40.  
  41. typedef struct {
  42.     int stepX, stepY;
  43.     Float tDX, tDY;
  44.     float minz, maxz;
  45.     int outX, outY;
  46.     Vector cp, pDX, pDY;
  47. } Trav2D;
  48.  
  49. hfTri *CreateHfTriangle(), *GetQueuedTri();
  50.  
  51. extern unsigned long HFTests, HFHits;
  52.  
  53. float rv();
  54.  
  55. Hf *
  56. FlCreate(seed,subdiv)
  57. long int seed;
  58. int subdiv;
  59. {
  60.     Hf *hf;
  61.     FILE *fp;
  62.     float val, *maxptr, *minptr;
  63.     int i, j, n;
  64.     int in,i2,x,y;
  65.     float rc, temp;
  66.  
  67.  
  68.     seednrand(seed);
  69.  
  70.     hf = (Hf *)Malloc(sizeof(Hf));
  71.     /*
  72.      * Make the following an option someday.
  73.      */
  74.     hf->BestSize = BESTSIZE;
  75.     /*
  76.      * Store the inverse for faster computation.
  77.      */
  78.     hf->iBestSize = 1. / (float)hf->BestSize;
  79.     /*
  80.      * Get HF size.
  81.      */
  82.     n = pow(2.0,(double)subdiv);
  83.     hf->size = n+1;
  84.         
  85.     hf->data = (float **)share_malloc(hf->size * sizeof(float *));
  86.  
  87.     for(i = 0; i < hf->size; i++){
  88.         hf->data[i] = (float *)share_malloc(hf->size * sizeof(float));
  89.     }
  90.  
  91.  
  92.  
  93.         hf->data[0][0] = 0;
  94.         hf->data[0][n] = 0;
  95.         hf->data[n][0] = 0;
  96.         hf->data[n][n] = 0;
  97.         for (i=0; i<subdiv; i++){
  98.                 in = (double)n/pow(2.0,(double)i);
  99.                 i2 = in/2;
  100.                 rc = 1.0/pow(2.0,(double)i);
  101.                 for (x=0; x<n; x+=in){
  102.                         for (y=0; y<n; y+=in){
  103.                 hf->data[x+i2][y   ] = (hf->data[x   ][y   ] + hf->data[x+in][y   ])/2 + rv()*rc;
  104.                 hf->data[x   ][y+i2] = (hf->data[x   ][y   ] + hf->data[x   ][y+in])/2 + rv()*rc;
  105.                 hf->data[x+in][y+i2] = (hf->data[x+in][y   ] + hf->data[x+in][y+in])/2 + rv()*rc;
  106.                 hf->data[x+i2][y+in] = (hf->data[x   ][y+in] + hf->data[x+in][y+in])/2 + rv()*rc;
  107.                 hf->data[x+i2][y+i2] = (hf->data[x][y] + hf->data[x+in][y] + hf->data[x][y+in] + hf->data[x+in][y+in])/4 + rv()*rc;
  108.                         }
  109.                 }
  110.         }
  111. /*
  112.  *  This flips the points around so the surfaces created here are the same
  113.  *    as the surfaces created with the program "fracland".
  114.  */
  115.     for (x = 0; x < hf->size; x++){
  116.         for (y = x; y < hf->size; y++){
  117.             temp = hf->data[x][y];
  118.             hf->data[x][y] = hf->data[y][x];
  119.             hf->data[y][x] = temp;
  120.         }
  121.     }
  122.  
  123.     for (i = 0; i < hf->size; i++) {
  124.         for (j = 0; j < hf->size; j++) {
  125.             val = hf->data[i][j];
  126.             if (val <= HF_UNSET) {
  127.                 hf->data[i][j] = HF_UNSET;
  128.                 /*
  129.                  * Don't include the point in min/max
  130.                  * calculations.
  131.                  */
  132.                 continue;
  133.             }
  134.             if (val > hf->maxz)
  135.                 hf->maxz = val;
  136.             if (val < hf->minz)
  137.                 hf->minz = val;
  138.         }
  139.     }
  140.  
  141.     /*
  142.      * Allocate levels of grid.  hf->levels = log base BestSize of hf->size
  143.      */
  144.     for (i = hf->size, hf->levels = 0; i > hf->BestSize; i /= hf->BestSize,
  145.                 hf->levels++)
  146.             ;
  147.     hf->levels++;
  148.     hf->qsize = CACHESIZE;
  149.     hf->q = (hfTri **)Calloc((unsigned)hf->qsize, sizeof(hfTri *));
  150.     hf->qtail = 0;
  151.  
  152.     hf->lsize = (int *)share_malloc(hf->levels * sizeof(int));
  153.     hf->spacing = (float *)share_malloc(hf->levels * sizeof(float));
  154.     hf->boundsmax = (float ***)share_malloc(hf->levels * sizeof(float **));
  155.     hf->boundsmin = (float ***)share_malloc(hf->levels * sizeof(float **));
  156.  
  157.     hf->spacing[0] = hf->size -1;
  158.     hf->lsize[0] = (int)hf->spacing[0];
  159.     hf->boundsmax[0] = (float **)share_malloc(hf->lsize[0]*sizeof(float *));
  160.     hf->boundsmin[0] = (float **)share_malloc(hf->lsize[0]*sizeof(float *));
  161.     /*
  162.      * Compute initial bounding boxes
  163.      */
  164.     for (i = 0; i < hf->lsize[0]; i++) {
  165.         hf->boundsmax[0][i]=(float *)share_malloc(hf->lsize[0]*sizeof(float));
  166.         hf->boundsmin[0][i]=(float *)share_malloc(hf->lsize[0]*sizeof(float));
  167.         maxptr = hf->boundsmax[0][i];
  168.         minptr = hf->boundsmin[0][i];
  169.         for (j = 0; j < hf->lsize[0]; j++) {
  170.             *maxptr++ = maxalt(i, j, hf->data) + EPSILON;
  171.             *minptr++ = minalt(i, j, hf->data) - EPSILON;
  172.         }
  173.     }
  174.  
  175.     for (i = 1; i < hf->levels; i++) {
  176.         hf->spacing[i] = hf->spacing[i-1] * hf->iBestSize;
  177.         hf->lsize[i] = (int)hf->spacing[i];
  178.         if ((Float)hf->lsize[i] != hf->spacing[i])
  179.             hf->lsize[i]++;
  180.         hf->boundsmax[i]=(float **)share_malloc(hf->lsize[i]*sizeof(float *));
  181.         hf->boundsmin[i]=(float **)share_malloc(hf->lsize[i]*sizeof(float *));
  182.         for (j = 0; j < hf->lsize[i]; j++) {
  183.             hf->boundsmax[i][j] = (float *)share_malloc(hf->lsize[i] *
  184.                             sizeof(float));
  185.             hf->boundsmin[i][j] = (float *)share_malloc(hf->lsize[i] *
  186.                             sizeof(float));
  187.         }
  188.         integrate_grid(hf, i);
  189.     }
  190.  
  191.     hf->boundbox[LOW][X] = hf->boundbox[LOW][Y] = 0.0;
  192.     hf->boundbox[HIGH][X] = hf->boundbox[HIGH][Y] = 1.0;
  193.     hf->boundbox[LOW][Z] = hf->minz;
  194.     hf->boundbox[HIGH][Z] = hf->maxz;
  195.  
  196.     return hf;
  197. }
  198. /*
  199.  * Generate a random value from -0.5 to 0.5
  200.  */
  201.  
  202. float rv(){
  203.         float y;
  204.         y = (float) nrand() - 0.5;
  205.         return (y);
  206. }
  207.  
  208. Methods *
  209. FlMethods()
  210. {
  211.     if (iHfMethods == (Methods *)NULL) {
  212.         iHfMethods = MethodsCreate();
  213.         iHfMethods->create = (GeomCreateFunc *)FlCreate;
  214.         iHfMethods->methods = FlMethods;
  215.         iHfMethods->name = HfName;
  216.         iHfMethods->intersect = HfIntersect;
  217.         iHfMethods->normal = HfNormal;
  218.         iHfMethods->uv = HfUV;
  219.         iHfMethods->bounds = HfBounds;
  220.         iHfMethods->stats = HfStats;
  221.         iHfMethods->checkbounds = TRUE;
  222.         iHfMethods->closed = FALSE;
  223.     }
  224.  
  225.     return iHfMethods;
  226. }
  227.  
  228.  
  229. /*
  230.  * Build min/max altitude value arrays for the given grid level.
  231.  */
  232. static void
  233. integrate_grid(hf, level)
  234. Hf *hf;
  235. int level;
  236. {
  237.     int i, j, k, l, ii, ji;
  238.     float max_alt, min_alt;
  239.     float **maxinto, **mininto, **frommax, **frommin, *minptr, *maxptr;
  240.     int insize, fromsize, fact;
  241.  
  242.     maxinto = hf->boundsmax[level];
  243.     mininto = hf->boundsmin[level];
  244.     insize = hf->lsize[level];
  245.     frommax = hf->boundsmax[level-1];
  246.     frommin = hf->boundsmin[level-1];
  247.     fact = hf->BestSize;
  248.     fromsize = hf->lsize[level-1];
  249.  
  250.     ii = 0;
  251.  
  252.     for (i = 0; i < insize; i++) {
  253.         ji = 0;
  254.         for (j = 0; j < insize; j++) {
  255.             max_alt = HF_UNSET;
  256.             min_alt = -HF_UNSET;
  257.             for (k = 0; k <= fact; k++) {
  258.                 if (ii+k >= fromsize)
  259.                     continue;
  260.                 maxptr = &frommax[ii+k][ji];
  261.                 minptr = &frommin[ii+k][ji];
  262.                 for (l = 0; l <= fact; l++,maxptr++,minptr++) {
  263.                     if (ji+l >= fromsize)
  264.                         continue;
  265.                     if (*maxptr > max_alt)
  266.                         max_alt = *maxptr;
  267.                     if (*minptr < min_alt)
  268.                         min_alt = *minptr;
  269.                 }
  270.             }
  271.             maxinto[i][j] = max_alt + EPSILON;
  272.             mininto[i][j] = min_alt - EPSILON;
  273.             ji += fact;
  274.         }
  275.         ii += fact;
  276.     }
  277. }
  278. /*
  279.  * Return maximum height of cell indexed by y,x.  This could be done
  280.  * as a macro, but many C compliers will choke on it.
  281.  */
  282. static float
  283. minalt(y,x,data)
  284. int x, y;
  285. float **data;
  286. {
  287.     float  min_alt;
  288.  
  289.     min_alt = min(data[y][x], data[y+1][x]);
  290.     min_alt = min(min_alt, data[y][x+1]);
  291.     min_alt = min(min_alt, data[y+1][x+1]);
  292.     return min_alt;
  293. }
  294.  
  295. /*
  296.  * Return maximum cell height, as above.
  297.  */
  298. static float
  299. maxalt(y,x,data)
  300. int x, y;
  301. float **data;
  302. {
  303.     float  max_alt;
  304.  
  305.     max_alt = max(data[y][x], data[y+1][x]);
  306.     max_alt = max(max_alt, data[y][x+1]);
  307.     max_alt = max(max_alt, data[y+1][x+1]);
  308.     return max_alt;
  309. }
  310.