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