/****************************************************************************
*                   hfield.c
*
*    This file implements the height field shape primitive.  The shape is
*    implemented as a collection of triangles which are calculated as
*    needed.  The basic intersection routine first computes the rays
*    intersection with the box marking the limits of the shape, then
*    follows the line from one intersection point to the other, testing the
*    two triangles which form the pixel for an intersection with the ray at
*    each step.
*
*    height field added by Doug Muir with lots of advice and support
*    from David Buck and Drew Wells.
*
*  from Persistence of Vision(tm) Ray Tracer
*  Copyright 1996 Persistence of Vision Team
*---------------------------------------------------------------------------
*  NOTICE: This source code file is provided so that users may experiment
*  with enhancements to POV-Ray and to port the software to platforms other
*  than those supported by the POV-Ray Team.  There are strict rules under
*  which you are permitted to use this file.  The rules are in the file
*  named POVLEGAL.DOC which should be distributed with this file. If
*  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
*  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
*  Forum.  The latest version of POV-Ray may be found there as well.
*
* This program is based on the popular DKB raytracer version 2.12.
* DKBTrace was originally written by David K. Buck.
* DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
*
*****************************************************************************/

/****************************************************************************
*
*  Explanation:
*
*    -
*
*  Syntax:
*
*  ---
*
*  Aug 1994 : Merged functions for CSG height fields into functions for
*             non-CSG height fiels. Moved all height field map related
*             data into one data structure. Fixed memory problems. [DB]
*
*  Feb 1995 : Major rewrite of the height field intersection tests. [DB]
*
*****************************************************************************/

#include "frame.h"
#include "povray.h"
#include "vector.h"
#include "povproto.h"
#include "hfield.h"
#include "matrices.h"
#include "objects.h"



/*****************************************************************************
* Local preprocessor defines
******************************************************************************/

#define sign(x) (((x) >= 0.0) ? 1.0 : -1.0)

#define Get_Height(x, z, HField) ((DBL)(HField)->Data->Map[(z)][(x)])

/* Small offest. */

#define HFIELD_OFFSET 0.001



/*****************************************************************************
* Static functions
******************************************************************************/

static DBL normalize PARAMS((VECTOR A, VECTOR B));

static DBL stretch PARAMS((DBL x));

static void smooth_height_field PARAMS((HFIELD *HField, int xsize, int zsize));

static int intersect_pixel PARAMS((int x,int z,RAY *Ray,HFIELD *HField,DBL height1,DBL height2));

static int add_single_normal PARAMS((HF_VAL **data, int xsize, int zsize,
  int x0, int z0,int x1, int z1,int x2, int z2,VECTOR N));

static int  All_HField_Intersections PARAMS((OBJECT *Object,RAY *Ray,ISTACK *Depth_Stack));
static int  Inside_HField PARAMS((VECTOR IPoint,OBJECT *Object));
static void HField_Normal PARAMS((VECTOR Result,OBJECT *Object,INTERSECTION *Inter));
static void Translate_HField PARAMS((OBJECT *Object,VECTOR Vector, TRANSFORM *Trans));
static void Rotate_HField PARAMS((OBJECT *Object,VECTOR Vector, TRANSFORM *Trans));
static void Scale_HField PARAMS((OBJECT *Object,VECTOR Vector, TRANSFORM *Trans));
static void Invert_HField PARAMS((OBJECT *Object));
static void Transform_HField PARAMS((OBJECT *Object,TRANSFORM *Trans));
static void *Copy_HField PARAMS((OBJECT *Object));
static void Destroy_HField PARAMS((OBJECT *Object));

static int dda_traversal PARAMS((RAY *Ray, HFIELD *HField, VECTOR Start, HFIELD_BLOCK *Block));
static int block_traversal PARAMS((RAY *Ray, HFIELD *HField, VECTOR Start));
static void build_hfield_blocks PARAMS((HFIELD *HField));



/*****************************************************************************
* Local variables
******************************************************************************/

METHODS HField_Methods =
{
  All_HField_Intersections,
  Inside_HField, HField_Normal,
  Copy_HField, Translate_HField, Rotate_HField,
  Scale_HField, Transform_HField, Invert_HField, Destroy_HField
};

static ISTACK *HField_Stack;
static RAY *RRay;
static DBL mindist, maxdist;


/*****************************************************************************
*
* FUNCTION
*
*   All_HField_Intersections
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   Feb 1995 : Modified to work with new intersection functions. [DB]
*
******************************************************************************/

static int All_HField_Intersections(Object, Ray, Depth_Stack)
OBJECT *Object;
RAY *Ray;
ISTACK *Depth_Stack;
{
  int Side1, Side2;
  VECTOR Start;
  RAY Temp_Ray;
  DBL depth1, depth2;
  HFIELD *HField = (HFIELD *) Object;

  Increase_Counter(stats[Ray_HField_Tests]);

  MInvTransPoint(Temp_Ray.Initial, Ray->Initial, HField->Trans);
  MInvTransDirection(Temp_Ray.Direction, Ray->Direction, HField->Trans);

#ifdef HFIELD_EXTRA_STATS
  Increase_Counter(stats[Ray_HField_Box_Tests]);
#endif

  if (!Intersect_Box(&Temp_Ray,HField->bounding_box,&depth1,&depth2,&Side1,&Side2))
  {
    return(FALSE);
  }

#ifdef HFIELD_EXTRA_STATS
  Increase_Counter(stats[Ray_HField_Box_Tests_Succeeded]);
#endif

  HField->Data->cache_pos = 0;

  if (depth1 < EPSILON)
  {
    depth1 = EPSILON;

    if (depth1 > depth2)
    {
      return(FALSE);
    }
  }

  VEvaluateRay(Start, Temp_Ray.Initial, depth1, Temp_Ray.Direction);

  mindist = depth1;
  maxdist = depth2;

  HField_Stack = Depth_Stack;

  RRay = Ray;

  if (block_traversal(&Temp_Ray, HField, Start))
  {
    Increase_Counter(stats[Ray_HField_Tests_Succeeded]);

    return(TRUE);
  }
  else
  {
    return(FALSE);
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   Inside_HField
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   -
*
******************************************************************************/

static int Inside_HField (IPoint, Object)
VECTOR IPoint;
OBJECT *Object;
{
  HFIELD *HField = (HFIELD *) Object;
  int px, pz;
  DBL x,z,y1,y2,y3,water, dot1, dot2;
  VECTOR Local_Origin, H_Normal, Test;

  MInvTransPoint(Test, IPoint, HField->Trans);

  water = HField->bounding_box->bounds[0][Y];

  if ((Test[X] < 0.0) || (Test[X] >= HField->bounding_box->bounds[1][X]) ||
      (Test[Z] < 0.0) || (Test[Z] >= HField->bounding_box->bounds[1][Z]))
  {
    return(Test_Flag(HField, INVERTED_FLAG));
  }

  if (Test[Y] >= HField->bounding_box->bounds[1][Y])
  {
    return(Test_Flag(HField, INVERTED_FLAG));
  }

  if (Test[Y] < water)
  {
    return(!Test_Flag(HField, INVERTED_FLAG));
  }

  px = (int)Test[X];
  pz = (int)Test[Z];

  x = Test[X] - (DBL)px;
  z = Test[Z] - (DBL)pz;

  if ((x+z) < 1.0)
  {
    y1 = max(Get_Height(px,   pz,   HField), water);
    y2 = max(Get_Height(px+1, pz,   HField), water);
    y3 = max(Get_Height(px,   pz+1, HField), water);

    Make_Vector(Local_Origin,(DBL)px,y1,(DBL)pz);

    Make_Vector(H_Normal, y1-y2, 1.0, y1-y3);
  }
  else
  {
    px = (int)ceil(Test[X]);
    pz = (int)ceil(Test[Z]);

    y1 = max(Get_Height(px,   pz,   HField), water);
    y2 = max(Get_Height(px-1, pz,   HField), water);
    y3 = max(Get_Height(px,   pz-1, HField), water);

    Make_Vector(Local_Origin,(DBL)px,y1,(DBL)pz);

    Make_Vector(H_Normal, y2-y1, 1.0, y3-y1);
  }

  VDot(dot1, Test, H_Normal);
  VDot(dot2, Local_Origin, H_Normal);

  if (dot1 < dot2)
  {
    return(!Test_Flag(HField, INVERTED_FLAG));
  }

  return(Test_Flag(HField, INVERTED_FLAG));
}



/*****************************************************************************
*
* FUNCTION
*
*   HField_Normal
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   -
*
******************************************************************************/

static void HField_Normal (Result, Object, Inter)
OBJECT *Object;
VECTOR Result;
INTERSECTION *Inter;
{
  HFIELD *HField = (HFIELD *) Object;
  int px,pz, i;
  DBL x,z,y1,y2,y3,u,v;
  VECTOR Local_Origin;
  VECTOR n[5];

  MInvTransPoint(Local_Origin, Inter->IPoint, HField->Trans);

  for (i = 0; i < HField->Data->cache_pos; i++)
  {
    if ((Local_Origin[X] == HField->Data->Normal_Cache[i].fx) &&
       (Local_Origin[Z] == HField->Data->Normal_Cache[i].fz))
    {
      Assign_Vector(Result,HField->Data->Normal_Cache[i].normal);

      MTransNormal(Result,Result,HField->Trans);

      VNormalize(Result,Result);

      return;
    }
  }

  px = (int)Local_Origin[X];
  pz = (int)Local_Origin[Z];

  x = Local_Origin[X] - (DBL)px;
  z = Local_Origin[Z] - (DBL)pz;

  if (Test_Flag(HField, SMOOTHED_FLAG))
  {
    n[0][X] = HField->Data->Normals[pz][px][0];
    n[0][Y] = HField->Data->Normals[pz][px][1];
    n[0][Z] = HField->Data->Normals[pz][px][2];
    n[1][X] = HField->Data->Normals[pz][px+1][0];
    n[1][Y] = HField->Data->Normals[pz][px+1][1];
    n[1][Z] = HField->Data->Normals[pz][px+1][2];
    n[2][X] = HField->Data->Normals[pz+1][px][0];
    n[2][Y] = HField->Data->Normals[pz+1][px][1];
    n[2][Z] = HField->Data->Normals[pz+1][px][2];
    n[3][X] = HField->Data->Normals[pz+1][px+1][0];
    n[3][Y] = HField->Data->Normals[pz+1][px+1][1];
    n[3][Z] = HField->Data->Normals[pz+1][px+1][2];

    x = stretch(x);
    z = stretch(z);

    u = (1.0 - x);
    v = (1.0 - z);

    Result[X] = v*(u*n[0][X] + x*n[1][X]) + z*(u*n[2][X] + x*n[3][X]);
    Result[Y] = v*(u*n[0][Y] + x*n[1][Y]) + z*(u*n[2][Y] + x*n[3][Y]);
    Result[Z] = v*(u*n[0][Z] + x*n[1][Z]) + z*(u*n[2][Z] + x*n[3][Z]);
  }
  else
  {
    if ((x+z) <= 1.0)
    {
      /* Lower triangle. */

      y1 = Get_Height(px,   pz,   HField);
      y2 = Get_Height(px+1, pz,   HField);
      y3 = Get_Height(px,   pz+1, HField);

      Make_Vector(Result, y1-y2, 1.0, y1-y3);
    }
    else
    {
      /* Upper triangle. */

      y1 = Get_Height(px+1, pz+1, HField);
      y2 = Get_Height(px,   pz+1, HField);
      y3 = Get_Height(px+1, pz,   HField);

      Make_Vector(Result, y2-y1, 1.0, y3-y1);
    }
  }

  MTransNormal(Result, Result, HField->Trans);

  VNormalize(Result, Result);
}



/*****************************************************************************
*
* FUNCTION
*
*   stretch
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   -
*
******************************************************************************/

static DBL stretch (x)
DBL x;
{
  if (x <= 0.5)
  {
    x = 2.0 * x * x;
  }
  else
  {
    x = 1.0 - (2.0 * (1.0 - x) * (1.0 - x));
  }

  return(x);
}



/*****************************************************************************
*
* FUNCTION
*
*   normalize
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   -
*
******************************************************************************/

static DBL normalize(A, B)
VECTOR A, B;
{
  VLength(VTemp, B);

  if (fabs(VTemp) > EPSILON)
  {
    VInverseScale(A, B, VTemp);
  }
  else
  {
    Make_Vector(A, 0.0, 1.0, 0.0);
  }

  return(VTemp);
}



/*****************************************************************************
*
* FUNCTION
*
*   intersect_pixel
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   -
*
******************************************************************************/

static int intersect_pixel(x, z, Ray, HField, height1, height2)
int x;
int z;
RAY *Ray;
HFIELD *HField;
DBL height1, height2;
{
  int Found;
  DBL dot, depth1, depth2;
  DBL s, t, y1, y2, y3, y4;
  DBL min_y2_y3, max_y2_y3;
  DBL max_height, min_height;
  VECTOR P, N, V1;

#ifdef HFIELD_EXTRA_STATS
  Increase_Counter(stats[Ray_HField_Cell_Tests]);
#endif

  y1 = Get_Height(x,   z,   HField);
  y2 = Get_Height(x+1, z,   HField);
  y3 = Get_Height(x,   z+1, HField);
  y4 = Get_Height(x+1, z+1, HField);

  /* Do we hit this cell at all? */

  if (y2 < y3)
  {
    min_y2_y3 = y2;
    max_y2_y3 = y3;
  }
  else
  {
    min_y2_y3 = y3;
    max_y2_y3 = y2;
  }

  min_height = min(min(y1, y4), min_y2_y3);
  max_height = max(max(y1, y4), max_y2_y3);

  if ((max_height < height1) || (min_height > height2))
  {
    return(FALSE);
  }

#ifdef HFIELD_EXTRA_STATS
  Increase_Counter(stats[Ray_HField_Cell_Tests_Succeeded]);
#endif

  Found = FALSE;

  /* Check if we'll hit first triangle. */

  min_height = min(y1, min_y2_y3);
  max_height = max(y1, max_y2_y3);

  if ((max_height >= height1) && (min_height <= height2))
  {
#ifdef HFIELD_EXTRA_STATS
    Increase_Counter(stats[Ray_HField_Triangle_Tests]);
#endif

    /* Set up triangle. */

    Make_Vector(P, (DBL)x, y1, (DBL)z);

    /*
     * Calculate the normal vector from:
     *
     * N = V2 x V1, with V1 = <1, y2-y1, 0>, V2 = <0, y3-y1, 1>.
     */

    Make_Vector(N, y1-y2, 1.0, y1-y3);

    /* Now intersect the triangle. */

    VDot(dot, N, Ray->Direction);

    if ((dot > EPSILON) || (dot < -EPSILON))
    {
      VSub(V1, P, Ray->Initial);

      VDot(depth1, N, V1);

      depth1 /= dot;

      if ((depth1 >= mindist) && (depth1 <= maxdist))
      {
        s = Ray->Initial[X] + depth1 * Ray->Direction[X] - (DBL)x;
        t = Ray->Initial[Z] + depth1 * Ray->Direction[Z] - (DBL)z;

        if ((s >= -0.0001) && (t >= -0.0001) && ((s+t) <= 1.0001))
        {
#ifdef HFIELD_EXTRA_STATS
          Increase_Counter(stats[Ray_HField_Triangle_Tests_Succeeded]);
#endif

          VEvaluateRay(P, RRay->Initial, depth1, RRay->Direction);

          if (Point_In_Clip(P, HField->Clip))
          {
            push_entry(depth1, P, (OBJECT *)HField, HField_Stack);

            Found = TRUE;

            /* Cache normal. */

            if (!Test_Flag(HField, SMOOTHED_FLAG))
            {
              if (HField->Data->cache_pos < HF_CACHE_SIZE)
              {
                Assign_Vector(HField->Data->Normal_Cache[HField->Data->cache_pos].normal, N);

                HField->Data->Normal_Cache[HField->Data->cache_pos].fx = x + s;
                HField->Data->Normal_Cache[HField->Data->cache_pos].fz = z + t;

                HField->Data->cache_pos++;
              }
            }
          }
        }
      }
    }
  }

  /* Check if we'll hit second triangle. */

  min_height = min(y4, min_y2_y3);
  max_height = max(y4, max_y2_y3);

  if ((max_height >= height1) && (min_height <= height2))
  {
#ifdef HFIELD_EXTRA_STATS
    Increase_Counter(stats[Ray_HField_Triangle_Tests]);
#endif

    /* Set up triangle. */

    Make_Vector(P, (DBL)(x+1), y4, (DBL)(z+1));

    /*
     * Calculate the normal vector from:
     *
     * N = V2 x V1, with V1 = <-1, y3-y4, 0>, V2 = <0, y2-y4, -1>.
     */

    Make_Vector(N, y3-y4, 1.0, y2-y4);

    /* Now intersect the triangle. */

    VDot(dot, N, Ray->Direction);

    if ((dot > EPSILON) || (dot < -EPSILON))
    {
      VSub(V1, P, Ray->Initial);

      VDot(depth2, N, V1);

      depth2 /= dot;

      if ((depth2 >= mindist) && (depth2 <= maxdist))
      {
        s = Ray->Initial[X] + depth2 * Ray->Direction[X] - (DBL)x;
        t = Ray->Initial[Z] + depth2 * Ray->Direction[Z] - (DBL)z;

        if ((s <= 1.0001) && (t <= 1.0001) && ((s+t) >= 0.9999))
        {
#ifdef HFIELD_EXTRA_STATS
          Increase_Counter(stats[Ray_HField_Triangle_Tests_Succeeded]);
#endif

          VEvaluateRay(P, RRay->Initial, depth2, RRay->Direction);

          if (Point_In_Clip(P, HField->Clip))
          {
            push_entry(depth2, P, (OBJECT *)HField, HField_Stack);

            Found = TRUE;

            /* Cache normal. */

            if (!Test_Flag(HField, SMOOTHED_FLAG))
            {
              if (HField->Data->cache_pos < HF_CACHE_SIZE)
              {
                Assign_Vector(HField->Data->Normal_Cache[HField->Data->cache_pos].normal, N);

                HField->Data->Normal_Cache[HField->Data->cache_pos].fx = x + s;
                HField->Data->Normal_Cache[HField->Data->cache_pos].fz = z + t;

                HField->Data->cache_pos++;
              }
            }
          }
        }
      }
    }
  }

  return(Found);
}



/*****************************************************************************
*
* FUNCTION
*
*   add_single_normal
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   -
*
******************************************************************************/

static int add_single_normal(data, xsize, zsize, x0, z0, x1, z1, x2, z2, N)
HF_VAL **data;
int xsize,zsize,x0,z0,x1,z1,x2,z2;
VECTOR N;
{
  VECTOR v0, v1, v2, t0, t1, Nt;

  if ((x0 < 0) || (z0 < 0) ||
      (x1 < 0) || (z1 < 0) ||
      (x2 < 0) || (z2 < 0) ||
      (x0 > xsize) || (z0 > zsize) ||
      (x1 > xsize) || (z1 > zsize) ||
      (x2 > xsize) || (z2 > zsize))
  {
    return(0);
  }
  else
  {
    Make_Vector(v0, x0, (DBL)data[z0][x0], z0);
    Make_Vector(v1, x1, (DBL)data[z1][x1], z1);
    Make_Vector(v2, x2, (DBL)data[z2][x2], z2);

    VSub(t0, v2, v0);
    VSub(t1, v1, v0);

    VCross(Nt, t0, t1);

    normalize(Nt, Nt);

    if (Nt[Y] < 0.0)
    {
      VScaleEq(Nt, -1.0);
    }

    VAddEq(N, Nt);

    return(1);
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   smooth_height_field
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   Given a height field that only contains an elevation grid, this
*   routine will walk through the data and produce averaged normals
*   for all points on the grid.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void smooth_height_field(HField, xsize, zsize)
HFIELD *HField;
int xsize;
int zsize;
{
  int i, j, k;
  VECTOR N;
  HF_VAL **map = HField->Data->Map;

  /* First off, allocate all the memory needed to store the normal information */

  HField->Data->Normals_Height = zsize+1;

  HField->Data->Normals = (HF_Normals **)POV_MALLOC((zsize+1)*sizeof(HF_Normals *), "height field normals");

  for (i = 0; i <= zsize; i++)
  {
    HField->Data->Normals[i] = (HF_Normals *)POV_MALLOC((xsize+1)*sizeof(HF_Normals), "height field normals");
  }

  /*
   * For now we will do it the hard way - by generating the normals
   * individually for each elevation point.
   */

  for (i = 0; i < zsize; i++)
  {
    COOPERATE_1

    for (j = 0; j < xsize; j++)
    {
      Make_Vector(N, 0.0, 0.0, 0.0);

      k = 0;

      k += add_single_normal(map, xsize, zsize, j, i, j+1, i, j, i+1, N);
      k += add_single_normal(map, xsize, zsize, j, i, j, i+1, j-1, i, N);
      k += add_single_normal(map, xsize, zsize, j, i, j-1, i, j, i-1, N);
      k += add_single_normal(map, xsize, zsize, j, i, j, i-1, j+1, i, N);

      if (k == 0)
      {
        Error ("Failed to find any normals at: (%d, %d).\n", i, j);
      }

      normalize(N, N);

      HField->Data->Normals[i][j][0] = (short)(32767 * N[X]);
      HField->Data->Normals[i][j][1] = (short)(32767 * N[Y]);
      HField->Data->Normals[i][j][2] = (short)(32767 * N[Z]);
    }
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   Compute_HField
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*
* DESCRIPTION
*
*   Copy image data into height field map. Create bounding blocks
*   for the block traversal. Calculate normals for smoothed height fields.
*
* CHANGES
*
*   Feb 1995 : Modified to work with new intersection functions. [DB]
*
******************************************************************************/

void Compute_HField(HField, Image)
HFIELD *HField;
IMAGE *Image;
{
  int x, z, max_x, max_z;
  int temp1 = 0, temp2 = 0;
  HF_VAL min_y, max_y, temp_y;

  /* Get height field map size. */

  max_x = Image->iwidth;

  if (Image->File_Type == POT_FILE)
  {
    max_x = max_x / 2;
  }

  max_z = Image->iheight;

  /* Allocate memory for map. */

  HField->Data->Map = (HF_VAL **)POV_MALLOC(max_z*sizeof(HF_VAL *), "height field");

  for (z = 0; z < max_z; z++)
  {
    HField->Data->Map[z] = (HF_VAL *)POV_MALLOC(max_x*sizeof(HF_VAL), "height field");
  }

  /* Copy map. */

  min_y = 65535L;
  max_y = 0;

  for (z = 0; z < max_z; z++)
  {
    COOPERATE_1
    for (x = 0; x < max_x; x++)
    {

      switch (Image->File_Type)
      {
        case GIF_FILE:

          temp1 = Image->data.map_lines[max_z - z - 1][x];
          temp2 = 0;

          break;

        case POT_FILE:

          temp1 = Image->data.map_lines[max_z - z - 1][x];
          temp2 = Image->data.map_lines[max_z - z - 1][x + max_x];

          break;


        case PPM_FILE:

          temp1 = Image->data.rgb_lines[max_z - z - 1].red[x];
          temp2 = Image->data.rgb_lines[max_z - z - 1].green[x];

          break;

        case PGM_FILE:
        case TGA_FILE:
        case PNG_FILE:

          if (Image->Colour_Map == NULL)
          {
            temp1 = Image->data.rgb_lines[max_z - z - 1].red[x];
            temp2 = Image->data.rgb_lines[max_z - z - 1].green[x];
          }
          else
          {
            temp1 = Image->data.map_lines[max_z - z - 1][x];
            temp2 = 0;
          }

          break;

        default:

          Error("Unknown image type in Compute_HField().\n");
      }

      temp_y = (HF_VAL)(256*temp1 + temp2);

      HField->Data->Map[z][x] = temp_y;

      min_y = min(min_y, temp_y);
      max_y = max(max_y, temp_y);
    }
  }

  /* Resize bounding box. */

  HField->Data->min_y = min_y;
  HField->Data->max_y = max_y;

  HField->bounding_box->bounds[0][Y] = max((DBL)min_y, HField->bounding_box->bounds[0][Y]) - HFIELD_OFFSET;
  HField->bounding_box->bounds[1][Y] = (DBL)max_y + HFIELD_OFFSET;

  /* Compute smoothed height field. */

  if (Test_Flag(HField, SMOOTHED_FLAG))
  {
    smooth_height_field(HField, max_x-1, max_z-1);
  }

  HField->Data->max_x = max_x-2;
  HField->Data->max_z = max_z-2;

  build_hfield_blocks(HField);
}



/*****************************************************************************
*
* FUNCTION
*
*   build_hfield_blocks
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Create the bounding block hierarchy used by the block traversal.
*
* CHANGES
*
*   Feb 1995 : Creation.
*
******************************************************************************/

static void build_hfield_blocks(HField)
HFIELD *HField;
{
  int x, z, nx, nz, wx, wz;
  int i, j;
  int xmin, xmax, zmin, zmax;
  DBL y, ymin, ymax, water;

  /* Get block size. */

  nx = max(1, (int)(sqrt((DBL)HField->Data->max_x)));
  nz = max(1, (int)(sqrt((DBL)HField->Data->max_z)));

  /* Get dimensions of sub-block. */

  wx = (int)ceil((DBL)(HField->Data->max_x + 2) / (DBL)nx);
  wz = (int)ceil((DBL)(HField->Data->max_z + 2) / (DBL)nz);

  /* Increase number of sub-blocks if necessary. */

  if (nx * wx < HField->Data->max_x + 2)
  {
    nx++;
  }

  if (nz * wz < HField->Data->max_z + 2)
  {
    nz++;
  }

  if (!Test_Flag(HField, HIERARCHY_FLAG) || ((nx == 1) && (nz == 1)))
  {
    /* We don't want a bounding hierarchy. Just use one block. */

    HField->Data->Block = (HFIELD_BLOCK **)POV_MALLOC(sizeof(HFIELD_BLOCK *), "height field blocks");

    HField->Data->Block[0] = (HFIELD_BLOCK *)POV_MALLOC(sizeof(HFIELD_BLOCK), "height field blocks");

    HField->Data->Block[0][0].xmin = 0;
    HField->Data->Block[0][0].xmax = HField->Data->max_x;
    HField->Data->Block[0][0].zmin = 0;
    HField->Data->Block[0][0].zmax = HField->Data->max_z;

    HField->Data->Block[0][0].ymin = HField->bounding_box->bounds[0][Y];
    HField->Data->Block[0][0].ymax = HField->bounding_box->bounds[1][Y];

    HField->Data->block_max_x = 1;
    HField->Data->block_max_z = 1;

    HField->Data->block_width_x = HField->Data->max_x + 2;
    HField->Data->block_width_z = HField->Data->max_y + 2;

/*
    Debug_Info("\nHeight field: %d x %d (1 x 1 blocks)", HField->Data->max_x+2, HField->Data->max_z+2);
*/

    return;
  }

/*
  Debug_Info("\nHeight field: %d x %d (%d x %d blocks)", HField->Data->max_x+2, HField->Data->max_z+2, nx, nz);
*/

  /* Allocate memory for blocks. */

  HField->Data->Block = (HFIELD_BLOCK **)POV_MALLOC(nz*sizeof(HFIELD_BLOCK *), "height field blocks");

  /* Store block information. */

  HField->Data->block_max_x = nx;
  HField->Data->block_max_z = nz;

  HField->Data->block_width_x = wx;
  HField->Data->block_width_z = wz;

  water = HField->bounding_box->bounds[0][Y];

  for (z = 0; z < nz; z++)
  {
    COOPERATE_1

    HField->Data->Block[z] = (HFIELD_BLOCK *)POV_MALLOC(nx*sizeof(HFIELD_BLOCK), "height field blocks");

    for (x = 0; x < nx; x++)
    {
      /* Get block's borders. */

      xmin = x * wx;
      zmin = z * wz;

      xmax = min((x + 1) * wx - 1, HField->Data->max_x);
      zmax = min((z + 1) * wz - 1, HField->Data->max_z);

      /* Find min. and max. height in current block. */

      ymin = BOUND_HUGE;
      ymax = -BOUND_HUGE;

      for (i = xmin; i <= xmax+1; i++)
      {
        for (j = zmin; j <= zmax+1; j++)
        {
          y = Get_Height(i, j, HField);

          ymin = min(ymin, y);
          ymax = max(ymax, y);
        }
      }

      /* Store block's borders. */

      HField->Data->Block[z][x].xmin = xmin;
      HField->Data->Block[z][x].xmax = xmax;
      HField->Data->Block[z][x].zmin = zmin;
      HField->Data->Block[z][x].zmax = zmax;

      HField->Data->Block[z][x].ymin = max(ymin, water) - HFIELD_OFFSET;
      HField->Data->Block[z][x].ymax = ymax + HFIELD_OFFSET;
    }
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   Translate_HField
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Translate_HField (Object, Vector, Trans)
OBJECT *Object;
VECTOR Vector;
TRANSFORM *Trans;
{
  Transform_HField(Object, Trans);
}



/*****************************************************************************
*
* FUNCTION
*
*   Rotate_HField
*
* INPUT
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Rotate_HField (Object, Vector, Trans)
OBJECT *Object;
VECTOR Vector;
TRANSFORM *Trans;
{
  Transform_HField(Object, Trans);
}



/*****************************************************************************
*
* FUNCTION
*
*   Scale_HField
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Scale_HField (Object, Vector, Trans)
OBJECT *Object;
VECTOR Vector;
TRANSFORM *Trans;
{
  Transform_HField(Object, Trans);
}



/*****************************************************************************
*
* FUNCTION
*
*   Invert_HField
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Invert_HField (Object)
OBJECT *Object;
{
  Invert_Flag(Object, INVERTED_FLAG);
}



/*****************************************************************************
*
* FUNCTION
*
*   Transform_HField
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   -
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Transform_HField (Object, Trans)
OBJECT *Object;
TRANSFORM *Trans;
{
  Compose_Transforms(((HFIELD *)Object)->Trans, Trans);

  Compute_HField_BBox((HFIELD *)Object);
}



/*****************************************************************************
*
* FUNCTION
*
*   Create_HField
*
* INPUT
*
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   Allocate and intialize a height field.
*
* CHANGES
*
*   Feb 1995 : Modified to work with new intersection functions. [DB]
*
******************************************************************************/

HFIELD *Create_HField()
{
  HFIELD *New;

  /* Allocate height field. */

  New = (HFIELD *)POV_MALLOC(sizeof(HFIELD), "height field");

  INIT_OBJECT_FIELDS(New, HFIELD_OBJECT, &HField_Methods)

  /* Always uses Trans so always create one. */

  New->Trans = Create_Transform();

  New->bounding_box = Create_Box();

  /* Allocate height field data. */

  New->Data = (HFIELD_DATA *)POV_MALLOC(sizeof(HFIELD_DATA), "height field");

  New->Data->References = 1;

  New->Data->cache_pos = 0;

  New->Data->Normals_Height = 0;

  New->Data->Map     = NULL;
  New->Data->Normals = NULL;

  New->Data->max_x = 0;
  New->Data->max_z = 0;

  New->Data->block_max_x = 0;
  New->Data->block_max_z = 0;

  New->Data->block_width_x = 0;
  New->Data->block_width_z = 0;

  Set_Flag(New, HIERARCHY_FLAG);

  return(New);
}



/*****************************************************************************
*
* FUNCTION
*
*   Copy_HField
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   NOTE: The height field data is not copied, only the number of references
*         is counted, so that Destray_HField() knows if it can be destroyed.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void *Copy_HField(Object)
OBJECT *Object;
{
  HFIELD *New;

  New = Create_HField();

  /* Destroy obsolete things created in Create_HField(). */

  Destroy_Transform(New->Trans);

  Destroy_Box((OBJECT *)(New->bounding_box));

  POV_FREE (New->Data);

  /* Copy height field. */

  *New = *((HFIELD *)Object);

  New->Trans = Copy_Transform(((HFIELD *)Object)->Trans);

  New->bounding_box = Copy_Box((OBJECT *)(((HFIELD *)Object)->bounding_box));

  New->Data->References++;

  return(New);
}



/*****************************************************************************
*
* FUNCTION
*
*   Destroy_HField
*
* INPUT
*   
* OUTPUT
*   
* RETURNS
*   
* AUTHOR
*
*   Doug Muir, David Buck, Drew Wells
*   
* DESCRIPTION
*
*   NOTE: The height field data is destroyed if it's no longer
*         used by any copy.
*
* CHANGES
*
*   Feb 1995 : Modified to work with new intersection functions. [DB]
*
******************************************************************************/

static void Destroy_HField (Object)
OBJECT *Object;
{
  int i;
  HFIELD *HField = (HFIELD *)Object;

  Destroy_Transform(HField->Trans);

  Destroy_Box((OBJECT *)(HField->bounding_box));

  if (--(HField->Data->References) == 0)
  {
    if (HField->Data->Map != NULL)
    {
      for (i = 0; i < HField->Data->max_z+2; i++)
      {
        if (HField->Data->Map[i] != NULL)
        {
          POV_FREE (HField->Data->Map[i]);
        }
      }

      POV_FREE (HField->Data->Map);
    }

    if (HField->Data->Normals != NULL)
    {
      for (i = 0; i < HField->Data->Normals_Height; i++)
      {
        POV_FREE (HField->Data->Normals[i]);
      }

      POV_FREE (HField->Data->Normals);
    }

    if (HField->Data->Block != NULL)
    {
      for (i = 0; i < HField->Data->block_max_z; i++)
      {
        POV_FREE(HField->Data->Block[i]);
      }

      POV_FREE(HField->Data->Block);
    }

    POV_FREE (HField->Data);
  }

  POV_FREE (Object);
}



/*****************************************************************************
*
* FUNCTION
*
*   Compute_HField_BBox
*
* INPUT
*
*   HField - Height field
*   
* OUTPUT
*
*   HField
*   
* RETURNS
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Calculate the bounding box of a height field.
*
* CHANGES
*
*   Aug 1994 : Creation.
*
******************************************************************************/

void Compute_HField_BBox(HField)
HFIELD *HField;
{
  Assign_BBox_Vect(HField->BBox.Lower_Left, HField->bounding_box->bounds[0]);

  VSub (HField->BBox.Lengths, HField->bounding_box->bounds[1], HField->bounding_box->bounds[0]);

  if (HField->Trans != NULL)
  {
    Recompute_BBox(&HField->BBox, HField->Trans);
  }
}



/*****************************************************************************
*
* FUNCTION
*
*   dda_traversal
*
* INPUT
*
*   Ray    - Current ray
*   HField - Height field
*   Start  - Start point for the walk
*   Block  - Sub-block of the height field to traverse
*   
* OUTPUT
*   
* RETURNS
*
*   int - TRUE if intersection was found
*   
* AUTHOR
*
*   Dieter Bayer
*   
* DESCRIPTION
*
*   Traverse the grid cell of the height field using a modified DDA.
*
*   Based on the following article:
*
*     Musgrave, F. Kenton, "Grid Tracing: Fast Ray Tracing for Height
*     Fields", Research Report YALEU-DCS-RR-39, Yale University, July 1988
*
*   You should note that there are (n-1) x (m-1) grid cells in a height
*   field of (image) size n x m. A grid cell (i,j), 0 <= i <= n-1,
*   0 <= j <= m-1, extends from x = i to x = i + 1 - epsilon and
*   y = j to y = j + 1 -epsilon.
*
* CHANGES
*
*   Feb 1995 : Creation.
*
******************************************************************************/

static int dda_traversal(Ray, HField, Start, Block)
RAY *Ray;
HFIELD *HField;
VECTOR Start;
HFIELD_BLOCK *Block;
{
  int found;
  int xmin, xmax, zmin, zmax;
  int x, z, signx, signz;
  DBL ymin, ymax, y1, y2;
  DBL px, pz, dx, dy, dz;
  DBL delta, error, x0, z0;
  DBL neary, fary, deltay;

  /* Setup DDA. */

  found = FALSE;

  px = Start[X];
  pz = Start[Z];

  /* Get dimensions of current block. */

  xmin = Block->xmin;
  xmax = min(Block->xmax + 1, HField->Data->max_x);
  zmin = Block->zmin;
  zmax = min(Block->zmax + 1, HField->Data->max_z);

  ymin = min(Start[Y], Block->ymin) - EPSILON;
  ymax = max(Start[Y], Block->ymax) + EPSILON;

  /* Check for illegal grid values (caused by numerical inaccuracies). */

  if (px < (DBL)xmin)
  {
    if (px < (DBL)xmin - HFIELD_OFFSET)
    {
      Debug_Info("Illegal grid value in dda_traversal().\n");

      return(FALSE);
    }
    else
    {
      px = (DBL)xmin;
    }
  }
  else
  {
    if (px > (DBL)xmax + 1.0 - EPSILON)
    {
      if (px > (DBL)xmax + 1.0 + EPSILON)
      {
        Debug_Info("Illegal grid value in dda_traversal().\n");

        return(FALSE);
      }
      else
      {
        px = (DBL)xmax + 1.0 - EPSILON;
      }
    }
  }

  if (pz < (DBL)zmin)
  {
    if (pz < (DBL)zmin - HFIELD_OFFSET)
    {
      Debug_Info("Illegal grid value in dda_traversal().\n");

      return(FALSE);
    }
    else
    {
      pz = (DBL)zmin;
    }
  }
  else
  {
    if (pz > (DBL)zmax + 1.0 - EPSILON)
    {
      if (pz > (DBL)zmax + 1.0 + EPSILON)
      {
        Debug_Info("Illegal grid value in dda_traversal().\n");

        return(FALSE);
      }
      else
      {
        pz = (DBL)zmax + 1.0 - EPSILON;
      }
    }
  }

  dx = Ray->Direction[X];
  dy = Ray->Direction[Y];
  dz = Ray->Direction[Z];

  /*
   * Here comes the DDA algorithm.
   */

  /* Choose algorithm depending on the driving axis. */

  if (fabs(dx) >= fabs(dz))
  {
    /*
     * X-axis is driving axis.
     */

    delta = fabs(dz / dx);

    x = (int)px;
    z = (int)pz;

    x0 = px - floor(px);
    z0 = pz - floor(pz);

    signx = sign(dx);
    signz = sign(dz);

    /* Get initial error. */

    if (dx >= 0.0)
    {
      if (dz >= 0.0)
      {
        error = z0 + delta * (1.0 - x0) - 1.0;
      }
      else
      {
        error = -(z0 - delta * (1.0 - x0));
      }
    }
    else
    {
      if (dz >= 0.0)
      {
        error = z0 + delta * x0 - 1.0;
      }
      else
      {
        error = -(z0 - delta * x0);
      }
    }

    /* Get y differential. */

    deltay = dy / fabs(dx);

    if (dx >= 0.0)
    {
      neary = Start[Y] - x0 * deltay;

      fary = neary + deltay;
    }
    else
    {
      neary = Start[Y] - (1.0 - x0) * deltay;

      fary = neary + deltay;
    }

    /* Step through the cells. */

    do
    {
      if (neary < fary)
      {
        y1 = neary;
        y2 = fary;
      }
      else
      {
        y1 = fary;
        y2 = neary;
      }

      if (intersect_pixel(x, z, Ray, HField, y1, y2))
      {
        if (HField->Type & IS_CHILD_OBJECT)
        {
          found = TRUE;
        }
        else
        {
          return(TRUE);
        }
      }

      if (error > EPSILON)
      {
        z += signz;

        if ((z < zmin) || (z > zmax))
        {
          break;
        }
        else
        {
          if (intersect_pixel(x, z, Ray, HField, y1, y2))
          {
            if (HField->Type & IS_CHILD_OBJECT)
            {
              found = TRUE;
            }
            else
            {
              return(TRUE);
            }
          }
        }

        error--;
      }
      else
      {
        if (error > -EPSILON)
        {
          z += signz;

          error--;
        }
      }

      x += signx;

      error += delta;

      neary = fary;

      fary += deltay;
    }
    while ((neary >= ymin) && (neary <= ymax) && (x >= xmin) && (x <= xmax) && (z >= zmin) && (z <= zmax));
  }
  else
  {
    /*
     * Z-axis is driving axis.
     */

    delta = fabs(dx / dz);

    x = (int)px;
    z = (int)pz;

    x0 = px - floor(px);
    z0 = pz - floor(pz);

    signx = sign(dx);
    signz = sign(dz);

    /* Get initial error. */

    if (dz >= 0.0)
    {
      if (dx >= 0.0)
      {
        error = x0 + delta * (1.0 - z0) - 1.0;
      }
      else
      {
        error = -(x0 - delta * (1.0 - z0));
      }
    }
    else
    {
      if (dx >= 0.0)
      {
        error = x0 + delta * z0 - 1.0;
      }
      else
      {
        error = -(x0 - delta * z0);
      }
    }

    /* Get y differential. */

    deltay = dy / fabs(dz);

    if (dz >= 0.0)
    {
      neary = Start[Y] - z0 * deltay;

      fary = neary + deltay;
    }
    else
    {
      neary = Start[Y] - (1.0 - z0) * deltay;

      fary = neary + deltay;
    }

    /* Step through the cells. */

    do
    {
      if (neary < fary)
      {
        y1 = neary;
        y2 = fary;
      }
      else
      {
        y1 = fary;
        y2 = neary;
      }

      if (intersect_pixel(x, z, Ray, HField, y1, y2))
      {
        if (HField->Type & IS_CHILD_OBJECT)
        {
          found = TRUE;
        }
        else
        {
          return(TRUE);
        }
      }

      if (error > EPSILON)
      {
        x += signx;

        if ((x < xmin) || (x > xmax))
        {
          break;
        }
        else
        {
          if (intersect_pixel(x, z, Ray, HField, y1, y2))
          {
            if (HField->Type & IS_CHILD_OBJECT)
            {
              found = TRUE;
            }
            else
            {
              return(TRUE);
            }
          }
        }

        error--;
      }
      else
      {
        if (error > -EPSILON)
        {
          x += signx;

          error--;
        }
      }

      z += signz;

      error += delta;

      neary = fary;

      fary += deltay;
    }
    while ((neary >= ymin-EPSILON) && (neary <= ymax+EPSILON) &&
           (x >= xmin) && (x <= xmax) &&
           (z >= zmin) && (z <= zmax));
  }

  return(found);
}



/*****************************************************************************
*
* FUNCTION
*
*   block_traversal
*
* INPUT
*
*   Ray    - Current ray
*   HField - Height field
*   Start  - Start point for the walk
*
* OUTPUT
*
* RETURNS
*
*   int - TRUE if intersection was found
*   
* AUTHOR
*
*   Dieter Bayer
*
* DESCRIPTION
*
*   Traverse the blocks of the height field.
*
* CHANGES
*
*   Feb 1995 : Creation.
*
*   Aug 1996 : Fixed bug as reported by Dean M. Phillips:
*              "I found a bug in the height field code which resulted
*              in "Illegal grid value in dda_traversal()." messages
*              along with dark vertical lines in the height fields.
*              This turned out to be caused by overlooking the units in
*              which some boundary tests in two different places were
*              made. It was easy to fix.
*
******************************************************************************/

static int block_traversal(Ray, HField, Start)
RAY *Ray;
HFIELD *HField;
VECTOR Start;
{
  int xmax, zmax;
  int x, z, nx, nz, signx, signz;
  int found = FALSE;
  int dx_zero, dz_zero;
  DBL px, pz, dx, dy, dz;
  DBL maxdv;
  DBL ymin, ymax, y1, y2;
  DBL neary, fary;
  DBL k1, k2, dist;
  VECTOR nearP, farP;
  HFIELD_BLOCK *Block;

  px = Start[X];
  pz = Start[Z];

  dx = Ray->Direction[X];
  dy = Ray->Direction[Y];
  dz = Ray->Direction[Z];

  maxdv = (dx > dz) ? dx : dz;

  /* First test for 'perpendicular' rays. */

  if ((fabs(dx) < EPSILON) && (fabs(dz) < EPSILON))
  {
    x = (int)px;
    z = (int)pz;

    neary = Start[Y];

    if (dy >= 0.0)
    {
      fary = 65536.0;
    }
    else
    {
      fary = 0.0;
    }

    return(intersect_pixel(x, z, Ray, HField, min(neary, fary), max(neary, fary)));
  }

  /* If we don't have blocks we just step through the grid. */

  if ((HField->Data->block_max_x <= 1) && (HField->Data->block_max_z <= 1))
  {
    return(dda_traversal(Ray, HField, Start, &HField->Data->Block[0][0]));
  }

  /* Get dimensions of grid. */

  xmax = HField->Data->block_max_x;
  zmax = HField->Data->block_max_z;

  ymin = (DBL)HField->Data->min_y - EPSILON;
  ymax = (DBL)HField->Data->max_y + EPSILON;

  dx_zero = (fabs(dx) < EPSILON);
  dz_zero = (fabs(dz) < EPSILON);

  signx = sign(dx);
  signz = sign(dz);

  /* Walk on the block grid. */

  px /= HField->Data->block_width_x;
  pz /= HField->Data->block_width_z;

  x = (int)px;
  z = (int)pz;

  Assign_Vector(nearP, Start);

  neary = Start[Y];

  /*
   * Here comes the block walk algorithm.
   */

  do
  {
#ifdef HFIELD_EXTRA_STATS
    Increase_Counter(stats[Ray_HField_Block_Tests]);
#endif

    /* Get current block. */

    Block = &HField->Data->Block[z][x];

    /* Intersect ray with bounding planes. */

    if (dx_zero)
    {
      k1 = BOUND_HUGE;
    }
    else
    {
      if (signx >= 0)
      {
        k1 = ((DBL)Block->xmax + 1.0 - Ray->Initial[X]) / dx;
      }
      else
      {
        k1 = ((DBL)Block->xmin - Ray->Initial[X]) / dx;
      }
    }

    if (dz_zero)
    {
      k2 = BOUND_HUGE;
    }
    else
    {
      if (signz >= 0)
      {
        k2 = ((DBL)Block->zmax + 1.0 - Ray->Initial[Z]) / dz;
      }
      else
      {
        k2 = ((DBL)Block->zmin - Ray->Initial[Z]) / dz;
      }
    }

    /* Figure out the indices of the next block. */

    if ((k1 < k2 - EPSILON / maxdv) && (k1 > 0.0))
    {
      /* Step along the x-axis. */

      dist = k1;

      nx = x + signx;
      nz = z;
    }
    else
    {
      if ((k1 < k2 + EPSILON / maxdv) && (k1 > 0.0))
      {
        /* Step along both axis (very rare case). */

        dist = k1;

        nx = x + signx;
        nz = z + signz;
      }
      else
      {
        /* Step along the z-axis. */

        dist = k2;

        nx = x;
        nz = z + signz;
      }
    }

    /* Get point where ray leaves current block. */

    VEvaluateRay(farP, Ray->Initial, dist, Ray->Direction);

    fary = farP[Y];

    if (neary < fary)
    {
      y1 = neary;
      y2 = fary;
    }
    else
    {
      y1 = fary;
      y2 = neary;
    }

    /* Can we hit current block at all? */

    if ((y1 <= (DBL)Block->ymax + EPSILON) && (y2 >= (DBL)Block->ymin - EPSILON))
    {
      /* Test current block. */

#ifdef HFIELD_EXTRA_STATS
      Increase_Counter(stats[Ray_HField_Block_Tests_Succeeded]);
#endif

      if (dda_traversal(Ray, HField, nearP, &HField->Data->Block[z][x]))
      {
        if (HField->Type & IS_CHILD_OBJECT)
        {
          found = TRUE;
        }
        else
        {
          return(TRUE);
        }
      }
    }

    /* Step to next block. */

    x = nx;
    z = nz;

    Assign_Vector(nearP, farP);

    neary = fary;
  }
  while ((x >= 0) && (x < xmax) && (z >= 0) && (z < zmax) && (neary >= ymin) && (neary <= ymax));

  return(found);
}