home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Graphics / Graphics.zip / povsrc31.zip / hfield.c < prev    next >
C/C++ Source or Header  |  2000-05-27  |  49KB  |  2,301 lines
  1. /****************************************************************************
  2. *                   hfield.c
  3. *
  4. *    This file implements the height field shape primitive.  The shape is
  5. *    implemented as a collection of triangles which are calculated as
  6. *    needed.  The basic intersection routine first computes the rays
  7. *    intersection with the box marking the limits of the shape, then
  8. *    follows the line from one intersection point to the other, testing the
  9. *    two triangles which form the pixel for an intersection with the ray at
  10. *    each step.
  11. *
  12. *    height field added by Doug Muir with lots of advice and support
  13. *    from David Buck and Drew Wells.
  14. *
  15. *  from Persistence of Vision(tm) Ray Tracer
  16. *  Copyright 1996,1999 Persistence of Vision Team
  17. *---------------------------------------------------------------------------
  18. *  NOTICE: This source code file is provided so that users may experiment
  19. *  with enhancements to POV-Ray and to port the software to platforms other
  20. *  than those supported by the POV-Ray Team.  There are strict rules under
  21. *  which you are permitted to use this file.  The rules are in the file
  22. *  named POVLEGAL.DOC which should be distributed with this file.
  23. *  If POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  24. *  Team Coordinator by email to team-coord@povray.org or visit us on the web at
  25. *  http://www.povray.org. The latest version of POV-Ray may be found at this site.
  26. *
  27. * This program is based on the popular DKB raytracer version 2.12.
  28. * DKBTrace was originally written by David K. Buck.
  29. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  30. *
  31. *****************************************************************************/
  32.  
  33. /****************************************************************************
  34. *
  35. *  Explanation:
  36. *
  37. *    -
  38. *
  39. *  Syntax:
  40. *
  41. *  ---
  42. *
  43. *  Aug 1994 : Merged functions for CSG height fields into functions for
  44. *             non-CSG height fiels. Moved all height field map related
  45. *             data into one data structure. Fixed memory problems. [DB]
  46. *
  47. *  Feb 1995 : Major rewrite of the height field intersection tests. [DB]
  48. *
  49. *****************************************************************************/
  50.  
  51. #include "frame.h"
  52. #include "povray.h"
  53. #include "vector.h"
  54. #include "povproto.h"
  55. #include "hfield.h"
  56. #include "matrices.h"
  57. #include "objects.h"
  58.  
  59.  
  60.  
  61. /*****************************************************************************
  62. * Local preprocessor defines
  63. ******************************************************************************/
  64.  
  65. #define sign(x) (((x) >= 0.0) ? 1.0 : -1.0)
  66.  
  67. #define Get_Height(x, z, HField) ((DBL)(HField)->Data->Map[(z)][(x)])
  68.  
  69. /* Small offest. */
  70.  
  71. #define HFIELD_OFFSET 0.001
  72.  
  73.  
  74.  
  75. /*****************************************************************************
  76. * Static functions
  77. ******************************************************************************/
  78.  
  79. static DBL normalize (VECTOR A, VECTOR B);
  80.  
  81. static DBL stretch (DBL x);
  82.  
  83. static void smooth_height_field (HFIELD *HField, int xsize, int zsize);
  84.  
  85. static int intersect_pixel (int x,int z,RAY *Ray,HFIELD *HField,DBL height1,DBL height2);
  86.  
  87. static int add_single_normal (HF_VAL **data, int xsize, int zsize,
  88.   int x0, int z0,int x1, int z1,int x2, int z2,VECTOR N);
  89.  
  90. static int  All_HField_Intersections (OBJECT *Object,RAY *Ray,ISTACK *Depth_Stack);
  91. static int  Inside_HField (VECTOR IPoint,OBJECT *Object);
  92. static void HField_Normal (VECTOR Result,OBJECT *Object,INTERSECTION *Inter);
  93. static void Translate_HField (OBJECT *Object,VECTOR Vector, TRANSFORM *Trans);
  94. static void Rotate_HField (OBJECT *Object,VECTOR Vector, TRANSFORM *Trans);
  95. static void Scale_HField (OBJECT *Object,VECTOR Vector, TRANSFORM *Trans);
  96. static void Invert_HField (OBJECT *Object);
  97. static void Transform_HField (OBJECT *Object,TRANSFORM *Trans);
  98. static HFIELD *Copy_HField (OBJECT *Object);
  99. static void Destroy_HField (OBJECT *Object);
  100.  
  101. static int dda_traversal (RAY *Ray, HFIELD *HField, VECTOR Start, HFIELD_BLOCK *Block);
  102. static int block_traversal (RAY *Ray, HFIELD *HField, VECTOR Start);
  103. static void build_hfield_blocks (HFIELD *HField);
  104.  
  105.  
  106.  
  107. /*****************************************************************************
  108. * Local variables
  109. ******************************************************************************/
  110.  
  111. METHODS HField_Methods =
  112. {
  113.   All_HField_Intersections,
  114.   Inside_HField, HField_Normal, Default_UVCoord,
  115.   (COPY_METHOD)Copy_HField, Translate_HField, Rotate_HField,
  116.   Scale_HField, Transform_HField, Invert_HField, Destroy_HField
  117. };
  118.  
  119. static ISTACK *HField_Stack;
  120. static RAY *RRay;
  121. static DBL mindist, maxdist;
  122.  
  123.  
  124. /*****************************************************************************
  125. *
  126. * FUNCTION
  127. *
  128. *   All_HField_Intersections
  129. *
  130. * INPUT
  131. *   
  132. * OUTPUT
  133. *   
  134. * RETURNS
  135. *   
  136. * AUTHOR
  137. *
  138. *   Doug Muir, David Buck, Drew Wells
  139. *   
  140. * DESCRIPTION
  141. *
  142. *   -
  143. *
  144. * CHANGES
  145. *
  146. *   Feb 1995 : Modified to work with new intersection functions. [DB]
  147. *
  148. ******************************************************************************/
  149.  
  150. static int All_HField_Intersections(OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack)
  151. {
  152.   int Side1, Side2;
  153.   VECTOR Start;
  154.   RAY Temp_Ray;
  155.   DBL depth1, depth2;
  156.   HFIELD *HField = (HFIELD *) Object;
  157.  
  158.   Increase_Counter(stats[Ray_HField_Tests]);
  159.  
  160.   MInvTransPoint(Temp_Ray.Initial, Ray->Initial, HField->Trans);
  161.   MInvTransDirection(Temp_Ray.Direction, Ray->Direction, HField->Trans);
  162.  
  163. #ifdef HFIELD_EXTRA_STATS
  164.   Increase_Counter(stats[Ray_HField_Box_Tests]);
  165. #endif
  166.  
  167.   if (!Intersect_Box(&Temp_Ray,HField->bounding_box,&depth1,&depth2,&Side1,&Side2))
  168.   {
  169.     return(FALSE);
  170.   }
  171.  
  172. #ifdef HFIELD_EXTRA_STATS
  173.   Increase_Counter(stats[Ray_HField_Box_Tests_Succeeded]);
  174. #endif
  175.  
  176.   HField->Data->cache_pos = 0;
  177.  
  178.   if (depth1 < EPSILON)
  179.   {
  180.     depth1 = EPSILON;
  181.  
  182.     if (depth1 > depth2)
  183.     {
  184.       return(FALSE);
  185.     }
  186.   }
  187.  
  188.   VEvaluateRay(Start, Temp_Ray.Initial, depth1, Temp_Ray.Direction);
  189.  
  190.   mindist = depth1;
  191.   maxdist = depth2;
  192.  
  193.   HField_Stack = Depth_Stack;
  194.  
  195.   RRay = Ray;
  196.  
  197.   if (block_traversal(&Temp_Ray, HField, Start))
  198.   {
  199.     Increase_Counter(stats[Ray_HField_Tests_Succeeded]);
  200.  
  201.     return(TRUE);
  202.   }
  203.   else
  204.   {
  205.     return(FALSE);
  206.   }
  207. }
  208.  
  209.  
  210.  
  211. /*****************************************************************************
  212. *
  213. * FUNCTION
  214. *
  215. *   Inside_HField
  216. *
  217. * INPUT
  218. *   
  219. * OUTPUT
  220. *   
  221. * RETURNS
  222. *   
  223. * AUTHOR
  224. *
  225. *   Doug Muir, David Buck, Drew Wells
  226. *   
  227. * DESCRIPTION
  228. *
  229. *   -
  230. *
  231. * CHANGES
  232. *
  233. *   -
  234. *
  235. ******************************************************************************/
  236.  
  237. static int Inside_HField (VECTOR IPoint, OBJECT *Object)
  238. {
  239.   HFIELD *HField = (HFIELD *) Object;
  240.   int px, pz;
  241.   DBL x,z,y1,y2,y3,water, dot1, dot2;
  242.   VECTOR Local_Origin, H_Normal, Test;
  243.  
  244.   MInvTransPoint(Test, IPoint, HField->Trans);
  245.  
  246.   water = HField->bounding_box->bounds[0][Y];
  247.  
  248.   if ((Test[X] < 0.0) || (Test[X] >= HField->bounding_box->bounds[1][X]) ||
  249.       (Test[Z] < 0.0) || (Test[Z] >= HField->bounding_box->bounds[1][Z]))
  250.   {
  251.     return(Test_Flag(HField, INVERTED_FLAG));
  252.   }
  253.  
  254.   if (Test[Y] >= HField->bounding_box->bounds[1][Y])
  255.   {
  256.     return(Test_Flag(HField, INVERTED_FLAG));
  257.   }
  258.  
  259.   if (Test[Y] < water)
  260.   {
  261.     return(!Test_Flag(HField, INVERTED_FLAG));
  262.   }
  263.  
  264.   px = (int)Test[X];
  265.   pz = (int)Test[Z];
  266.  
  267.   x = Test[X] - (DBL)px;
  268.   z = Test[Z] - (DBL)pz;
  269.  
  270.   if ((x+z) < 1.0)
  271.   {
  272.     y1 = max(Get_Height(px,   pz,   HField), water);
  273.     y2 = max(Get_Height(px+1, pz,   HField), water);
  274.     y3 = max(Get_Height(px,   pz+1, HField), water);
  275.  
  276.     Make_Vector(Local_Origin,(DBL)px,y1,(DBL)pz);
  277.  
  278.     Make_Vector(H_Normal, y1-y2, 1.0, y1-y3);
  279.   }
  280.   else
  281.   {
  282.     px = (int)ceil(Test[X]);
  283.     pz = (int)ceil(Test[Z]);
  284.  
  285.     y1 = max(Get_Height(px,   pz,   HField), water);
  286.     y2 = max(Get_Height(px-1, pz,   HField), water);
  287.     y3 = max(Get_Height(px,   pz-1, HField), water);
  288.  
  289.     Make_Vector(Local_Origin,(DBL)px,y1,(DBL)pz);
  290.  
  291.     Make_Vector(H_Normal, y2-y1, 1.0, y3-y1);
  292.   }
  293.  
  294.   VDot(dot1, Test, H_Normal);
  295.   VDot(dot2, Local_Origin, H_Normal);
  296.  
  297.   if (dot1 < dot2)
  298.   {
  299.     return(!Test_Flag(HField, INVERTED_FLAG));
  300.   }
  301.  
  302.   return(Test_Flag(HField, INVERTED_FLAG));
  303. }
  304.  
  305.  
  306.  
  307. /*****************************************************************************
  308. *
  309. * FUNCTION
  310. *
  311. *   HField_Normal
  312. *
  313. * INPUT
  314. *   
  315. * OUTPUT
  316. *   
  317. * RETURNS
  318. *   
  319. * AUTHOR
  320. *
  321. *   Doug Muir, David Buck, Drew Wells
  322. *   
  323. * DESCRIPTION
  324. *
  325. *   -
  326. *
  327. * CHANGES
  328. *
  329. *   June 1999: Modified by Lummox JR to allow different types of smoothing
  330. *              via a new smooth_type variable. Type 0 is no smoothing, 1 is
  331. *              standard, 2 is biquadric, and 3 is biquadric with the
  332. *              stretch function applied.
  333. *
  334. *
  335. ******************************************************************************/
  336.  
  337. static void HField_Normal (VECTOR Result, OBJECT *Object, INTERSECTION *Inter)
  338. {
  339.   HFIELD *HField = (HFIELD *) Object;
  340.   int px,pz, i;
  341.   DBL x,z,y1,y2,y3,u,v,biquad[9];
  342.   VECTOR Local_Origin;
  343.   VECTOR n[5];
  344.  
  345.   MInvTransPoint(Local_Origin, Inter->IPoint, HField->Trans);
  346.  
  347.   for (i = 0; i < HField->Data->cache_pos; i++)
  348.   {
  349.     if ((Local_Origin[X] == HField->Data->Normal_Cache[i].fx) &&
  350.        (Local_Origin[Z] == HField->Data->Normal_Cache[i].fz))
  351.     {
  352.       Assign_Vector(Result,HField->Data->Normal_Cache[i].normal);
  353.  
  354.       MTransNormal(Result,Result,HField->Trans);
  355.  
  356.       VNormalize(Result,Result);
  357.  
  358.       return;
  359.     }
  360.   }
  361.  
  362.   px = (int)Local_Origin[X];
  363.   pz = (int)Local_Origin[Z];
  364.  
  365.   /* NK crash bugfix here */
  366.   if (px>HField->Data->max_x) px=HField->Data->max_x;
  367.   if (pz>HField->Data->max_z) pz=HField->Data->max_z;
  368.   /* NK ---- */
  369.  
  370.   x = Local_Origin[X] - (DBL)px;
  371.   z = Local_Origin[Z] - (DBL)pz;
  372.  
  373.   if (Test_Flag(HField, SMOOTHED_FLAG))
  374.   {
  375.     n[0][X] = HField->Data->Normals[pz][px][0];
  376.     n[0][Y] = HField->Data->Normals[pz][px][1];
  377.     n[0][Z] = HField->Data->Normals[pz][px][2];
  378.     n[1][X] = HField->Data->Normals[pz][px+1][0];
  379.     n[1][Y] = HField->Data->Normals[pz][px+1][1];
  380.     n[1][Z] = HField->Data->Normals[pz][px+1][2];
  381.     n[2][X] = HField->Data->Normals[pz+1][px][0];
  382.     n[2][Y] = HField->Data->Normals[pz+1][px][1];
  383.     n[2][Z] = HField->Data->Normals[pz+1][px][2];
  384.     n[3][X] = HField->Data->Normals[pz+1][px+1][0];
  385.     n[3][Y] = HField->Data->Normals[pz+1][px+1][1];
  386.     n[3][Z] = HField->Data->Normals[pz+1][px+1][2];
  387.  
  388.     if(HField->smooth_type==1 || HField->smooth_type==3) {
  389.       x = stretch(x);
  390.       z = stretch(z);
  391.     }
  392.  
  393.     if(HField->smooth_type<=1)
  394.     {
  395.       u = (1.0 - x);
  396.       v = (1.0 - z);
  397.  
  398.       Result[X] = v*(u*n[0][X] + x*n[1][X]) + z*(u*n[2][X] + x*n[3][X]);
  399.       Result[Y] = v*(u*n[0][Y] + x*n[1][Y]) + z*(u*n[2][Y] + x*n[3][Y]);
  400.       Result[Z] = v*(u*n[0][Z] + x*n[1][Z]) + z*(u*n[2][Z] + x*n[3][Z]);
  401.     }
  402.     else
  403.     {
  404.       /*
  405.         Biquadric interpolation
  406.         added by Lummox JR, June 1999
  407.  
  408.         Assume y=x^2(az^2+bz+c)+x(dz^2+ez+f)+(gz^2+hz+i)
  409.  
  410.         Normal vector is:
  411.             <2x(az^2+bz+c)+(dz^2+ez+f),1,2z(ax^2+dx+g)+(bx^2+ex+h)>
  412.         Given normals at the four corners, find a through h. Assume a=1,
  413.         since it is arbitrary.
  414.        */
  415.       biquad[4]=n[0][X]/n[0][Y];                /* f*/
  416.       biquad[6]=n[0][Z]/n[0][Y];                /* h*/
  417.       biquad[1]=(n[1][X]/n[1][Y]-biquad[4])/2;  /* c*/
  418.       biquad[5]=(n[2][Z]/n[2][Y]-biquad[6])/2;  /* g*/
  419.       biquad[7]=(n[1][Z]/n[1][Y]-biquad[6])/2;  /* b+e*/
  420.       biquad[8]=(n[2][X]/n[2][Y]-biquad[4])/2;  /* d+e*/
  421.       biquad[0]=(n[3][X]/n[3][Y]-biquad[8]-biquad[4])/2-biquad[1]-1;     /* b  (a+b-1)*/
  422.       biquad[2]=(n[3][Z]/n[3][Y]-biquad[7]-biquad[6])/2-biquad[5]-1;     /* d  (a+d-1)*/
  423.       biquad[3]=(biquad[7]+biquad[8]-biquad[0]-biquad[2])/2;
  424.  
  425.       Result[X] = 2.0*x*(z*z+biquad[0]*z+biquad[1])+biquad[2]*z*z+biquad[3]*z+biquad[4];
  426.       Result[Y] = 1;
  427.       Result[Z] = 2.0*z*(x*x+biquad[2]*x+biquad[5])+biquad[0]*x*x+biquad[3]*x+biquad[6];
  428.     }
  429.   }
  430.   else
  431.   {
  432.     if ((x+z) <= 1.0)
  433.     {
  434.       /* Lower triangle. */
  435.  
  436.       y1 = Get_Height(px,   pz,   HField);
  437.       y2 = Get_Height(px+1, pz,   HField);
  438.       y3 = Get_Height(px,   pz+1, HField);
  439.  
  440.       Make_Vector(Result, y1-y2, 1.0, y1-y3);
  441.     }
  442.     else
  443.     {
  444.       /* Upper triangle. */
  445.  
  446.       y1 = Get_Height(px+1, pz+1, HField);
  447.       y2 = Get_Height(px,   pz+1, HField);
  448.       y3 = Get_Height(px+1, pz,   HField);
  449.  
  450.       Make_Vector(Result, y2-y1, 1.0, y3-y1);
  451.     }
  452.   }
  453.  
  454.   MTransNormal(Result, Result, HField->Trans);
  455.  
  456.   VNormalize(Result, Result);
  457. }
  458.  
  459.  
  460.  
  461. /*****************************************************************************
  462. *
  463. * FUNCTION
  464. *
  465. *   stretch
  466. *
  467. * INPUT
  468. *   
  469. * OUTPUT
  470. *   
  471. * RETURNS
  472. *   
  473. * AUTHOR
  474. *
  475. *   Doug Muir, David Buck, Drew Wells
  476. *   
  477. * DESCRIPTION
  478. *
  479. *   -
  480. *
  481. * CHANGES
  482. *
  483. *   -
  484. *
  485. ******************************************************************************/
  486.  
  487. static DBL stretch (DBL x)
  488. {
  489.   if (x <= 0.5)
  490.   {
  491.     x = 2.0 * x * x;
  492.   }
  493.   else
  494.   {
  495.     x = 1.0 - (2.0 * (1.0 - x) * (1.0 - x));
  496.   }
  497.  
  498.   return(x);
  499. }
  500.  
  501.  
  502.  
  503. /*****************************************************************************
  504. *
  505. * FUNCTION
  506. *
  507. *   normalize
  508. *
  509. * INPUT
  510. *   
  511. * OUTPUT
  512. *   
  513. * RETURNS
  514. *   
  515. * AUTHOR
  516. *
  517. *   Doug Muir, David Buck, Drew Wells
  518. *   
  519. * DESCRIPTION
  520. *
  521. *   -
  522. *
  523. * CHANGES
  524. *
  525. *   -
  526. *
  527. ******************************************************************************/
  528.  
  529. static DBL normalize(VECTOR A, VECTOR  B)
  530. {
  531.   VLength(VTemp, B);
  532.  
  533.   if (fabs(VTemp) > EPSILON)
  534.   {
  535.     VInverseScale(A, B, VTemp);
  536.   }
  537.   else
  538.   {
  539.     Make_Vector(A, 0.0, 1.0, 0.0);
  540.   }
  541.  
  542.   return(VTemp);
  543. }
  544.  
  545.  
  546.  
  547. /*****************************************************************************
  548. *
  549. * FUNCTION
  550. *
  551. *   intersect_pixel
  552. *
  553. * INPUT
  554. *   
  555. * OUTPUT
  556. *   
  557. * RETURNS
  558. *   
  559. * AUTHOR
  560. *
  561. *   Doug Muir, David Buck, Drew Wells
  562. *   
  563. * DESCRIPTION
  564. *
  565. *   -
  566. *
  567. * CHANGES
  568. *
  569. *   -
  570. *
  571. ******************************************************************************/
  572.  
  573. static int intersect_pixel(int x, int z, RAY *Ray, HFIELD *HField, DBL height1, DBL  height2)
  574. {
  575.   int Found;
  576.   DBL dot, depth1, depth2;
  577.   DBL s, t, y1, y2, y3, y4;
  578.   DBL min_y2_y3, max_y2_y3;
  579.   DBL max_height, min_height;
  580.   VECTOR P, N, V1;
  581.  
  582. #ifdef HFIELD_EXTRA_STATS
  583.   Increase_Counter(stats[Ray_HField_Cell_Tests]);
  584. #endif
  585.  
  586.   /* NK crash bugfix here */
  587.   if (z>HField->Data->max_z) z = HField->Data->max_z;
  588.   if (x>HField->Data->max_x) x = HField->Data->max_x;
  589.   /* NK ---- */
  590.  
  591.  
  592.   y1 = Get_Height(x,   z,   HField);
  593.   y2 = Get_Height(x+1, z,   HField);
  594.   y3 = Get_Height(x,   z+1, HField);
  595.   y4 = Get_Height(x+1, z+1, HField);
  596.  
  597.   /* Do we hit this cell at all? */
  598.  
  599.   if (y2 < y3)
  600.   {
  601.     min_y2_y3 = y2;
  602.     max_y2_y3 = y3;
  603.   }
  604.   else
  605.   {
  606.     min_y2_y3 = y3;
  607.     max_y2_y3 = y2;
  608.   }
  609.  
  610.   min_height = min(min(y1, y4), min_y2_y3);
  611.   max_height = max(max(y1, y4), max_y2_y3);
  612.  
  613.   if ((max_height < height1) || (min_height > height2))
  614.   {
  615.     return(FALSE);
  616.   }
  617.  
  618. #ifdef HFIELD_EXTRA_STATS
  619.   Increase_Counter(stats[Ray_HField_Cell_Tests_Succeeded]);
  620. #endif
  621.  
  622.   Found = FALSE;
  623.  
  624.   /* Check if we'll hit first triangle. */
  625.  
  626.   min_height = min(y1, min_y2_y3);
  627.   max_height = max(y1, max_y2_y3);
  628.  
  629.   if ((max_height >= height1) && (min_height <= height2))
  630.   {
  631. #ifdef HFIELD_EXTRA_STATS
  632.     Increase_Counter(stats[Ray_HField_Triangle_Tests]);
  633. #endif
  634.  
  635.     /* Set up triangle. */
  636.  
  637.     Make_Vector(P, (DBL)x, y1, (DBL)z);
  638.  
  639.     /*
  640.      * Calculate the normal vector from:
  641.      *
  642.      * N = V2 x V1, with V1 = <1, y2-y1, 0>, V2 = <0, y3-y1, 1>.
  643.      */
  644.  
  645.     Make_Vector(N, y1-y2, 1.0, y1-y3);
  646.  
  647.     /* Now intersect the triangle. */
  648.  
  649.     VDot(dot, N, Ray->Direction);
  650.  
  651.     if ((dot > EPSILON) || (dot < -EPSILON))
  652.     {
  653.       VSub(V1, P, Ray->Initial);
  654.  
  655.       VDot(depth1, N, V1);
  656.  
  657.       depth1 /= dot;
  658.  
  659.       if ((depth1 >= mindist) && (depth1 <= maxdist))
  660.       {
  661.         s = Ray->Initial[X] + depth1 * Ray->Direction[X] - (DBL)x;
  662.         t = Ray->Initial[Z] + depth1 * Ray->Direction[Z] - (DBL)z;
  663.  
  664.         if ((s >= -0.0001) && (t >= -0.0001) && ((s+t) <= 1.0001))
  665.         {
  666. #ifdef HFIELD_EXTRA_STATS
  667.           Increase_Counter(stats[Ray_HField_Triangle_Tests_Succeeded]);
  668. #endif
  669.  
  670.           VEvaluateRay(P, RRay->Initial, depth1, RRay->Direction);
  671.  
  672.           if (Point_In_Clip(P, HField->Clip))
  673.           {
  674.             push_entry(depth1, P, (OBJECT *)HField, HField_Stack);
  675.  
  676.             Found = TRUE;
  677.  
  678.             /* Cache normal. */
  679.  
  680.             if (!Test_Flag(HField, SMOOTHED_FLAG))
  681.             {
  682.               if (HField->Data->cache_pos < HF_CACHE_SIZE)
  683.               {
  684.                 Assign_Vector(HField->Data->Normal_Cache[HField->Data->cache_pos].normal, N);
  685.  
  686.                 HField->Data->Normal_Cache[HField->Data->cache_pos].fx = x + s;
  687.                 HField->Data->Normal_Cache[HField->Data->cache_pos].fz = z + t;
  688.  
  689.                 HField->Data->cache_pos++;
  690.               }
  691.             }
  692.           }
  693.         }
  694.       }
  695.     }
  696.   }
  697.  
  698.   /* Check if we'll hit second triangle. */
  699.  
  700.   min_height = min(y4, min_y2_y3);
  701.   max_height = max(y4, max_y2_y3);
  702.  
  703.   if ((max_height >= height1) && (min_height <= height2))
  704.   {
  705. #ifdef HFIELD_EXTRA_STATS
  706.     Increase_Counter(stats[Ray_HField_Triangle_Tests]);
  707. #endif
  708.  
  709.     /* Set up triangle. */
  710.  
  711.     Make_Vector(P, (DBL)(x+1), y4, (DBL)(z+1));
  712.  
  713.     /*
  714.      * Calculate the normal vector from:
  715.      *
  716.      * N = V2 x V1, with V1 = <-1, y3-y4, 0>, V2 = <0, y2-y4, -1>.
  717.      */
  718.  
  719.     Make_Vector(N, y3-y4, 1.0, y2-y4);
  720.  
  721.     /* Now intersect the triangle. */
  722.  
  723.     VDot(dot, N, Ray->Direction);
  724.  
  725.     if ((dot > EPSILON) || (dot < -EPSILON))
  726.     {
  727.       VSub(V1, P, Ray->Initial);
  728.  
  729.       VDot(depth2, N, V1);
  730.  
  731.       depth2 /= dot;
  732.  
  733.       if ((depth2 >= mindist) && (depth2 <= maxdist))
  734.       {
  735.         s = Ray->Initial[X] + depth2 * Ray->Direction[X] - (DBL)x;
  736.         t = Ray->Initial[Z] + depth2 * Ray->Direction[Z] - (DBL)z;
  737.  
  738.         if ((s <= 1.0001) && (t <= 1.0001) && ((s+t) >= 0.9999))
  739.         {
  740. #ifdef HFIELD_EXTRA_STATS
  741.           Increase_Counter(stats[Ray_HField_Triangle_Tests_Succeeded]);
  742. #endif
  743.  
  744.           VEvaluateRay(P, RRay->Initial, depth2, RRay->Direction);
  745.  
  746.           if (Point_In_Clip(P, HField->Clip))
  747.           {
  748.             push_entry(depth2, P, (OBJECT *)HField, HField_Stack);
  749.  
  750.             Found = TRUE;
  751.  
  752.             /* Cache normal. */
  753.  
  754.             if (!Test_Flag(HField, SMOOTHED_FLAG))
  755.             {
  756.               if (HField->Data->cache_pos < HF_CACHE_SIZE)
  757.               {
  758.                 Assign_Vector(HField->Data->Normal_Cache[HField->Data->cache_pos].normal, N);
  759.  
  760.                 HField->Data->Normal_Cache[HField->Data->cache_pos].fx = x + s;
  761.                 HField->Data->Normal_Cache[HField->Data->cache_pos].fz = z + t;
  762.  
  763.                 HField->Data->cache_pos++;
  764.               }
  765.             }
  766.           }
  767.         }
  768.       }
  769.     }
  770.   }
  771.  
  772.   return(Found);
  773. }
  774.  
  775.  
  776.  
  777. /*****************************************************************************
  778. *
  779. * FUNCTION
  780. *
  781. *   add_single_normal
  782. *
  783. * INPUT
  784. *   
  785. * OUTPUT
  786. *   
  787. * RETURNS
  788. *   
  789. * AUTHOR
  790. *
  791. *   Doug Muir, David Buck, Drew Wells
  792. *   
  793. * DESCRIPTION
  794. *
  795. *   -
  796. *
  797. * CHANGES
  798. *
  799. *   -
  800. *
  801. ******************************************************************************/
  802.  
  803. static int add_single_normal(HF_VAL **data, int xsize, int zsize, int x0, int z0, int x1, int z1, int x2, int z2, VECTOR N)
  804. {
  805.   VECTOR v0, v1, v2, t0, t1, Nt;
  806.  
  807.   if ((x0 < 0) || (z0 < 0) ||
  808.       (x1 < 0) || (z1 < 0) ||
  809.       (x2 < 0) || (z2 < 0) ||
  810.       (x0 > xsize) || (z0 > zsize) ||
  811.       (x1 > xsize) || (z1 > zsize) ||
  812.       (x2 > xsize) || (z2 > zsize))
  813.   {
  814.     return(0);
  815.   }
  816.   else
  817.   {
  818.     Make_Vector(v0, x0, (DBL)data[z0][x0], z0);
  819.     Make_Vector(v1, x1, (DBL)data[z1][x1], z1);
  820.     Make_Vector(v2, x2, (DBL)data[z2][x2], z2);
  821.  
  822.     VSub(t0, v2, v0);
  823.     VSub(t1, v1, v0);
  824.  
  825.     VCross(Nt, t0, t1);
  826.  
  827.     normalize(Nt, Nt);
  828.  
  829.     if (Nt[Y] < 0.0)
  830.     {
  831.       VScaleEq(Nt, -1.0);
  832.     }
  833.  
  834.     VAddEq(N, Nt);
  835.  
  836.     return(1);
  837.   }
  838. }
  839.  
  840.  
  841.  
  842. /*****************************************************************************
  843. *
  844. * FUNCTION
  845. *
  846. *   smooth_height_field
  847. *
  848. * INPUT
  849. *   
  850. * OUTPUT
  851. *   
  852. * RETURNS
  853. *   
  854. * AUTHOR
  855. *
  856. *   Doug Muir, David Buck, Drew Wells
  857. *   
  858. * DESCRIPTION
  859. *
  860. *   Given a height field that only contains an elevation grid, this
  861. *   routine will walk through the data and produce averaged normals
  862. *   for all points on the grid.
  863. *
  864. * CHANGES
  865. *
  866. *   -
  867. *
  868. ******************************************************************************/
  869.  
  870. static void smooth_height_field(HFIELD *HField, int xsize, int zsize)
  871. {
  872.   int i, j, k;
  873.   VECTOR N;
  874.   HF_VAL **map = HField->Data->Map;
  875.  
  876.   /* First off, allocate all the memory needed to store the normal information */
  877.  
  878.   HField->Data->Normals_Height = zsize+1;
  879.  
  880.   HField->Data->Normals = (HF_Normals **)POV_MALLOC((zsize+1)*sizeof(HF_Normals *), "height field normals");
  881.  
  882.   for (i = 0; i <= zsize; i++)
  883.   {
  884.     HField->Data->Normals[i] = (HF_Normals *)POV_MALLOC((xsize+1)*sizeof(HF_Normals), "height field normals");
  885.   }
  886.  
  887.   /*
  888.    * For now we will do it the hard way - by generating the normals
  889.    * individually for each elevation point.
  890.    */
  891.  
  892.   for (i = 0; i < zsize; i++)
  893.   {
  894.     COOPERATE_0
  895.  
  896.     for (j = 0; j < xsize; j++)
  897.     {
  898.       Make_Vector(N, 0.0, 0.0, 0.0);
  899.  
  900.       k = 0;
  901.  
  902.       k += add_single_normal(map, xsize, zsize, j, i, j+1, i, j, i+1, N);
  903.       k += add_single_normal(map, xsize, zsize, j, i, j, i+1, j-1, i, N);
  904.       k += add_single_normal(map, xsize, zsize, j, i, j-1, i, j, i-1, N);
  905.       k += add_single_normal(map, xsize, zsize, j, i, j, i-1, j+1, i, N);
  906.  
  907.       if (k == 0)
  908.       {
  909.         Error ("Failed to find any normals at: (%d, %d).\n", i, j);
  910.       }
  911.  
  912.       normalize(N, N);
  913.  
  914.       HField->Data->Normals[i][j][0] = (short)(32767 * N[X]);
  915.       HField->Data->Normals[i][j][1] = (short)(32767 * N[Y]);
  916.       HField->Data->Normals[i][j][2] = (short)(32767 * N[Z]);
  917.     }
  918.   }
  919. }
  920.  
  921.  
  922.  
  923. /*****************************************************************************
  924. *
  925. * FUNCTION
  926. *
  927. *   Compute_HField
  928. *
  929. * INPUT
  930. *
  931. * OUTPUT
  932. *
  933. * RETURNS
  934. *
  935. * AUTHOR
  936. *
  937. *   Doug Muir, David Buck, Drew Wells
  938. *
  939. * DESCRIPTION
  940. *
  941. *   Copy image data into height field map. Create bounding blocks
  942. *   for the block traversal. Calculate normals for smoothed height fields.
  943. *
  944. * CHANGES
  945. *
  946. *   Feb 1995 : Modified to work with new intersection functions. [DB]
  947. *
  948. ******************************************************************************/
  949.  
  950. void Compute_HField(HFIELD *HField, IMAGE *Image)
  951. {
  952.   int x, z, max_x, max_z;
  953.   int temp1 = 0, temp2 = 0;
  954.   HF_VAL min_y, max_y, temp_y;
  955.  
  956.   /* Get height field map size. */
  957.  
  958.   max_x = Image->iwidth;
  959.  
  960.   if (Image->File_Type == POT_FILE)
  961.   {
  962.     max_x = max_x / 2;
  963.   }
  964.  
  965.   max_z = Image->iheight;
  966.  
  967.   /* Allocate memory for map. */
  968.  
  969.   HField->Data->Map = (HF_VAL **)POV_MALLOC(max_z*sizeof(HF_VAL *), "height field");
  970.  
  971.   for (z = 0; z < max_z; z++)
  972.   {
  973.     HField->Data->Map[z] = (HF_VAL *)POV_MALLOC(max_x*sizeof(HF_VAL), "height field");
  974.   }
  975.  
  976.   /* Copy map. */
  977.  
  978.   min_y = 65535L;
  979.   max_y = 0;
  980.  
  981.   for (z = 0; z < max_z; z++)
  982.   {
  983.     COOPERATE_0
  984.     for (x = 0; x < max_x; x++)
  985.     {
  986.  
  987.       switch (Image->File_Type)
  988.       {
  989. #ifndef BurnAllGifs
  990.         case GIF_FILE:
  991.  
  992.           temp1 = Image->data.map_lines[max_z - z - 1][x];
  993.           temp2 = 0;
  994.  
  995.           break;
  996. #endif
  997.         case POT_FILE:
  998.  
  999.           temp1 = Image->data.map_lines[max_z - z - 1][x];
  1000.           temp2 = Image->data.map_lines[max_z - z - 1][x + max_x];
  1001.  
  1002.           break;
  1003.  
  1004.  
  1005.         case PPM_FILE:
  1006.  
  1007.           temp1 = Image->data.rgb_lines[max_z - z - 1].red[x];
  1008.           temp2 = Image->data.rgb_lines[max_z - z - 1].green[x];
  1009.  
  1010.           break;
  1011.  
  1012.         case PGM_FILE:
  1013.         case TGA_FILE:
  1014.         case PNG_FILE:
  1015.  
  1016.           if (Image->Colour_Map == NULL)
  1017.           {
  1018.             temp1 = Image->data.rgb_lines[max_z - z - 1].red[x];
  1019.             temp2 = Image->data.rgb_lines[max_z - z - 1].green[x];
  1020.           }
  1021.           else
  1022.           {
  1023.             temp1 = Image->data.map_lines[max_z - z - 1][x];
  1024.             temp2 = 0;
  1025.           }
  1026.  
  1027.           break;
  1028.  
  1029.         default:
  1030.  
  1031.           Error("Unknown image type in Compute_HField().\n");
  1032.       }
  1033.  
  1034.       temp_y = (HF_VAL)(256*temp1 + temp2);
  1035.  
  1036.       HField->Data->Map[z][x] = temp_y;
  1037.  
  1038.       min_y = min(min_y, temp_y);
  1039.       max_y = max(max_y, temp_y);
  1040.     }
  1041.   }
  1042.  
  1043.   /* Resize bounding box. */
  1044.  
  1045.   HField->Data->min_y = min_y;
  1046.   HField->Data->max_y = max_y;
  1047.  
  1048.   HField->bounding_box->bounds[0][Y] = max((DBL)min_y, HField->bounding_box->bounds[0][Y]) - HFIELD_OFFSET;
  1049.   HField->bounding_box->bounds[1][Y] = (DBL)max_y + HFIELD_OFFSET;
  1050.  
  1051.   /* Compute smoothed height field. */
  1052.  
  1053.   if (Test_Flag(HField, SMOOTHED_FLAG))
  1054.   {
  1055.     smooth_height_field(HField, max_x-1, max_z-1);
  1056.   }
  1057.  
  1058.   HField->Data->max_x = max_x-2;
  1059.   HField->Data->max_z = max_z-2;
  1060.  
  1061.   build_hfield_blocks(HField);
  1062. }
  1063.  
  1064.  
  1065.  
  1066. /*****************************************************************************
  1067. *
  1068. * FUNCTION
  1069. *
  1070. *   build_hfield_blocks
  1071. *
  1072. * INPUT
  1073. *   
  1074. * OUTPUT
  1075. *   
  1076. * RETURNS
  1077. *   
  1078. * AUTHOR
  1079. *
  1080. *   Dieter Bayer
  1081. *   
  1082. * DESCRIPTION
  1083. *
  1084. *   Create the bounding block hierarchy used by the block traversal.
  1085. *
  1086. * CHANGES
  1087. *
  1088. *   Feb 1995 : Creation.
  1089. *
  1090. ******************************************************************************/
  1091.  
  1092. static void build_hfield_blocks(HFIELD *HField)
  1093. {
  1094.   int x, z, nx, nz, wx, wz;
  1095.   int i, j;
  1096.   int xmin, xmax, zmin, zmax;
  1097.   DBL y, ymin, ymax, water;
  1098.  
  1099.   /* Get block size. */
  1100.  
  1101.   nx = max(1, (int)(sqrt((DBL)HField->Data->max_x)));
  1102.   nz = max(1, (int)(sqrt((DBL)HField->Data->max_z)));
  1103.  
  1104.   /* Get dimensions of sub-block. */
  1105.  
  1106.   wx = (int)ceil((DBL)(HField->Data->max_x + 2) / (DBL)nx);
  1107.   wz = (int)ceil((DBL)(HField->Data->max_z + 2) / (DBL)nz);
  1108.  
  1109.   /* Increase number of sub-blocks if necessary. */
  1110.  
  1111.   if (nx * wx < HField->Data->max_x + 2)
  1112.   {
  1113.     nx++;
  1114.   }
  1115.  
  1116.   if (nz * wz < HField->Data->max_z + 2)
  1117.   {
  1118.     nz++;
  1119.   }
  1120.  
  1121.   if (!Test_Flag(HField, HIERARCHY_FLAG) || ((nx == 1) && (nz == 1)))
  1122.   {
  1123.     /* We don't want a bounding hierarchy. Just use one block. */
  1124.  
  1125.     HField->Data->Block = (HFIELD_BLOCK **)POV_MALLOC(sizeof(HFIELD_BLOCK *), "height field blocks");
  1126.  
  1127.     HField->Data->Block[0] = (HFIELD_BLOCK *)POV_MALLOC(sizeof(HFIELD_BLOCK), "height field blocks");
  1128.  
  1129.     HField->Data->Block[0][0].xmin = 0;
  1130.     HField->Data->Block[0][0].xmax = HField->Data->max_x;
  1131.     HField->Data->Block[0][0].zmin = 0;
  1132.     HField->Data->Block[0][0].zmax = HField->Data->max_z;
  1133.  
  1134.     HField->Data->Block[0][0].ymin = HField->bounding_box->bounds[0][Y];
  1135.     HField->Data->Block[0][0].ymax = HField->bounding_box->bounds[1][Y];
  1136.  
  1137.     HField->Data->block_max_x = 1;
  1138.     HField->Data->block_max_z = 1;
  1139.  
  1140.     HField->Data->block_width_x = HField->Data->max_x + 2;
  1141.     HField->Data->block_width_z = HField->Data->max_y + 2;
  1142.  
  1143. /*
  1144.     Debug_Info("\nHeight field: %d x %d (1 x 1 blocks)", HField->Data->max_x+2, HField->Data->max_z+2);
  1145. */
  1146.  
  1147.     return;
  1148.   }
  1149.  
  1150. /*
  1151.   Debug_Info("\nHeight field: %d x %d (%d x %d blocks)", HField->Data->max_x+2, HField->Data->max_z+2, nx, nz);
  1152. */
  1153.  
  1154.   /* Allocate memory for blocks. */
  1155.  
  1156.   HField->Data->Block = (HFIELD_BLOCK **)POV_MALLOC(nz*sizeof(HFIELD_BLOCK *), "height field blocks");
  1157.  
  1158.   /* Store block information. */
  1159.  
  1160.   HField->Data->block_max_x = nx;
  1161.   HField->Data->block_max_z = nz;
  1162.  
  1163.   HField->Data->block_width_x = wx;
  1164.   HField->Data->block_width_z = wz;
  1165.  
  1166.   water = HField->bounding_box->bounds[0][Y];
  1167.  
  1168.   for (z = 0; z < nz; z++)
  1169.   {
  1170.     COOPERATE_1
  1171.  
  1172.     HField->Data->Block[z] = (HFIELD_BLOCK *)POV_MALLOC(nx*sizeof(HFIELD_BLOCK), "height field blocks");
  1173.  
  1174.     for (x = 0; x < nx; x++)
  1175.     {
  1176.       /* Get block's borders. */
  1177.  
  1178.       xmin = x * wx;
  1179.       zmin = z * wz;
  1180.  
  1181.       xmax = min((x + 1) * wx - 1, HField->Data->max_x);
  1182.       zmax = min((z + 1) * wz - 1, HField->Data->max_z);
  1183.  
  1184.       /* Find min. and max. height in current block. */
  1185.  
  1186.       ymin = BOUND_HUGE;
  1187.       ymax = -BOUND_HUGE;
  1188.  
  1189.       for (i = xmin; i <= xmax+1; i++)
  1190.       {
  1191.         for (j = zmin; j <= zmax+1; j++)
  1192.         {
  1193.           y = Get_Height(i, j, HField);
  1194.  
  1195.           ymin = min(ymin, y);
  1196.           ymax = max(ymax, y);
  1197.         }
  1198.       }
  1199.  
  1200.       /* Store block's borders. */
  1201.  
  1202.       HField->Data->Block[z][x].xmin = xmin;
  1203.       HField->Data->Block[z][x].xmax = xmax;
  1204.       HField->Data->Block[z][x].zmin = zmin;
  1205.       HField->Data->Block[z][x].zmax = zmax;
  1206.  
  1207.       HField->Data->Block[z][x].ymin = max(ymin, water) - HFIELD_OFFSET;
  1208.       HField->Data->Block[z][x].ymax = ymax + HFIELD_OFFSET;
  1209.     }
  1210.   }
  1211. }
  1212.  
  1213.  
  1214.  
  1215. /*****************************************************************************
  1216. *
  1217. * FUNCTION
  1218. *
  1219. *   Translate_HField
  1220. *
  1221. * INPUT
  1222. *   
  1223. * OUTPUT
  1224. *   
  1225. * RETURNS
  1226. *   
  1227. * AUTHOR
  1228. *
  1229. *   Doug Muir, David Buck, Drew Wells
  1230. *
  1231. * DESCRIPTION
  1232. *
  1233. *   -
  1234. *
  1235. * CHANGES
  1236. *
  1237. *   -
  1238. *
  1239. ******************************************************************************/
  1240.  
  1241. static void Translate_HField (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
  1242. {
  1243.   Transform_HField(Object, Trans);
  1244. }
  1245.  
  1246.  
  1247.  
  1248. /*****************************************************************************
  1249. *
  1250. * FUNCTION
  1251. *
  1252. *   Rotate_HField
  1253. *
  1254. * INPUT
  1255. *
  1256. * OUTPUT
  1257. *
  1258. * RETURNS
  1259. *
  1260. * AUTHOR
  1261. *
  1262. *   Doug Muir, David Buck, Drew Wells
  1263. *
  1264. * DESCRIPTION
  1265. *
  1266. *   -
  1267. *
  1268. * CHANGES
  1269. *
  1270. *   -
  1271. *
  1272. ******************************************************************************/
  1273.  
  1274. static void Rotate_HField (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
  1275. {
  1276.   Transform_HField(Object, Trans);
  1277. }
  1278.  
  1279.  
  1280.  
  1281. /*****************************************************************************
  1282. *
  1283. * FUNCTION
  1284. *
  1285. *   Scale_HField
  1286. *
  1287. * INPUT
  1288. *   
  1289. * OUTPUT
  1290. *   
  1291. * RETURNS
  1292. *   
  1293. * AUTHOR
  1294. *
  1295. *   Doug Muir, David Buck, Drew Wells
  1296. *   
  1297. * DESCRIPTION
  1298. *
  1299. *   -
  1300. *
  1301. * CHANGES
  1302. *
  1303. *   -
  1304. *
  1305. ******************************************************************************/
  1306.  
  1307. static void Scale_HField (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
  1308. {
  1309.   Transform_HField(Object, Trans);
  1310. }
  1311.  
  1312.  
  1313.  
  1314. /*****************************************************************************
  1315. *
  1316. * FUNCTION
  1317. *
  1318. *   Invert_HField
  1319. *
  1320. * INPUT
  1321. *   
  1322. * OUTPUT
  1323. *   
  1324. * RETURNS
  1325. *   
  1326. * AUTHOR
  1327. *
  1328. *   Doug Muir, David Buck, Drew Wells
  1329. *   
  1330. * DESCRIPTION
  1331. *
  1332. *   -
  1333. *
  1334. * CHANGES
  1335. *
  1336. *   -
  1337. *
  1338. ******************************************************************************/
  1339.  
  1340. static void Invert_HField (OBJECT *Object)
  1341. {
  1342.   Invert_Flag(Object, INVERTED_FLAG);
  1343. }
  1344.  
  1345.  
  1346.  
  1347. /*****************************************************************************
  1348. *
  1349. * FUNCTION
  1350. *
  1351. *   Transform_HField
  1352. *
  1353. * INPUT
  1354. *   
  1355. * OUTPUT
  1356. *   
  1357. * RETURNS
  1358. *   
  1359. * AUTHOR
  1360. *
  1361. *   Doug Muir, David Buck, Drew Wells
  1362. *   
  1363. * DESCRIPTION
  1364. *
  1365. *   -
  1366. *
  1367. * CHANGES
  1368. *
  1369. *   -
  1370. *
  1371. ******************************************************************************/
  1372.  
  1373. static void Transform_HField (OBJECT *Object, TRANSFORM *Trans)
  1374. {
  1375.   Compose_Transforms(((HFIELD *)Object)->Trans, Trans);
  1376.  
  1377.   Compute_HField_BBox((HFIELD *)Object);
  1378. }
  1379.  
  1380.  
  1381.  
  1382. /*****************************************************************************
  1383. *
  1384. * FUNCTION
  1385. *
  1386. *   Create_HField
  1387. *
  1388. * INPUT
  1389. *
  1390. * OUTPUT
  1391. *   
  1392. * RETURNS
  1393. *   
  1394. * AUTHOR
  1395. *
  1396. *   Doug Muir, David Buck, Drew Wells
  1397. *   
  1398. * DESCRIPTION
  1399. *
  1400. *   Allocate and intialize a height field.
  1401. *
  1402. * CHANGES
  1403. *
  1404. *   Feb 1995 : Modified to work with new intersection functions. [DB]
  1405. *
  1406. *   June 1999 : Modified by Lummox JR to set smooth_type to an initial value
  1407. *               (0 default)
  1408. *
  1409. ******************************************************************************/
  1410.  
  1411. HFIELD *Create_HField()
  1412. {
  1413.   HFIELD *New;
  1414.  
  1415.   /* Allocate height field. */
  1416.  
  1417.   New = (HFIELD *)POV_MALLOC(sizeof(HFIELD), "height field");
  1418.  
  1419.   INIT_OBJECT_FIELDS(New, HFIELD_OBJECT, &HField_Methods)
  1420.  
  1421.   /* Always uses Trans so always create one. */
  1422.  
  1423.   New->Trans = Create_Transform();
  1424.  
  1425.   New->bounding_box = Create_Box();
  1426.  
  1427.   /* Allocate height field data. */
  1428.  
  1429.   New->Data = (HFIELD_DATA *)POV_MALLOC(sizeof(HFIELD_DATA), "height field");
  1430.  
  1431.   New->Data->References = 1;
  1432.  
  1433.   New->Data->cache_pos = 0;
  1434.  
  1435.   New->Data->Normals_Height = 0;
  1436.  
  1437.   New->Data->Map     = NULL;
  1438.   New->Data->Normals = NULL;
  1439.  
  1440.   New->Data->max_x = 0;
  1441.   New->Data->max_z = 0;
  1442.  
  1443.   New->Data->block_max_x = 0;
  1444.   New->Data->block_max_z = 0;
  1445.  
  1446.   New->Data->block_width_x = 0;
  1447.   New->Data->block_width_z = 0;
  1448.  
  1449.   Set_Flag(New, HIERARCHY_FLAG);
  1450.  
  1451.   New->smooth_type = 0;
  1452.  
  1453.   return(New);
  1454. }
  1455.  
  1456.  
  1457.  
  1458. /*****************************************************************************
  1459. *
  1460. * FUNCTION
  1461. *
  1462. *   Copy_HField
  1463. *
  1464. * INPUT
  1465. *   
  1466. * OUTPUT
  1467. *   
  1468. * RETURNS
  1469. *   
  1470. * AUTHOR
  1471. *
  1472. *   Doug Muir, David Buck, Drew Wells
  1473. *   
  1474. * DESCRIPTION
  1475. *
  1476. *   NOTE: The height field data is not copied, only the number of references
  1477. *         is counted, so that Destray_HField() knows if it can be destroyed.
  1478. *
  1479. * CHANGES
  1480. *
  1481. ******************************************************************************/
  1482.  
  1483. static HFIELD *Copy_HField(OBJECT *Object)
  1484. {
  1485.   HFIELD *New;
  1486.  
  1487.   New = Create_HField();
  1488.  
  1489.   /* Destroy obsolete things created in Create_HField(). */
  1490.  
  1491.   Destroy_Transform(New->Trans);
  1492.  
  1493.   Destroy_Box((OBJECT *)(New->bounding_box));
  1494.  
  1495.   POV_FREE (New->Data);
  1496.  
  1497.   /* Copy height field. */
  1498.  
  1499.   *New = *((HFIELD *)Object);
  1500.  
  1501.   New->Trans = Copy_Transform(((HFIELD *)Object)->Trans);
  1502.  
  1503.   New->bounding_box = Copy_Box((OBJECT *)(((HFIELD *)Object)->bounding_box));
  1504.  
  1505.   New->Data->References++;
  1506.  
  1507.   return(New);
  1508. }
  1509.  
  1510.  
  1511.  
  1512. /*****************************************************************************
  1513. *
  1514. * FUNCTION
  1515. *
  1516. *   Destroy_HField
  1517. *
  1518. * INPUT
  1519. *   
  1520. * OUTPUT
  1521. *   
  1522. * RETURNS
  1523. *   
  1524. * AUTHOR
  1525. *
  1526. *   Doug Muir, David Buck, Drew Wells
  1527. *   
  1528. * DESCRIPTION
  1529. *
  1530. *   NOTE: The height field data is destroyed if it's no longer
  1531. *         used by any copy.
  1532. *
  1533. * CHANGES
  1534. *
  1535. *   Feb 1995 : Modified to work with new intersection functions. [DB]
  1536. *
  1537. ******************************************************************************/
  1538.  
  1539. static void Destroy_HField (OBJECT *Object)
  1540. {
  1541.   int i;
  1542.   HFIELD *HField = (HFIELD *)Object;
  1543.  
  1544.   Destroy_Transform(HField->Trans);
  1545.  
  1546.   Destroy_Box((OBJECT *)(HField->bounding_box));
  1547.  
  1548.   if (--(HField->Data->References) == 0)
  1549.   {
  1550.     if (HField->Data->Map != NULL)
  1551.     {
  1552.       for (i = 0; i < HField->Data->max_z+2; i++)
  1553.       {
  1554.         if (HField->Data->Map[i] != NULL)
  1555.         {
  1556.           POV_FREE (HField->Data->Map[i]);
  1557.         }
  1558.       }
  1559.  
  1560.       POV_FREE (HField->Data->Map);
  1561.     }
  1562.  
  1563.     if (HField->Data->Normals != NULL)
  1564.     {
  1565.       for (i = 0; i < HField->Data->Normals_Height; i++)
  1566.       {
  1567.         POV_FREE (HField->Data->Normals[i]);
  1568.       }
  1569.  
  1570.       POV_FREE (HField->Data->Normals);
  1571.     }
  1572.  
  1573.     if (HField->Data->Block != NULL)
  1574.     {
  1575.       for (i = 0; i < HField->Data->block_max_z; i++)
  1576.       {
  1577.         POV_FREE(HField->Data->Block[i]);
  1578.       }
  1579.  
  1580.       POV_FREE(HField->Data->Block);
  1581.     }
  1582.  
  1583.     POV_FREE (HField->Data);
  1584.   }
  1585.  
  1586.   POV_FREE (Object);
  1587. }
  1588.  
  1589.  
  1590.  
  1591. /*****************************************************************************
  1592. *
  1593. * FUNCTION
  1594. *
  1595. *   Compute_HField_BBox
  1596. *
  1597. * INPUT
  1598. *
  1599. *   HField - Height field
  1600. *   
  1601. * OUTPUT
  1602. *
  1603. *   HField
  1604. *   
  1605. * RETURNS
  1606. *   
  1607. * AUTHOR
  1608. *
  1609. *   Dieter Bayer
  1610. *   
  1611. * DESCRIPTION
  1612. *
  1613. *   Calculate the bounding box of a height field.
  1614. *
  1615. * CHANGES
  1616. *
  1617. *   Aug 1994 : Creation.
  1618. *
  1619. ******************************************************************************/
  1620.  
  1621. void Compute_HField_BBox(HFIELD *HField)
  1622. {
  1623.   Assign_BBox_Vect(HField->BBox.Lower_Left, HField->bounding_box->bounds[0]);
  1624.  
  1625.   VSub (HField->BBox.Lengths, HField->bounding_box->bounds[1], HField->bounding_box->bounds[0]);
  1626.  
  1627.   if (HField->Trans != NULL)
  1628.   {
  1629.     Recompute_BBox(&HField->BBox, HField->Trans);
  1630.   }
  1631. }
  1632.  
  1633.  
  1634.  
  1635. /*****************************************************************************
  1636. *
  1637. * FUNCTION
  1638. *
  1639. *   dda_traversal
  1640. *
  1641. * INPUT
  1642. *
  1643. *   Ray    - Current ray
  1644. *   HField - Height field
  1645. *   Start  - Start point for the walk
  1646. *   Block  - Sub-block of the height field to traverse
  1647. *   
  1648. * OUTPUT
  1649. *   
  1650. * RETURNS
  1651. *
  1652. *   int - TRUE if intersection was found
  1653. *   
  1654. * AUTHOR
  1655. *
  1656. *   Dieter Bayer
  1657. *   
  1658. * DESCRIPTION
  1659. *
  1660. *   Traverse the grid cell of the height field using a modified DDA.
  1661. *
  1662. *   Based on the following article:
  1663. *
  1664. *     Musgrave, F. Kenton, "Grid Tracing: Fast Ray Tracing for Height
  1665. *     Fields", Research Report YALEU-DCS-RR-39, Yale University, July 1988
  1666. *
  1667. *   You should note that there are (n-1) x (m-1) grid cells in a height
  1668. *   field of (image) size n x m. A grid cell (i,j), 0 <= i <= n-1,
  1669. *   0 <= j <= m-1, extends from x = i to x = i + 1 - epsilon and
  1670. *   y = j to y = j + 1 -epsilon.
  1671. *
  1672. * CHANGES
  1673. *
  1674. *   Feb 1995 : Creation.
  1675. *
  1676. ******************************************************************************/
  1677.  
  1678. static int dda_traversal(RAY *Ray, HFIELD *HField, VECTOR Start, HFIELD_BLOCK *Block)
  1679. {
  1680.   int found;
  1681.   int xmin, xmax, zmin, zmax;
  1682.   int x, z, signx, signz;
  1683.   DBL ymin, ymax, y1, y2;
  1684.   DBL px, pz, dx, dy, dz;
  1685.   DBL delta, error, x0, z0;
  1686.   DBL neary, fary, deltay;
  1687.  
  1688.   /* Setup DDA. */
  1689.  
  1690.   found = FALSE;
  1691.  
  1692.   px = Start[X];
  1693.   pz = Start[Z];
  1694.  
  1695.   /* Get dimensions of current block. */
  1696.  
  1697.   xmin = Block->xmin;
  1698.   xmax = min(Block->xmax + 1, HField->Data->max_x);
  1699.   zmin = Block->zmin;
  1700.   zmax = min(Block->zmax + 1, HField->Data->max_z);
  1701.  
  1702.   ymin = min(Start[Y], Block->ymin) - EPSILON;
  1703.   ymax = max(Start[Y], Block->ymax) + EPSILON;
  1704.  
  1705.   /* Check for illegal grid values (caused by numerical inaccuracies). */
  1706.  
  1707.   if (px < (DBL)xmin)
  1708.   {
  1709.     if (px < (DBL)xmin - HFIELD_OFFSET)
  1710.     {
  1711.       Debug_Info("Illegal grid value in dda_traversal().\n");
  1712.  
  1713.       return(FALSE);
  1714.     }
  1715.     else
  1716.     {
  1717.       px = (DBL)xmin;
  1718.     }
  1719.   }
  1720.   else
  1721.   {
  1722.     if (px > (DBL)xmax + 1.0 - EPSILON)
  1723.     {
  1724.       if (px > (DBL)xmax + 1.0 + EPSILON)
  1725.       {
  1726.         Debug_Info("Illegal grid value in dda_traversal().\n");
  1727.  
  1728.         return(FALSE);
  1729.       }
  1730.       else
  1731.       {
  1732.         px = (DBL)xmax + 1.0 - EPSILON;
  1733.       }
  1734.     }
  1735.   }
  1736.  
  1737.   if (pz < (DBL)zmin)
  1738.   {
  1739.     if (pz < (DBL)zmin - HFIELD_OFFSET)
  1740.     {
  1741.       Debug_Info("Illegal grid value in dda_traversal().\n");
  1742.  
  1743.       return(FALSE);
  1744.     }
  1745.     else
  1746.     {
  1747.       pz = (DBL)zmin;
  1748.     }
  1749.   }
  1750.   else
  1751.   {
  1752.     if (pz > (DBL)zmax + 1.0 - EPSILON)
  1753.     {
  1754.       if (pz > (DBL)zmax + 1.0 + EPSILON)
  1755.       {
  1756.         Debug_Info("Illegal grid value in dda_traversal().\n");
  1757.  
  1758.         return(FALSE);
  1759.       }
  1760.       else
  1761.       {
  1762.         pz = (DBL)zmax + 1.0 - EPSILON;
  1763.       }
  1764.     }
  1765.   }
  1766.  
  1767.   dx = Ray->Direction[X];
  1768.   dy = Ray->Direction[Y];
  1769.   dz = Ray->Direction[Z];
  1770.  
  1771.   /*
  1772.    * Here comes the DDA algorithm.
  1773.    */
  1774.  
  1775.   /* Choose algorithm depending on the driving axis. */
  1776.  
  1777.   if (fabs(dx) >= fabs(dz))
  1778.   {
  1779.     /*
  1780.      * X-axis is driving axis.
  1781.      */
  1782.  
  1783.     delta = fabs(dz / dx);
  1784.  
  1785.     x = (int)px;
  1786.     z = (int)pz;
  1787.  
  1788.     x0 = px - floor(px);
  1789.     z0 = pz - floor(pz);
  1790.  
  1791.     signx = sign(dx);
  1792.     signz = sign(dz);
  1793.  
  1794.     /* Get initial error. */
  1795.  
  1796.     if (dx >= 0.0)
  1797.     {
  1798.       if (dz >= 0.0)
  1799.       {
  1800.         error = z0 + delta * (1.0 - x0) - 1.0;
  1801.       }
  1802.       else
  1803.       {
  1804.         error = -(z0 - delta * (1.0 - x0));
  1805.       }
  1806.     }
  1807.     else
  1808.     {
  1809.       if (dz >= 0.0)
  1810.       {
  1811.         error = z0 + delta * x0 - 1.0;
  1812.       }
  1813.       else
  1814.       {
  1815.         error = -(z0 - delta * x0);
  1816.       }
  1817.     }
  1818.  
  1819.     /* Get y differential. */
  1820.  
  1821.     deltay = dy / fabs(dx);
  1822.  
  1823.     if (dx >= 0.0)
  1824.     {
  1825.       neary = Start[Y] - x0 * deltay;
  1826.  
  1827.       fary = neary + deltay;
  1828.     }
  1829.     else
  1830.     {
  1831.       neary = Start[Y] - (1.0 - x0) * deltay;
  1832.  
  1833.       fary = neary + deltay;
  1834.     }
  1835.  
  1836.     /* Step through the cells. */
  1837.  
  1838.     do
  1839.     {
  1840.       if (neary < fary)
  1841.       {
  1842.         y1 = neary;
  1843.         y2 = fary;
  1844.       }
  1845.       else
  1846.       {
  1847.         y1 = fary;
  1848.         y2 = neary;
  1849.       }
  1850.  
  1851.       if (intersect_pixel(x, z, Ray, HField, y1, y2))
  1852.       {
  1853.         if (HField->Type & IS_CHILD_OBJECT)
  1854.         {
  1855.           found = TRUE;
  1856.         }
  1857.         else
  1858.         {
  1859.           return(TRUE);
  1860.         }
  1861.       }
  1862.  
  1863.       if (error > EPSILON)
  1864.       {
  1865.         z += signz;
  1866.  
  1867.         if ((z < zmin) || (z > zmax))
  1868.         {
  1869.           break;
  1870.         }
  1871.         else
  1872.         {
  1873.           if (intersect_pixel(x, z, Ray, HField, y1, y2))
  1874.           {
  1875.             if (HField->Type & IS_CHILD_OBJECT)
  1876.             {
  1877.               found = TRUE;
  1878.             }
  1879.             else
  1880.             {
  1881.               return(TRUE);
  1882.             }
  1883.           }
  1884.         }
  1885.  
  1886.         error--;
  1887.       }
  1888.       else
  1889.       {
  1890.         if (error > -EPSILON)
  1891.         {
  1892.           z += signz;
  1893.  
  1894.           error--;
  1895.         }
  1896.       }
  1897.  
  1898.       x += signx;
  1899.  
  1900.       error += delta;
  1901.  
  1902.       neary = fary;
  1903.  
  1904.       fary += deltay;
  1905.     }
  1906.     while ((neary >= ymin) && (neary <= ymax) && (x >= xmin) && (x <= xmax) && (z >= zmin) && (z <= zmax));
  1907.   }
  1908.   else
  1909.   {
  1910.     /*
  1911.      * Z-axis is driving axis.
  1912.      */
  1913.  
  1914.     delta = fabs(dx / dz);
  1915.  
  1916.     x = (int)px;
  1917.     z = (int)pz;
  1918.  
  1919.     x0 = px - floor(px);
  1920.     z0 = pz - floor(pz);
  1921.  
  1922.     signx = sign(dx);
  1923.     signz = sign(dz);
  1924.  
  1925.     /* Get initial error. */
  1926.  
  1927.     if (dz >= 0.0)
  1928.     {
  1929.       if (dx >= 0.0)
  1930.       {
  1931.         error = x0 + delta * (1.0 - z0) - 1.0;
  1932.       }
  1933.       else
  1934.       {
  1935.         error = -(x0 - delta * (1.0 - z0));
  1936.       }
  1937.     }
  1938.     else
  1939.     {
  1940.       if (dx >= 0.0)
  1941.       {
  1942.         error = x0 + delta * z0 - 1.0;
  1943.       }
  1944.       else
  1945.       {
  1946.         error = -(x0 - delta * z0);
  1947.       }
  1948.     }
  1949.  
  1950.     /* Get y differential. */
  1951.  
  1952.     deltay = dy / fabs(dz);
  1953.  
  1954.     if (dz >= 0.0)
  1955.     {
  1956.       neary = Start[Y] - z0 * deltay;
  1957.  
  1958.       fary = neary + deltay;
  1959.     }
  1960.     else
  1961.     {
  1962.       neary = Start[Y] - (1.0 - z0) * deltay;
  1963.  
  1964.       fary = neary + deltay;
  1965.     }
  1966.  
  1967.     /* Step through the cells. */
  1968.  
  1969.     do
  1970.     {
  1971.       if (neary < fary)
  1972.       {
  1973.         y1 = neary;
  1974.         y2 = fary;
  1975.       }
  1976.       else
  1977.       {
  1978.         y1 = fary;
  1979.         y2 = neary;
  1980.       }
  1981.  
  1982.       if (intersect_pixel(x, z, Ray, HField, y1, y2))
  1983.       {
  1984.         if (HField->Type & IS_CHILD_OBJECT)
  1985.         {
  1986.           found = TRUE;
  1987.         }
  1988.         else
  1989.         {
  1990.           return(TRUE);
  1991.         }
  1992.       }
  1993.  
  1994.       if (error > EPSILON)
  1995.       {
  1996.         x += signx;
  1997.  
  1998.         if ((x < xmin) || (x > xmax))
  1999.         {
  2000.           break;
  2001.         }
  2002.         else
  2003.         {
  2004.           if (intersect_pixel(x, z, Ray, HField, y1, y2))
  2005.           {
  2006.             if (HField->Type & IS_CHILD_OBJECT)
  2007.             {
  2008.               found = TRUE;
  2009.             }
  2010.             else
  2011.             {
  2012.               return(TRUE);
  2013.             }
  2014.           }
  2015.         }
  2016.  
  2017.         error--;
  2018.       }
  2019.       else
  2020.       {
  2021.         if (error > -EPSILON)
  2022.         {
  2023.           x += signx;
  2024.  
  2025.           error--;
  2026.         }
  2027.       }
  2028.  
  2029.       z += signz;
  2030.  
  2031.       error += delta;
  2032.  
  2033.       neary = fary;
  2034.  
  2035.       fary += deltay;
  2036.     }
  2037.     while ((neary >= ymin-EPSILON) && (neary <= ymax+EPSILON) &&
  2038.            (x >= xmin) && (x <= xmax) &&
  2039.            (z >= zmin) && (z <= zmax));
  2040.   }
  2041.  
  2042.   return(found);
  2043. }
  2044.  
  2045.  
  2046.  
  2047. /*****************************************************************************
  2048. *
  2049. * FUNCTION
  2050. *
  2051. *   block_traversal
  2052. *
  2053. * INPUT
  2054. *
  2055. *   Ray    - Current ray
  2056. *   HField - Height field
  2057. *   Start  - Start point for the walk
  2058. *
  2059. * OUTPUT
  2060. *
  2061. * RETURNS
  2062. *
  2063. *   int - TRUE if intersection was found
  2064. *   
  2065. * AUTHOR
  2066. *
  2067. *   Dieter Bayer
  2068. *
  2069. * DESCRIPTION
  2070. *
  2071. *   Traverse the blocks of the height field.
  2072. *
  2073. * CHANGES
  2074. *
  2075. *   Feb 1995 : Creation.
  2076. *
  2077. *   Aug 1996 : Fixed bug as reported by Dean M. Phillips:
  2078. *              "I found a bug in the height field code which resulted
  2079. *              in "Illegal grid value in dda_traversal()." messages
  2080. *              along with dark vertical lines in the height fields.
  2081. *              This turned out to be caused by overlooking the units in
  2082. *              which some boundary tests in two different places were
  2083. *              made. It was easy to fix.
  2084. *
  2085. ******************************************************************************/
  2086.  
  2087. static int block_traversal(RAY *Ray, HFIELD *HField, VECTOR Start)
  2088. {
  2089.   int xmax, zmax;
  2090.   int x, z, nx, nz, signx, signz;
  2091.   int found = FALSE;
  2092.   int dx_zero, dz_zero;
  2093.   DBL px, pz, dx, dy, dz;
  2094.   DBL maxdv;
  2095.   DBL ymin, ymax, y1, y2;
  2096.   DBL neary, fary;
  2097.   DBL k1, k2, dist;
  2098.   VECTOR nearP, farP;
  2099.   HFIELD_BLOCK *Block;
  2100.  
  2101.   px = Start[X];
  2102.   pz = Start[Z];
  2103.  
  2104.   dx = Ray->Direction[X];
  2105.   dy = Ray->Direction[Y];
  2106.   dz = Ray->Direction[Z];
  2107.  
  2108.   maxdv = (dx > dz) ? dx : dz;
  2109.  
  2110.   /* First test for 'perpendicular' rays. */
  2111.  
  2112.   if ((fabs(dx) < EPSILON) && (fabs(dz) < EPSILON))
  2113.   {
  2114.     x = (int)px;
  2115.     z = (int)pz;
  2116.  
  2117.     neary = Start[Y];
  2118.  
  2119.     if (dy >= 0.0)
  2120.     {
  2121.       fary = 65536.0;
  2122.     }
  2123.     else
  2124.     {
  2125.       fary = 0.0;
  2126.     }
  2127.  
  2128.     return(intersect_pixel(x, z, Ray, HField, min(neary, fary), max(neary, fary)));
  2129.   }
  2130.  
  2131.   /* If we don't have blocks we just step through the grid. */
  2132.  
  2133.   if ((HField->Data->block_max_x <= 1) && (HField->Data->block_max_z <= 1))
  2134.   {
  2135.     return(dda_traversal(Ray, HField, Start, &HField->Data->Block[0][0]));
  2136.   }
  2137.  
  2138.   /* Get dimensions of grid. */
  2139.  
  2140.   xmax = HField->Data->block_max_x;
  2141.   zmax = HField->Data->block_max_z;
  2142.  
  2143.   ymin = (DBL)HField->Data->min_y - EPSILON;
  2144.   ymax = (DBL)HField->Data->max_y + EPSILON;
  2145.  
  2146.   dx_zero = (fabs(dx) < EPSILON);
  2147.   dz_zero = (fabs(dz) < EPSILON);
  2148.  
  2149.   signx = sign(dx);
  2150.   signz = sign(dz);
  2151.  
  2152.   /* Walk on the block grid. */
  2153.  
  2154.   px /= HField->Data->block_width_x;
  2155.   pz /= HField->Data->block_width_z;
  2156.  
  2157.   x = (int)px;
  2158.   z = (int)pz;
  2159.  
  2160.   Assign_Vector(nearP, Start);
  2161.  
  2162.   neary = Start[Y];
  2163.  
  2164.   /*
  2165.    * Here comes the block walk algorithm.
  2166.    */
  2167.  
  2168.   do
  2169.   {
  2170. #ifdef HFIELD_EXTRA_STATS
  2171.     Increase_Counter(stats[Ray_HField_Block_Tests]);
  2172. #endif
  2173.  
  2174.     /* Get current block. */
  2175.  
  2176.     Block = &HField->Data->Block[z][x];
  2177.  
  2178.     /* Intersect ray with bounding planes. */
  2179.  
  2180.     if (dx_zero)
  2181.     {
  2182.       k1 = BOUND_HUGE;
  2183.     }
  2184.     else
  2185.     {
  2186.       if (signx >= 0)
  2187.       {
  2188.         k1 = ((DBL)Block->xmax + 1.0 - Ray->Initial[X]) / dx;
  2189.       }
  2190.       else
  2191.       {
  2192.         k1 = ((DBL)Block->xmin - Ray->Initial[X]) / dx;
  2193.       }
  2194.     }
  2195.  
  2196.     if (dz_zero)
  2197.     {
  2198.       k2 = BOUND_HUGE;
  2199.     }
  2200.     else
  2201.     {
  2202.       if (signz >= 0)
  2203.       {
  2204.         k2 = ((DBL)Block->zmax + 1.0 - Ray->Initial[Z]) / dz;
  2205.       }
  2206.       else
  2207.       {
  2208.         k2 = ((DBL)Block->zmin - Ray->Initial[Z]) / dz;
  2209.       }
  2210.     }
  2211.  
  2212.     /* Figure out the indices of the next block. */
  2213.  
  2214.     if (dz_zero || ((!dx_zero) && (k1<k2 - EPSILON / maxdv) && (k1>0.0)))
  2215. /*  if ((k1 < k2 - EPSILON / maxdv) && (k1 > 0.0)) */
  2216.     {
  2217.       /* Step along the x-axis. */
  2218.  
  2219.       dist = k1;
  2220.  
  2221.       nx = x + signx;
  2222.       nz = z;
  2223.     }
  2224.     else
  2225.     {
  2226.       if (dz_zero || ((!dx_zero) && (k1<k2 + EPSILON / maxdv) && (k1>0.0)))
  2227. /*    if ((k1 < k2 + EPSILON / maxdv) && (k1 > 0.0))  */
  2228.       {
  2229.         /* Step along both axis (very rare case). */
  2230.  
  2231.         dist = k1;
  2232.  
  2233.         nx = x + signx;
  2234.         nz = z + signz;
  2235.       }
  2236.       else
  2237.       {
  2238.         /* Step along the z-axis. */
  2239.  
  2240.         dist = k2;
  2241.  
  2242.         nx = x;
  2243.         nz = z + signz;
  2244.       }
  2245.     }
  2246.  
  2247.     /* Get point where ray leaves current block. */
  2248.  
  2249.     VEvaluateRay(farP, Ray->Initial, dist, Ray->Direction);
  2250.  
  2251.     fary = farP[Y];
  2252.  
  2253.     if (neary < fary)
  2254.     {
  2255.       y1 = neary;
  2256.       y2 = fary;
  2257.     }
  2258.     else
  2259.     {
  2260.       y1 = fary;
  2261.       y2 = neary;
  2262.     }
  2263.  
  2264.     /* Can we hit current block at all? */
  2265.  
  2266.     if ((y1 <= (DBL)Block->ymax + EPSILON) && (y2 >= (DBL)Block->ymin - EPSILON))
  2267.     {
  2268.       /* Test current block. */
  2269.  
  2270. #ifdef HFIELD_EXTRA_STATS
  2271.       Increase_Counter(stats[Ray_HField_Block_Tests_Succeeded]);
  2272. #endif
  2273.  
  2274.       if (dda_traversal(Ray, HField, nearP, &HField->Data->Block[z][x]))
  2275.       {
  2276.         if (HField->Type & IS_CHILD_OBJECT)
  2277.         {
  2278.           found = TRUE;
  2279.         }
  2280.         else
  2281.         {
  2282.           return(TRUE);
  2283.         }
  2284.       }
  2285.     }
  2286.  
  2287.     /* Step to next block. */
  2288.  
  2289.     x = nx;
  2290.     z = nz;
  2291.  
  2292.     Assign_Vector(nearP, farP);
  2293.  
  2294.     neary = fary;
  2295.   }
  2296.   while ((x >= 0) && (x < xmax) && (z >= 0) && (z < zmax) && (neary >= ymin) && (neary <= ymax));
  2297.  
  2298.   return(found);
  2299. }
  2300.  
  2301.