home *** CD-ROM | disk | FTP | other *** search
/ Otherware / Otherware_1_SB_Development.iso / mac / developm / source / macraysh.sit / Code / Source / hf.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-15  |  16.9 KB  |  738 lines
  1. /*
  2.  * hf.c
  3.  *
  4.  * Copyright (C) 1989, 1991, Craig E. Kolb
  5.  * All rights reserved.
  6.  *
  7.  * This software may be freely copied, modified, and redistributed
  8.  * provided that this copyright notice is preserved on all copies.
  9.  *
  10.  * You may not distribute this software, in whole or in part, as part of
  11.  * any commercial product without the express consent of the authors.
  12.  *
  13.  * There is no warranty or other guarantee of fitness of this software
  14.  * for any purpose.  It is provided solely "as is".
  15.  *
  16.  * $Id: hf.c,v 4.0 91/07/17 14:38:15 kolb Exp Locker: kolb $
  17.  *
  18.  * $Log:    hf.c,v $
  19.  * Revision 4.0  91/07/17  14:38:15  kolb
  20.  * Initial version.
  21.  * 
  22.  */
  23. #include "geom.h"
  24. #include "hf.h"
  25.  
  26. Methods *iHfMethods = NULL;
  27. static char hfName[] = "heighfield";
  28.  
  29. static void integrate_grid(), QueueTri();
  30. static int DDA2D(), CheckCell();
  31. static Float intHftri();
  32. static float minalt(), maxalt();
  33.  
  34. typedef struct {
  35.     int stepX, stepY;
  36.     Float tDX, tDY;
  37.     float minz, maxz;
  38.     int outX, outY;
  39.     Vector cp, pDX, pDY;
  40. } Trav2D;
  41.  
  42. hfTri *CreateHfTriangle(), *GetQueuedTri();
  43.  
  44. unsigned long HFTests, HFHits;
  45.  
  46. Hf *
  47. HfCreate(filename)
  48. char *filename;
  49. {
  50.     Hf *hf;
  51.     FILE *fp;
  52.     float val, *maxptr, *minptr;
  53.     int i, j;
  54.  
  55.     fp = fopen(filename, "r");
  56.     if (fp == (FILE *)NULL) {
  57.         RLerror(RL_ABORT, "Cannot open heightfield file \"%s\".",
  58.                 filename,0,0);
  59.         return (Hf *)NULL;
  60.     }
  61.  
  62.     hf = (Hf *)Malloc(sizeof(Hf));
  63.     /*
  64.      * Make the following an option someday.
  65.      */
  66.     hf->BestSize = BESTSIZE;
  67.     /*
  68.      * Store the inverse for faster computation.
  69.      */
  70.     hf->iBestSize = 1. / (float)hf->BestSize;
  71.     /*
  72.      * Get HF size.
  73.      */
  74.     if (fread((char *)&hf->size, sizeof(int), 1, fp) == 0) {
  75.         RLerror(RL_ABORT, "Cannot read height field size.",0,0,0); 
  76.         return (Hf *)NULL;
  77.     }
  78.         
  79.     hf->data = (float **)share_malloc(hf->size * sizeof(float *));
  80.     for (i = 0; i < hf->size; i++) {
  81.         hf->data[i] = (float *)share_malloc(hf->size * sizeof(float));
  82.         /*
  83.          * Read in row of HF data.
  84.          */
  85.         if (fread((char *)hf->data[i],sizeof(float),hf->size,fp)
  86.             != hf->size) {
  87. /*            RLerror(RL_ABORT, "Not enough heightfield data."); */
  88.             return (Hf *)NULL;
  89.         }
  90.         for (j = 0; j < hf->size; j++) {
  91.             val = hf->data[i][j];
  92.             if (val <= HF_UNSET) {
  93.                 hf->data[i][j] = HF_UNSET;
  94.                 /*
  95.                  * Don't include the point in min/max
  96.                  * calculations.
  97.                  */
  98.                 continue;
  99.             }
  100.             if (val > hf->maxz)
  101.                 hf->maxz = val;
  102.             if (val < hf->minz)
  103.                 hf->minz = val;
  104.         }
  105.     }
  106.     (void)fclose(fp);
  107.     /*
  108.      * Allocate levels of grid.  hf->levels = log base BestSize of hf->size
  109.      */
  110.     for (i = hf->size, hf->levels = 0; i > hf->BestSize; i /= hf->BestSize,
  111.                 hf->levels++)
  112.             ;
  113.     hf->levels++;
  114.     hf->qsize = CACHESIZE;
  115.     hf->q = (hfTri **)Calloc((unsigned)hf->qsize, sizeof(hfTri *));
  116.     hf->qtail = 0;
  117.  
  118.     hf->lsize = (int *)share_malloc(hf->levels * sizeof(int));
  119.     hf->spacing = (float *)share_malloc(hf->levels * sizeof(float));
  120.     hf->boundsmax = (float ***)share_malloc(hf->levels * sizeof(float **));
  121.     hf->boundsmin = (float ***)share_malloc(hf->levels * sizeof(float **));
  122.  
  123.     hf->spacing[0] = hf->size -1;
  124.     hf->lsize[0] = (int)hf->spacing[0];
  125.     hf->boundsmax[0] = (float **)share_malloc(hf->lsize[0]*sizeof(float *));
  126.     hf->boundsmin[0] = (float **)share_malloc(hf->lsize[0]*sizeof(float *));
  127.     /*
  128.      * Compute initial bounding boxes
  129.      */
  130.     for (i = 0; i < hf->lsize[0]; i++) {
  131.         hf->boundsmax[0][i]=(float *)share_malloc(hf->lsize[0]*sizeof(float));
  132.         hf->boundsmin[0][i]=(float *)share_malloc(hf->lsize[0]*sizeof(float));
  133.         maxptr = hf->boundsmax[0][i];
  134.         minptr = hf->boundsmin[0][i];
  135.         for (j = 0; j < hf->lsize[0]; j++) {
  136.             *maxptr++ = maxalt(i, j, hf->data) + EPSILON;
  137.             *minptr++ = minalt(i, j, hf->data) - EPSILON;
  138.         }
  139.     }
  140.  
  141.     for (i = 1; i < hf->levels; i++) {
  142.         hf->spacing[i] = hf->spacing[i-1] * hf->iBestSize;
  143.         hf->lsize[i] = (int)hf->spacing[i];
  144.         if ((Float)hf->lsize[i] != hf->spacing[i])
  145.             hf->lsize[i]++;
  146.         hf->boundsmax[i]=(float **)share_malloc(hf->lsize[i]*sizeof(float *));
  147.         hf->boundsmin[i]=(float **)share_malloc(hf->lsize[i]*sizeof(float *));
  148.         for (j = 0; j < hf->lsize[i]; j++) {
  149.             hf->boundsmax[i][j] = (float *)share_malloc(hf->lsize[i] *
  150.                             sizeof(float));
  151.             hf->boundsmin[i][j] = (float *)share_malloc(hf->lsize[i] *
  152.                             sizeof(float));
  153.         }
  154.         integrate_grid(hf, i);
  155.     }
  156.  
  157.     hf->boundbox[LOW][X] = hf->boundbox[LOW][Y] = 0;
  158.     hf->boundbox[HIGH][X] = hf->boundbox[HIGH][Y] = 1;
  159.     hf->boundbox[LOW][Z] = hf->minz;
  160.     hf->boundbox[HIGH][Z] = hf->maxz;
  161.  
  162.     return hf;
  163. }
  164.  
  165. Methods *
  166. HfMethods()
  167. {
  168.     if (iHfMethods == (Methods *)NULL) {
  169.         iHfMethods = MethodsCreate();
  170.         iHfMethods->create = (GeomCreateFunc *)HfCreate;
  171.         iHfMethods->methods = HfMethods;
  172.         iHfMethods->name = HfName;
  173.         iHfMethods->intersect = HfIntersect;
  174.         iHfMethods->normal = HfNormal;
  175.         iHfMethods->uv = HfUV;
  176.         iHfMethods->bounds = HfBounds;
  177.         iHfMethods->stats = HfStats;
  178.         iHfMethods->checkbounds = TRUE;
  179.         iHfMethods->closed = FALSE;
  180.     }
  181.     return iHfMethods;
  182. }
  183.  
  184. /*
  185.  * Intersect ray with height field.
  186.  */
  187. int
  188. HfIntersect(hf, ray, mindist, maxdist)
  189. Hf *hf;
  190. Ray *ray;
  191. Float mindist, *maxdist;
  192. {
  193.     Vector hitpos;
  194.     Float offset;
  195.     Trav2D trav;
  196.  
  197.     HFTests++;
  198.  
  199.     /*
  200.      * Find where we hit the hf cube.
  201.      */
  202.     VecAddScaled(ray->pos, mindist, ray->dir, &hitpos);
  203.     if (OutOfBounds(&hitpos, hf->boundbox)) {
  204.         offset = *maxdist;
  205.         if (!BoundsIntersect(ray, hf->boundbox, mindist, &offset))
  206.             return FALSE;
  207.         hitpos.x = ray->pos.x + ray->dir.x * offset;
  208.         hitpos.y = ray->pos.y + ray->dir.y * offset;
  209.         hitpos.z = ray->pos.z + ray->dir.z * offset;
  210.     } else
  211.         hitpos = ray->pos;
  212.     /*
  213.      * Find out in which cell "hitpoint" is.
  214.      */
  215.     if (equal(hitpos.x, 1.))
  216.         hitpos.x -= EPSILON;
  217.     if (equal(hitpos.y, 1.))
  218.         hitpos.y -= EPSILON;
  219.  
  220.     if (ray->dir.x < 0.) {
  221.         trav.stepX = trav.outX = -1;
  222.         trav.tDX = -1. / (ray->dir.x * hf->spacing[hf->levels -1]);
  223.     } else if (ray->dir.x > 0.) {
  224.         trav.stepX = 1;
  225.         trav.outX = hf->lsize[hf->levels -1];
  226.         /*
  227.          * (1./size) / ray
  228.          */
  229.         trav.tDX = 1. / (ray->dir.x * hf->spacing[hf->levels -1]);
  230.     }
  231.  
  232.     if (ray->dir.y < 0.) {
  233.         trav.stepY = trav.outY = -1;
  234.         trav.tDY = -1. / (ray->dir.y * hf->spacing[hf->levels -1]);
  235.     } else if (ray->dir.y > 0.) {
  236.         trav.stepY = 1;
  237.         trav.outY = hf->lsize[hf->levels -1];
  238.         trav.tDY = 1. / (ray->dir.y * hf->spacing[hf->levels -1]);
  239.     }
  240.  
  241.     trav.pDX.x = ray->dir.x * trav.tDX;
  242.     trav.pDX.y = ray->dir.y * trav.tDX;
  243.     trav.pDX.z = ray->dir.z * trav.tDX;
  244.     trav.pDY.x = ray->dir.x * trav.tDY;
  245.     trav.pDY.y = ray->dir.y * trav.tDY;
  246.     trav.pDY.z = ray->dir.z * trav.tDY;
  247.  
  248.     trav.cp = hitpos;
  249.     trav.minz = hf->minz;
  250.     trav.maxz = hf->maxz;
  251.     if (DDA2D(hf, &ray->pos, &ray->dir, hf->levels -1, &trav, maxdist)) {
  252.         HFHits++;
  253.         return TRUE;
  254.     }
  255.     return FALSE;
  256. }
  257.  
  258. /*
  259.  * Traverse the grid using a modified DDA algorithm.  If the extent of
  260.  * the ray over a cell intersects the bounding volume defined by the
  261.  * four corners of the cell, either recurse or perform ray/surface
  262.  * intersection test.
  263.  */
  264. static int
  265. DDA2D(hf, pos, ray, level, trav, maxdist)
  266. Hf *hf;
  267. Vector *pos, *ray;
  268. int level;
  269. Trav2D *trav;
  270. Float *maxdist;
  271. {
  272.     int x, y, size, posZ;
  273.     float **boundsmin, **boundsmax, spacing;
  274.     Float tX, tY;
  275.     Trav2D newtrav;
  276.     Vector nxp, nyp;
  277.  
  278.     size = hf->lsize[level];
  279.     spacing = hf->spacing[level];
  280.  
  281.     posZ = (ray->z > 0.);
  282.  
  283.     x = trav->cp.x * hf->spacing[level];
  284.     if (x == size)
  285.         x--;
  286.     y = trav->cp.y * hf->spacing[level];
  287.     if (y == size)
  288.         y--;
  289.     boundsmax = hf->boundsmax[level];
  290.     boundsmin = hf->boundsmin[level];
  291.  
  292.     if (trav->outX > size) trav->outX = size;
  293.     if (trav->outY > size) trav->outY = size;
  294.     if (trav->outX < 0) trav->outX = -1;
  295.     if (trav->outY < 0) trav->outY = -1;
  296.  
  297.     if (ray->x < 0.) {
  298.         tX = (x /spacing - trav->cp.x) / ray->x;
  299.     } else if (ray->x > 0.)
  300.         tX = ((x+1)/spacing - trav->cp.x) / ray->x;
  301.     else
  302.         tX = FAR_AWAY;
  303.     if (ray->y < 0.) {
  304.         tY = (y /spacing - trav->cp.y) / ray->y;
  305.     } else if (ray->y > 0.)
  306.         tY = ((y+1)/spacing - trav->cp.y) / ray->y;
  307.     else
  308.         tY = FAR_AWAY;
  309.  
  310.     nxp.x = trav->cp.x + tX * ray->x;
  311.     nxp.y = trav->cp.y + tX * ray->y;
  312.     nxp.z = trav->cp.z + tX * ray->z;
  313.  
  314.     nyp.x = trav->cp.x + tY * ray->x;
  315.     nyp.y = trav->cp.y + tY * ray->y;
  316.     nyp.z = trav->cp.z + tY * ray->z;
  317.  
  318.     do {
  319.         if (tX < tY) {
  320.             if ((posZ && trav->cp.z <= boundsmax[y][x] &&
  321.                  nxp.z >= boundsmin[y][x]) ||
  322.                 (!posZ && trav->cp.z >= boundsmin[y][x] &&
  323.                  nxp.z <= boundsmax[y][x])) {
  324.                 if (level) {
  325.                     /*
  326.                      * Recurse -- compute constants
  327.                      * needed for next level.
  328.                      * Nicely enough, this just
  329.                      * involves a few multiplications.
  330.                      */
  331.                     newtrav = *trav;
  332.                     newtrav.tDX *= hf->iBestSize;
  333.                     newtrav.tDY *= hf->iBestSize;
  334.                     newtrav.maxz = boundsmax[y][x];
  335.                     newtrav.minz = boundsmin[y][x];
  336.                     if (ray->x < 0.)
  337.                         newtrav.outX=hf->BestSize*x-1;
  338.                     else
  339.                         newtrav.outX=hf->BestSize*(x+1);
  340.                     if (ray->y < 0.)
  341.                         newtrav.outY=hf->BestSize*y-1;
  342.                     else
  343.                         newtrav.outY=hf->BestSize*(y+1);
  344.                     newtrav.pDX.x *= hf->iBestSize;
  345.                     newtrav.pDX.y *= hf->iBestSize;
  346.                     newtrav.pDX.z *= hf->iBestSize;
  347.                     newtrav.pDY.x *= hf->iBestSize;
  348.                     newtrav.pDY.y *= hf->iBestSize;
  349.                     newtrav.pDY.z *= hf->iBestSize;
  350.                     if (DDA2D(hf,pos,ray,level-1,
  351.                         &newtrav, maxdist))
  352.                         return TRUE;
  353.                 } else if (CheckCell(x,y,hf,ray,pos,maxdist))
  354.                     return TRUE;
  355.             }
  356.             x += trav->stepX;        /* Move in X */
  357.             if (*maxdist < tX || x == trav->outX)
  358.                 /* If outside, quit */
  359.                 return FALSE;
  360.             tX += trav->tDX;    /* Update position on ray */
  361.             trav->cp = nxp;        /* cur pos gets next pos */
  362.             nxp.x += trav->pDX.x;    /* Compute next pos */
  363.             nxp.y += trav->pDX.y;
  364.             nxp.z += trav->pDX.z;
  365.         } else {
  366.             if ((posZ && trav->cp.z <= boundsmax[y][x] &&
  367.                  nyp.z >= boundsmin[y][x]) ||
  368.                 (!posZ && trav->cp.z >= boundsmin[y][x] &&
  369.                  nyp.z <= boundsmax[y][x])) {
  370.                 if (level) {
  371.                     /* Recurse */
  372.                     newtrav = *trav;
  373.                     newtrav.tDX *= hf->iBestSize;
  374.                     newtrav.tDY *= hf->iBestSize;
  375.                     newtrav.maxz = boundsmax[y][x];
  376.                     newtrav.minz = boundsmin[y][x];
  377.                     if (ray->x < 0.)
  378.                         newtrav.outX=hf->BestSize*x-1;
  379.                     else
  380.                         newtrav.outX=hf->BestSize*(x+1);
  381.                     if (ray->y < 0.)
  382.                         newtrav.outY=hf->BestSize*y-1;
  383.                     else
  384.                         newtrav.outY=hf->BestSize*(y+1);
  385.                     newtrav.pDX.x *= hf->iBestSize;
  386.                     newtrav.pDX.y *= hf->iBestSize;
  387.                     newtrav.pDX.z *= hf->iBestSize;
  388.                     newtrav.pDY.x *= hf->iBestSize;
  389.                     newtrav.pDY.y *= hf->iBestSize;
  390.                     newtrav.pDY.z *= hf->iBestSize;
  391.                     if (DDA2D(hf,pos,ray,level-1,
  392.                             &newtrav, maxdist))
  393.                         return TRUE;
  394.                 } else if (CheckCell(x,y,hf,ray,pos,maxdist))
  395.                     return TRUE;
  396.             }
  397.             y += trav->stepY;
  398.             if (*maxdist < tY || y == trav->outY)
  399.                 return FALSE;
  400.             tY += trav->tDY;
  401.             trav->cp = nyp;
  402.             nyp.x += trav->pDY.x;
  403.             nyp.y += trav->pDY.y;
  404.             nyp.z += trav->pDY.z;
  405.         }
  406.     } while ((trav->cp.x <= 1. && trav->cp.y <= 1.) &&
  407.          ((posZ && trav->cp.z <= trav->maxz) ||
  408.          (!posZ && trav->cp.z >= trav->minz)));
  409.  
  410.         /*
  411.          * while ((we're inside the horizontal bounding box)
  412.          *        (usually caught by outX & outY, but
  413.          *         it's possible to go "too far" due to
  414.          *         the fact that our levels of grids do
  415.          *         not "nest" exactly if gridsize%BestSize != 0)
  416.          *      and
  417.          *      ((if ray->z is positive and we haven't gone through
  418.          *       the upper bounding plane) or
  419.          *      (if ray->z is negative and we haven't gone through
  420.          *       the lower bounding plane)));
  421.          */
  422.  
  423.     return FALSE;
  424. }
  425.  
  426. /*
  427.  * Check for ray/cell intersection
  428.  */
  429. static int
  430. CheckCell(x, y, hf, ray, pos, maxdist)
  431. int x, y;
  432. Hf *hf;
  433. Vector *ray, *pos;
  434. Float *maxdist;
  435. {
  436.     hfTri *tri1, *tri2;
  437.     Float d1, d2;
  438.  
  439.     d1 = d2 = FAR_AWAY;
  440.  
  441.     if (tri1 = CreateHfTriangle(hf, x, y, x+1, y, x, y+1, TRI1))
  442.         d1 = intHftri(ray, pos, tri1);
  443.     if (tri2 = CreateHfTriangle(hf, x+1, y, x+1, y+1, x, y+1, TRI2))
  444.         d2 = intHftri(ray, pos, tri2);
  445.  
  446.     if (d1 == FAR_AWAY && d2 == FAR_AWAY)
  447.         return FALSE;
  448.  
  449.     if (d1 < d2) {
  450.         if (d1 < *maxdist) {
  451.             hf->hittri = *tri1;
  452.             *maxdist = d1;
  453.             return TRUE;
  454.         }
  455.         return FALSE;
  456.     }
  457.  
  458.     if (d2 < *maxdist) {
  459.         hf->hittri = *tri2;
  460.         *maxdist = d2;
  461.         return TRUE;
  462.     }
  463.     return FALSE;
  464. }
  465.  
  466. static hfTri *
  467. CreateHfTriangle(hf, x1, y1, x2, y2, x3, y3, which)
  468. Hf *hf;
  469. int x1, y1, x2, y2, x3, y3, which;
  470. {
  471.     hfTri *tri;
  472.     Float xid, yid;
  473.     Vector tmp1, tmp2;
  474.  
  475.     /*
  476.      * Don't use triangles with "unset" vertices.
  477.      */
  478.     if (hf->data[y1][x1] == HF_UNSET ||
  479.         hf->data[y2][x2] == HF_UNSET ||
  480.         hf->data[y3][x3] == HF_UNSET)
  481.         return (hfTri *)0;
  482.  
  483.     xid = (Float)x1 / (Float)(hf->size -1);
  484.     yid = (Float)y1 / (Float)(hf->size -1);
  485.  
  486.     if ((tri = GetQueuedTri(hf, xid, yid, which)) != (hfTri *)0)
  487.         return tri;
  488.  
  489.     tri = (hfTri *)Malloc(sizeof(hfTri));
  490.  
  491.     tri->type = which;
  492.     tri->v1.x = xid;
  493.     tri->v1.y = yid;
  494.     tri->v1.z = hf->data[y1][x1];
  495.     tri->v2.x = (Float)x2 / (Float)(hf->size-1);
  496.     tri->v2.y = (Float)y2 / (Float)(hf->size-1);
  497.     tri->v2.z = hf->data[y2][x2];
  498.     tri->v3.x = (Float)x3 / (Float)(hf->size-1);
  499.     tri->v3.y = (Float)y3 / (Float)(hf->size-1);
  500.     tri->v3.z = hf->data[y3][x3];
  501.  
  502.     tmp1.x = tri->v2.x - tri->v1.x;
  503.     tmp1.y = tri->v2.y - tri->v1.y;
  504.     tmp1.z = tri->v2.z - tri->v1.z;
  505.     tmp2.x = tri->v3.x - tri->v1.x;
  506.     tmp2.y = tri->v3.y - tri->v1.y;
  507.     tmp2.z = tri->v3.z - tri->v1.z;
  508.  
  509.     (void)VecNormCross(&tmp1, &tmp2, &tri->norm);
  510.  
  511.     tri->d = -dotp(&tri->v1, &tri->norm);
  512.  
  513.     QueueTri(hf, tri);
  514.     return tri;
  515. }
  516.  
  517. /*
  518.  * Intersect ray with right isoscoles triangle, the hypotenuse of which
  519.  * has slope of -1.
  520.  */
  521. static Float
  522. intHftri(ray, pos, tri)
  523. hfTri *tri;
  524. Vector *pos, *ray;
  525. {
  526.     Float u, v, dist, xpos, ypos;
  527.  
  528.     u = dotp(&tri->norm, pos) + tri->d;
  529.     v = dotp(&tri->norm, ray);
  530.  
  531.     if ((u <= 0. || v > -EPSILON) && (u >= 0. && v < EPSILON))
  532.         return FAR_AWAY;
  533.  
  534.     dist = -u / v;
  535.  
  536.     if (dist < EPSILON)
  537.         return FAR_AWAY;
  538.  
  539.     xpos = pos->x + dist*ray->x;
  540.     ypos = pos->y + dist*ray->y;
  541.  
  542.     if (tri->type == TRI1 && xpos >= tri->v1.x && ypos >= tri->v1.y &&
  543.         xpos + ypos <= tri->v2.x + tri->v2.y)
  544.         return dist;
  545.     if (tri->type == TRI2 && xpos <= tri->v2.x && ypos <= tri->v2.y &&
  546.         xpos + ypos >= tri->v1.x + tri->v1.y)
  547.         return dist;
  548.     return FAR_AWAY;
  549. }
  550.  
  551. /*
  552.  * Compute normal to height field.
  553.  */
  554. /*ARGSUSED*/
  555. int
  556. HfNormal(hf, pos, nrm, gnrm)
  557. Hf *hf;
  558. Vector *pos, *nrm, *gnrm;
  559. {
  560.     *gnrm = *nrm = hf->hittri.norm;
  561.     return FALSE;
  562. }
  563.  
  564. /*ARGSUSED*/
  565. void
  566. HfUV(hf, pos, norm, uv, dpdu, dpdv)
  567. Hf *hf;
  568. Vector *pos, *norm, *dpdu, *dpdv;
  569. Vec2d *uv;
  570. {
  571.     uv->u = pos->x;
  572.     uv->v = pos->y;
  573.     if (dpdu) {
  574.         dpdu->x = 1.;
  575.         dpdu->y = dpdv->z = 0.;
  576.         dpdv->x = dpdv->z = 0.;
  577.         dpdv->y = 1.;
  578.     }
  579. }
  580.  
  581. /*
  582.  * Compute heightfield bounding box.
  583.  */
  584. void
  585. HfBounds(hf, bounds)
  586. Hf *hf;
  587. Float bounds[2][3];
  588. {
  589.     /*
  590.      * By default, height fields are centered at (0.5, 0.5, 0.)
  591.      */
  592.     bounds[LOW][X] = bounds[LOW][Y] = 0;
  593.     bounds[HIGH][X] = bounds[HIGH][Y] = 1;
  594.     bounds[LOW][Z] = hf->minz;
  595.     bounds[HIGH][Z] = hf->maxz;
  596. }
  597.  
  598. /*
  599.  * Build min/max altitude value arrays for the given grid level.
  600.  */
  601. static void
  602. integrate_grid(hf, level)
  603. Hf *hf;
  604. int level;
  605. {
  606.     int i, j, k, l, ii, ji;
  607.     float max_alt, min_alt;
  608.     float **maxinto, **mininto, **frommax, **frommin, *minptr, *maxptr;
  609.     int insize, fromsize, fact;
  610.  
  611.     maxinto = hf->boundsmax[level];
  612.     mininto = hf->boundsmin[level];
  613.     insize = hf->lsize[level];
  614.     frommax = hf->boundsmax[level-1];
  615.     frommin = hf->boundsmin[level-1];
  616.     fact = hf->BestSize;
  617.     fromsize = hf->lsize[level-1];
  618.  
  619.     ii = 0;
  620.  
  621.     for (i = 0; i < insize; i++) {
  622.         ji = 0;
  623.         for (j = 0; j < insize; j++) {
  624.             max_alt = HF_UNSET;
  625.             min_alt = -HF_UNSET;
  626.             for (k = 0; k <= fact; k++) {
  627.                 if (ii+k >= fromsize)
  628.                     continue;
  629.                 maxptr = &frommax[ii+k][ji];
  630.                 minptr = &frommin[ii+k][ji];
  631.                 for (l = 0; l <= fact; l++,maxptr++,minptr++) {
  632.                     if (ji+l >= fromsize)
  633.                         continue;
  634.                     if (*maxptr > max_alt)
  635.                         max_alt = *maxptr;
  636.                     if (*minptr < min_alt)
  637.                         min_alt = *minptr;
  638.                 }
  639.             }
  640.             maxinto[i][j] = max_alt + EPSILON;
  641.             mininto[i][j] = min_alt - EPSILON;
  642.             ji += fact;
  643.         }
  644.         ii += fact;
  645.     }
  646. }
  647.  
  648. /*
  649.  * Place the given triangle in the triangle cache.
  650.  */
  651. static void
  652. QueueTri(hf, tri)
  653. Hf *hf;
  654. hfTri *tri;
  655. {
  656.     if (hf->q[hf->qtail])        /* Free old triangle data */
  657.         Free((voidstar)hf->q[hf->qtail]);
  658.     hf->q[hf->qtail] = tri;        /* Put on tail */
  659.     hf->qtail = (hf->qtail + 1) % hf->qsize;    /* Increment tail */
  660. }
  661.  
  662. /*
  663.  * Search list of cached trianges to see if this triangle has been
  664.  * cached.  If so, return a pointer to it.  If not, return null pointer.
  665.  */
  666. static hfTri *
  667. GetQueuedTri(hf, x, y, flag)
  668. Hf *hf;
  669. Float x, y;
  670. int flag;
  671. {
  672.     register int i;
  673.     register hfTri **tmp;
  674.  
  675.     for (i = 0, tmp = hf->q; i < hf->qsize; i++, tmp++) {
  676.         if (*tmp && (*tmp)->v1.x == x && (*tmp)->v1.y == y &&
  677.             (*tmp)->type == flag)
  678.             return *tmp;    /* vertices & flag match, return it */
  679.     }
  680.  
  681.     return (hfTri *)0;
  682. }
  683.  
  684. /*
  685.  * Return maximum height of cell indexed by y,x.  This could be done
  686.  * as a macro, but many C compliers will choke on it.
  687.  */
  688. static float
  689. minalt(y,x,data)
  690. int x, y;
  691. float **data;
  692. {
  693.     float  min_alt;
  694.  
  695.     min_alt = min(data[y][x], data[y+1][x]);
  696.     min_alt = min(min_alt, data[y][x+1]);
  697.     min_alt = min(min_alt, data[y+1][x+1]);
  698.     return min_alt;
  699. }
  700.  
  701. /*
  702.  * Return maximum cell height, as above.
  703.  */
  704. static float
  705. maxalt(y,x,data)
  706. int x, y;
  707. float **data;
  708. {
  709.     float  max_alt;
  710.  
  711.     max_alt = max(data[y][x], data[y+1][x]);
  712.     max_alt = max(max_alt, data[y][x+1]);
  713.     max_alt = max(max_alt, data[y+1][x+1]);
  714.     return max_alt;
  715. }
  716.  
  717. char *
  718. HfName()
  719. {
  720.     return hfName;
  721. }
  722.  
  723. void
  724. HfStats(tests, hits)
  725. unsigned long *tests, *hits;
  726. {
  727.     *tests = HFTests;
  728.     *hits = HFHits;
  729. }
  730.  
  731. void
  732. HfMethodRegister(meth)
  733. UserMethodType meth;
  734. {
  735.     if (iHfMethods)
  736.         iHfMethods->user = meth;
  737. }
  738.