/****************************************************************************
*
*  ATTENTION!!!
*
*  THIS FILE HAS BEEN MODIFIED!!! IT IS NOT PART OF THE OFFICAL
*  POV-RAY 2.2 DISTRIBUTION!!!
*
*  THIS FILE IS PART OF "FASTER THAN POV-RAY" (VERSION 1.1),
*  A SPED-UP VERSION OF POV-RAY 2.2. USE AT YOUR OWN RISK!!!!!!
*
*  New files: addon0.c, addon1.c, addon2.c, addon3.c, addon.h
*
*  The additional modules were written by Dieter Bayer.
*
*  Send comments, suggestions, bugs, ideas ... to:
*
*  dieter@cip.e-technik.uni-erlangen.de
*
*  All changed/added lines are enclosed in #ifdef DB_CODE ... #endif
*
*  The vista projection was taken from:
*
*    A. Hashimoto, T. Akimoto, K. Mase, and Y. Suenaga, 
*    "Vista Ray-Tracing: High Speed Ray Tracing Using Perspective
*    Projection Image", New Advances in Computer Graphics, Proceedings
*    of CG International '89, R. A. Earnshaw, B. Wyvill (Eds.), 
*    Springer, ..., pp. 549-560
*
*  The idea for the light buffer was taken from:
*
*    E. Haines and D. Greenberg, "The Light Buffer: A Shadow-Testing 
*    Accelerator", IEEE CG&A, Vol. 6, No. 9, Sept. 1986, pp. 6-16
*
*****************************************************************************/

/****************************************************************************
*  addon1.c
*
*  This module was written by Dieter Bayer.
*
*  This module contains functions used only for the vista buffer.
*
*  01.03.1994 : Creation
*
*  29.04.1994 : Version 2.0
*
******************************************************************************/

#include <time.h>
#include "frame.h"
#include "vector.h"
#include "povproto.h"
#include "addon.h"

#ifdef DB_CODE

/* count project tests */
#define COUNT_PROJECT_TESTS

#define NUMBER_OF_TYPES      18
#define TYPE_BICUBIC_PATCH    0
#define TYPE_BLOB             1
#define TYPE_BOX              2
#define TYPE_CONE             3
#define TYPE_CSG_INTERSECTION 4
#define TYPE_CSG_MERGE        5
#define TYPE_CSG_UNION        6
#define TYPE_DISC             7
#define TYPE_ELLIPSOID        8
#define TYPE_HEIGHT_FIELD     9
#define TYPE_LIGHT_SOURCE    10
#define TYPE_PLANE           11
#define TYPE_POLY            12
#define TYPE_QUADRIC         13
#define TYPE_SMOOTH_TRIANGLE 14
#define TYPE_SPHERE          15
#define TYPE_TRIANGLE        16
#define TYPE_UNKNOWN         17

#define NUMBER_OF_FLAGS     9
#define FLAG_SUM            0
#define FLAG_INFINITE       1
#define FLAG_SCREEN_FILLING 2
#define FLAG_USED_IN_CSG    3
#define FLAG_INVISIBLE      4
#define FLAG_BOUNDED        5
#define FLAG_CLIPPED        6
#define FLAG_BOUND_OBJECT   7
#define FLAG_CLIP_OBJECT    8



/*****************************************************************************
* external variables
******************************************************************************/

extern long Number_Of_Rays;
extern FRAME Frame;
extern int Trace_Level;
extern int Use_Slabs;
extern DBL Max_Trace_Level;
extern time_t tstart, tstop;
extern unsigned int Options;
extern OBJECT *Root_Object;

extern METHODS Bicubic_Patch_Methods;
extern METHODS Blob_Methods;
extern METHODS Box_Methods;
extern METHODS Cone_Methods;
extern METHODS Csg_Height_Field_Methods;
extern METHODS CSG_Intersection_Methods;
extern METHODS CSG_Merge_Methods;
extern METHODS CSG_Union_Methods;
extern METHODS Disc_Methods;
extern METHODS Ellipsoid_Methods;
extern METHODS Height_Field_Methods;
extern METHODS Light_Source_Methods;
extern METHODS Plane_Methods;
extern METHODS Poly_Methods;
extern METHODS Quadric_Methods;
extern METHODS Smooth_Triangle_Methods;
extern METHODS Sphere_Methods;
extern METHODS Triangle_Methods;

extern unsigned MAXQUEUE;
extern unsigned long totalQueues, totalQueueResets, nChecked, nEnqueued;

extern size_t Mem_Vista_Buffer;



/*****************************************************************************
* Non-static variables
******************************************************************************/

PROJECT_TREE_NODE *Root_Vista;
unsigned long Project_Tests = 0, Project_Tests_Succeeded = 0;
unsigned Project_Algorithm = PROJECTIONS_WITH_LEAF_SLABS;




/*****************************************************************************
* Static variables
******************************************************************************/

static DBL distance;
static MATRIX WC2VC, WC2VCinv;
static VECTOR O, U, V, W;

/* Planes for 3d-clipping */

static VECTOR VIEW_VX1 = {-0.8944271910, 0.0, -0.4472135955};
static VECTOR VIEW_VX2 = { 0.8944271910, 0.0, -0.4472135955};
static VECTOR VIEW_VY1 = {0.0, -0.8944271910, -0.4472135955};
static VECTOR VIEW_VY2 = {0.0,  0.8944271910, -0.4472135955};
static DBL VIEW_DX1 = 0.4472135955;
static DBL VIEW_DX2 = 0.4472135955;
static DBL VIEW_DY1 = 0.4472135955;
static DBL VIEW_DY2 = 0.4472135955;

/* Queues */

static Qelem *Tree_Leaf_Queue;
static PROJECT_TREE_NODE **Tree_Node_Queue;



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

static void Project_rectangle PARAMS((PROJECT *Project, VECTOR P1, VECTOR P2, VECTOR P3, VECTOR P4, int *visible));
static void Project_triangle PARAMS((PROJECT *Project, VECTOR P1, VECTOR P2, VECTOR P3, int *visible));

static void Project_Bbox PARAMS((PROJECT *Project, VECTOR *Points, int *visible));
static void Project_Bounds PARAMS((PROJECT *Project, BBOX *Bounds, int *visible));

static void Get_Perspective_Projection PARAMS((OBJECT *Object, PROJECT *Project, int infinite));

static void Project_Object PARAMS((OBJECT *Object, PROJECT *Project));

static void Project_Bicubic_Patch PARAMS((PROJECT *Project, OBJECT *Object));
static void Project_Box PARAMS((PROJECT *Project, OBJECT *Object, int *visible));
static void Project_Height_Field PARAMS((PROJECT *Project, OBJECT *Object, int *visible));
static void Project_Sphere PARAMS((PROJECT *Project, OBJECT *Object));
static void Project_Triangle PARAMS((PROJECT *Project, OBJECT *Object, int *visible));
static void Project_Smooth_Triangle PARAMS((PROJECT *Project, OBJECT *Object, int *visible));

static void Get_Ellipse_Projection PARAMS((PROJECT *Project, DBL a20, DBL a02, DBL a11, DBL a10, DBL a01, DBL a00));

static void Transform_Point PARAMS((VECTOR *P));

static void Project_Bounding_Slab PARAMS((PROJECT *Project, PROJECT_TREE_NODE **Tree, OBJECT *Object));

static void Create_Rayinfo PARAMS((RAY *Ray, RAYINFO *rayinfo));

static int Intersect_Vista_Tree PARAMS((RAY *Ray, PROJECT_TREE_NODE *Tree, int x, INTERSECTION *Best_Intersection));




/*****************************************************************************
*
* FUNCTION      : Init_Project_Tree_Queues
*
* ARGUMENTS     : none
*
* MODIFIED ARGS : none
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Init queues for descending project tress, i.e. allocate memory needed.
*
* CHANGES
*
*   -
*
******************************************************************************/

void Init_Project_Tree_Queues()
{
  Tree_Leaf_Queue = (Qelem *)VB_malloc(MAXQUEUE*sizeof(Qelem));

  if (Tree_Leaf_Queue == NULL)
  {
    Fatal_MAError("priority queue for vista/light buffer");
  }

  Tree_Node_Queue = (PROJECT_TREE_NODE **)VB_malloc(MAXQUEUE*sizeof(PROJECT_TREE_NODE *));

  if (Tree_Node_Queue == NULL)
  {
    Fatal_MAError("priority queue for vista/light buffer");
  }
}



/*****************************************************************************
*
* FUNCTION      : Prune_Vista_Tree
*
* ARGUMENTS     : y - Current scanline number
*
* MODIFIED ARGS : none
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Prune vista tree, i.e. mark all nodes not on the current line inactive.
*
* CHANGES
*
*   -
*
******************************************************************************/

void Prune_Vista_Tree(y)
int y;
{
  unsigned size;
  unsigned short i;
  PROJECT_TREE_NODE *Node, *Sib;

  size = 0;

#ifdef COUNT_PROJECT_TESTS
  Project_Tests++;
#endif

  if ((y < Root_Vista->Project.y1) || (y > Root_Vista->Project.y2))
  {
    /* Root doesn't lie on current line --> prune root */

    Root_Vista->is_leaf |= PRUNE_TEMPORARY;
  }
  else
  {
    /* Root lies on current line --> unprune root */

#ifdef COUNT_PROJECT_TESTS
    Project_Tests_Succeeded++;
#endif

    Root_Vista->is_leaf &= PRUNE_TEMPORARY_INVERS;

    Tree_Node_Queue[size++] = Root_Vista;
  }

  while (size > 0)
  {
    Node = Tree_Node_Queue[--size];

    if (Node->is_leaf & TRUE)
    {
#ifdef COUNT_PROJECT_TESTS
      Project_Tests++;
#endif

      if ((y < Node->Project.y1) || (y > Node->Project.y2))
      {
	/* Leaf doesn't lie on current line --> prune leaf */

	Node->is_leaf |= PRUNE_TEMPORARY;
      }
      else
      {
	/* Leaf lies on current line --> unprune leaf */

#ifdef COUNT_PROJECT_TESTS
	Project_Tests_Succeeded++;
#endif

	Node->is_leaf &= PRUNE_TEMPORARY_INVERS;
      }
    }
    else
    {
      /* Check siblings of the node */

      for (i = 0; i < Node->Entries; i++)
      {
	Sib = Node->Entry[i];

#ifdef COUNT_PROJECT_TESTS
	Project_Tests++;
#endif

	if ((y < Sib->Project.y1) || (y > Sib->Project.y2))
	{
	  /* Sibling doesn't lie on current line --> prune sibling */

	  Sib->is_leaf |= PRUNE_TEMPORARY;
	}
	else
	{
	  /* Sibling lies on current line --> unprune sibling */

#ifdef COUNT_PROJECT_TESTS
	  Project_Tests_Succeeded++;
#endif

	  Sib->is_leaf &= PRUNE_TEMPORARY_INVERS;

	  /* Add sibling to list */

	  Tree_Node_Queue[size++] = Sib;
	}
      }
    }
  }
}



/*****************************************************************************
*
* FUNCTION      : Trace_Primary_Ray
*
* ARGUMENTS     : Ray    - Current ray
*                 Colour - Ray's colour
*                 x      - Current x-coordinate
*
* MODIFIED ARGS : colour
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Trace a primary ray using the vista tree.
*
* CHANGES
*
*   -
*
******************************************************************************/

void Trace_Primary_Ray (Ray, Colour, x)
RAY *Ray;
COLOUR *Colour;
int x;
{
  INTERSECTION Best_Intersection;

  COOPERATE

  Number_Of_Rays++;

  if (Trace_Level > (int)Max_Trace_Level)
    return;

  Colour->Red = Colour->Green = Colour->Blue = 0.0;

  Best_Intersection.Depth = BOUND_HUGE;

  /* What objects does this ray intersect? */

  if (Intersect_Vista_Tree(Ray, Root_Vista, x, &Best_Intersection))
  {
    Determine_Apparent_Colour (&Best_Intersection, Colour, Ray);
  }
  else
  {
    if (Frame.Fog_Distance > 0.0)
    {
      *Colour = Frame.Fog_Colour;
    }
    else
    {
      *Colour = Frame.Background_Colour;
    }
  }
}



/*****************************************************************************
*
* FUNCTION      : Create_Rayinfo
*
* ARGUMENTS     : Ray     - Current ray
*                 rayinfo - Rayinfo structure
*
* MODIFIED ARGS : rayinfo
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Calculate the rayinfo structure for a given ray. It's need for
*   intersection the ray with bounding boxes.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Create_Rayinfo(Ray, rayinfo)
RAY *Ray;
RAYINFO *rayinfo;
{
  DBL t;

  /* Create the direction vectors for this ray */

  rayinfo->slab_num = Ray->Initial;

  if ((rayinfo->nonzero.x = ((t = Ray->Direction.x) != 0.0)) != 0)
  {
    rayinfo->slab_den.x = 1.0 / t;
    rayinfo->positive.x = (Ray->Direction.x > 0.0);
  }

  if ((rayinfo->nonzero.y = ((t = Ray->Direction.y) != 0.0)) != 0)
  {
    rayinfo->slab_den.y = 1.0 / t;
    rayinfo->positive.y = (Ray->Direction.y > 0.0);
  }

  if ((rayinfo->nonzero.z = ((t = Ray->Direction.z) != 0.0)) != 0)
  {
    rayinfo->slab_den.z = 1.0 / t;
    rayinfo->positive.z = (Ray->Direction.z > 0.0);
  }
}



/*****************************************************************************
*
* FUNCTION      : Intersect_Vista_Tree
*
* ARGUMENTS     : Ray               - Primary ray
*                 Tree              - Vista tree's top-node
*                 x                 - Current x-coordinate
*                 Best_Intersection - Intersection found
*
* MODIFIED ARGS : Best_Intersection
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Intersect a PRIMARY ray with the vista tree
*   (tree pruning is used can be primary ray!!!).
*
* CHANGES
*
*   -
*
******************************************************************************/

static int Intersect_Vista_Tree(Ray, Tree, x, Best_Intersection)
RAY *Ray;
PROJECT_TREE_NODE *Tree;
int x;
INTERSECTION *Best_Intersection;
{
  INTERSECTION New_Intersection;
  unsigned short i;
  unsigned Qsize, size;
  int Found;
  RAYINFO rayinfo;
  DBL key;
  OBJECT *Object;
  PROJECT_TREE_NODE *Node;

  /* Start with an empty priority queue */

  Qsize = 0;

  Found = FALSE;

  totalQueueResets++;

  /* Which algorithm to use? */

  switch (Project_Algorithm)
  {
    /************************************************************************
    * Do not use any bounding slab tests. just use their projections.
    *************************************************************************/

    case PROJECTIONS_WITHOUT_SLABS :

      /* Check root */

#ifdef COUNT_PROJECT_TESTS
      Project_Tests++;
#endif

      if ((x >= Tree->Project.x1) && (x <= Tree->Project.x2))
      {
#ifdef COUNT_PROJECT_TESTS
	Project_Tests_Succeeded++;
#endif

	Tree_Node_Queue[Qsize++] = Tree;
      }

      while (Qsize > 0)
      {
	Tree = Tree_Node_Queue[--Qsize];

	switch (Tree->is_leaf)
	{
	  case FALSE:

	    /* Check siblings of the unpruned node in 2d */

	    for (i = 0; i < Tree->Entries; i++)
	    {
	      Node = Tree->Entry[i];

	      /* Check unpruned siblings only */

	      if (Node->is_leaf < PRUNE_CHECK)
	      {
#ifdef COUNT_PROJECT_TESTS
		Project_Tests++;
#endif

		if ((x >= Node->Project.x1) && (x <= Node->Project.x2))
		{
#ifdef COUNT_PROJECT_TESTS
		  Project_Tests_Succeeded++;
#endif

		  /* Add node to node queue */

		  if (Qsize >= MAXQUEUE)
		    Fatal_Error("Node queue overflow.\n");

		  Tree_Node_Queue[Qsize++] = Node;
		}
	      }
	    }

	    break;

	  case TRUE:

	    /* Unpruned leaf --> test object */

	    if (Intersection(&New_Intersection, ((PROJECT_TREE_LEAF *)Tree)->Object, Ray))
	    {
	      if (New_Intersection.Depth < Best_Intersection->Depth)
	      {
		*Best_Intersection = New_Intersection;
		Found = TRUE;
	      }
	    }

	    break;

       /* default:

	    The node/leaf is pruned and needn't be checked */

	}
      }

      break;



    /************************************************************************
    * Use bounding slab tests for leaves (i.e. objects) in the tree.
    *************************************************************************/

    case PROJECTIONS_WITH_LEAF_SLABS :

      /* Create the direction vectors for this ray */

      Create_Rayinfo(Ray, &rayinfo);

      /* Fill the priority queue with all possible candidates */

      size = 0;

      /* Check root */

#ifdef COUNT_PROJECT_TESTS
      Project_Tests++;
#endif

      if ((x >= Tree->Project.x1) && (x <= Tree->Project.x2))
      {
#ifdef COUNT_PROJECT_TESTS
	Project_Tests_Succeeded++;
#endif

	Tree_Node_Queue[size++] = Tree;
      }

      while (size > 0)
      {
	Tree = Tree_Node_Queue[--size];

	switch (Tree->is_leaf)
	{
	  case FALSE:

	    /* Check siblings of the unpruned node in 2d */

	    for (i = 0; i < Tree->Entries; i++)
	    {
	      Node = Tree->Entry[i];

	      /* Check unpruned siblings only */

	      if (Node->is_leaf < PRUNE_CHECK)
	      {
#ifdef COUNT_PROJECT_TESTS
		Project_Tests++;
#endif

		if ((x >= Node->Project.x1) && (x <= Node->Project.x2))
		{
#ifdef COUNT_PROJECT_TESTS
		  Project_Tests_Succeeded++;
#endif

		  /* add node to node queue */

		  if (size >= MAXQUEUE)
		    Fatal_Error("Node queue overflow.\n");

		  Tree_Node_Queue[size++] = Node;
		}
	      }
	    }

	break;

	  case TRUE:

	    /* Unpruned leaf --> test object's bounding box in 3d */

	    CheckAndEnqueue(Tree_Leaf_Queue, &Qsize, ((PROJECT_TREE_LEAF *)Tree)->Object,
	      &((PROJECT_TREE_LEAF *)Tree)->Object->Bounds, &rayinfo);

	    break;

       /* default:

	    The node/leaf is pruned and needn't be checked */

	}
      }

      /* Now test the candidates in the priority queue */

      while (Qsize > 0)
      {
	PriorityQueueDelete(Tree_Leaf_Queue, &Qsize, &key, &Object);

	if (key > Best_Intersection->Depth)
	  break;

	if (Intersection(&New_Intersection, Object, Ray))
	{
	  if (New_Intersection.Depth < Best_Intersection->Depth)
	  {
	    *Best_Intersection = New_Intersection;
	    Found = TRUE;
	  }
	}
      }

      break;



    /************************************************************************
    * Use bounding slab tests for all nodes that pass the projection test.
    *************************************************************************/

    case PROJECTIONS_WITH_NODE_SLABS :

      /* Create the direction vectors for this ray */

      Create_Rayinfo(Ray, &rayinfo);

      /* Check root */

#ifdef COUNT_PROJECT_TESTS
      Project_Tests++;
#endif

      if ((x >= Tree->Project.x1) && (x <= Tree->Project.x2))
      {
#ifdef COUNT_PROJECT_TESTS
	Project_Tests_Succeeded++;
#endif

	CheckAndEnqueue(Tree_Leaf_Queue, &Qsize, (OBJECT *)Tree,
	  &Tree->Object->Bounds, &rayinfo);
      }

      while (Qsize > 0)
      {
	PriorityQueueDelete(Tree_Leaf_Queue, &Qsize, &key, &Object);

	if (key > Best_Intersection->Depth)
	  break;

	Tree = (PROJECT_TREE_NODE *)Object;

	switch (Tree->is_leaf)
	{
	  case FALSE:

	    /* Check siblings of the unpruned node in 2d */

	    for (i = 0; i < Tree->Entries; i++)
	    {
	      Node = Tree->Entry[i];

	      /* Check unpruned siblings only */

	      if (Node->is_leaf < PRUNE_CHECK)
	      {
#ifdef COUNT_PROJECT_TESTS
		Project_Tests++;
#endif

		if ((x >= Node->Project.x1) && (x <= Node->Project.x2))
		{
#ifdef COUNT_PROJECT_TESTS
		  Project_Tests_Succeeded++;
#endif

		  CheckAndEnqueue(Tree_Leaf_Queue, &Qsize, (OBJECT *)Node,
		    &Node->Object->Bounds, &rayinfo);
		}
	      }
	    }

	    break;

	  case TRUE:

	    /* Unpruned leaf --> test object */

	    if (Intersection(&New_Intersection, ((PROJECT_TREE_LEAF *)Tree)->Object, Ray))
	    {
	      if (New_Intersection.Depth < Best_Intersection->Depth)
	      {
		*Best_Intersection = New_Intersection;
		Found = TRUE;
	      }
	    }

	    break;

       /* default:

	    The node/leaf is pruned and needn't be checked */

	}
      }

      break;



    default :
      Fatal_Error("Unkown algorithm.\n");
  }

  return(Found);
}



/*****************************************************************************
*
* FUNCTION      : Intersect_Light_Tree
*
* ARGUMENTS     : Ray               - Shadow ray
*                 Tree              - Light tree's top-node
*                 x                 - X-coordinate of the shadow ray
*                 y                 - Y-coordinate of the shadow ray
*                 Best_Intersection - Intersection found
*                 Best_Object       - Object found
*
* MODIFIED ARGS : Best_Intersection, Best_Object
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Intersect a shadow ray with the light tree
*   (same as for the vista tree but without pruning).
*
* CHANGES
*
*   -
*
******************************************************************************/

int Intersect_Light_Tree(Ray, Tree, x, y, Best_Intersection, Best_Object)
RAY *Ray;
PROJECT_TREE_NODE *Tree;
int x, y;
INTERSECTION *Best_Intersection;
OBJECT **Best_Object;
{
  INTERSECTION New_Intersection;
  unsigned short i;
  unsigned Qsize, size;
  int Found;
  RAYINFO rayinfo;
  DBL key;
  OBJECT *Object;
  PROJECT_TREE_NODE *Node;

  /* Start with an empty priority queue */

  Qsize = 0;

  Found = FALSE;

  totalQueueResets++;

  /* which algorithm to use? */

  switch (Project_Algorithm)
  {
    /************************************************************************
    * Do not use any bounding slab tests. just use their projections.
    *************************************************************************/

    case PROJECTIONS_WITHOUT_SLABS :

      /* Check root */

#ifdef COUNT_PROJECT_TESTS
      Project_Tests++;
#endif

      if ((x >= Tree->Project.x1) && (x <= Tree->Project.x2) &&
	  (y >= Tree->Project.y1) && (y <= Tree->Project.y2))
      {
#ifdef COUNT_PROJECT_TESTS
	Project_Tests_Succeeded++;
#endif

	Tree_Node_Queue[Qsize++] = Tree;
      }

      while (Qsize > 0)
      {
	Tree = Tree_Node_Queue[--Qsize];

	if (Tree->is_leaf)
	{
	  /* Leaf --> test object */

	  if (Intersection(&New_Intersection, ((PROJECT_TREE_LEAF *)Tree)->Object, Ray))
	  {
	    if (New_Intersection.Depth < Best_Intersection->Depth)
	    {
	      *Best_Intersection = New_Intersection;
	      *Best_Object = ((PROJECT_TREE_LEAF *)Tree)->Object;
	      Found = TRUE;
	    }
	  }
	}
	else
	{
	  /* Check siblings of the node in 2d */

	  for (i = 0; i < Tree->Entries; i++)
	  {
	    Node = Tree->Entry[i];

#ifdef COUNT_PROJECT_TESTS
	    Project_Tests++;
#endif

	    if ((x >= Node->Project.x1) && (x <= Node->Project.x2) &&
		(y >= Node->Project.y1) && (y <= Node->Project.y2))
	    {
#ifdef COUNT_PROJECT_TESTS
	      Project_Tests_Succeeded++;
#endif

	      /* Add node to node queue */

	      if (Qsize >= MAXQUEUE)
		Fatal_Error("Node queue overflow.\n");

	      Tree_Node_Queue[Qsize++] = Node;
	    }
	  }
	}
      }

      break;



    /************************************************************************
    * Use bounding slab tests for leaves (i.e. objects) in the tree.
    *************************************************************************/

    case PROJECTIONS_WITH_LEAF_SLABS :

      /* Create the direction vectors for this ray */

      Create_Rayinfo(Ray, &rayinfo);

      /* Fill the priority queue with all possible candidates */

      size = 0;

      /* Check root */

#ifdef COUNT_PROJECT_TESTS
      Project_Tests++;
#endif

      if ((x >= Tree->Project.x1) && (x <= Tree->Project.x2) &&
	  (y >= Tree->Project.y1) && (y <= Tree->Project.y2))
      {
#ifdef COUNT_PROJECT_TESTS
	Project_Tests_Succeeded++;
#endif

	Tree_Node_Queue[size++] = Tree;
      }

      while (size > 0)
      {
	Tree = Tree_Node_Queue[--size];

	if (Tree->is_leaf)
	{
	  /* Leaf --> test object's bounding box in 3d */

	  CheckAndEnqueue(Tree_Leaf_Queue, &Qsize, ((PROJECT_TREE_LEAF *)Tree)->Object,
	    &((PROJECT_TREE_LEAF *)Tree)->Object->Bounds, &rayinfo);
	}
	else
	{
	  /* Check siblings of the node in 2d */

	  for (i = 0; i < Tree->Entries; i++)
	  {
	    Node = Tree->Entry[i];

#ifdef COUNT_PROJECT_TESTS
	    Project_Tests++;
#endif

	    if ((x >= Node->Project.x1) && (x <= Node->Project.x2) &&
		(y >= Node->Project.y1) && (y <= Node->Project.y2))
	    {
#ifdef COUNT_PROJECT_TESTS
	      Project_Tests_Succeeded++;
#endif

	      /* Add node to node queue */

	      if (size >= MAXQUEUE)
		Fatal_Error("Node queue overflow.\n");

	      Tree_Node_Queue[size++] = Node;
	    }
	  }
	}
      }

      /* Now test the candidates in the priority queue */

      while (Qsize > 0)
      {
	PriorityQueueDelete(Tree_Leaf_Queue, &Qsize, &key, &Object);

	if (key > Best_Intersection->Depth)
	  break;

	if (Intersection(&New_Intersection, Object, Ray))
	{
	  if (New_Intersection.Depth < Best_Intersection->Depth)
	  {
	    *Best_Intersection = New_Intersection;
	    *Best_Object = Object;
	    Found = TRUE;
	  }
	}
      }

      break;



    /************************************************************************
    * Use bounding slab tests for all nodes that pass the projection test.
    *************************************************************************/

    case PROJECTIONS_WITH_NODE_SLABS :

      /* Create the direction vectors for this ray */

      Create_Rayinfo(Ray, &rayinfo);

      /* Check root */

#ifdef COUNT_PROJECT_TESTS
      Project_Tests++;
#endif

      if ((x >= Tree->Project.x1) && (x <= Tree->Project.x2) &&
	  (y >= Tree->Project.y1) && (y <= Tree->Project.y2))
      {
#ifdef COUNT_PROJECT_TESTS
	Project_Tests_Succeeded++;
#endif

	CheckAndEnqueue(Tree_Leaf_Queue, &Qsize, (OBJECT *)Tree,
	  &Tree->Object->Bounds, &rayinfo);
      }

      while (Qsize > 0)
      {
	PriorityQueueDelete(Tree_Leaf_Queue, &Qsize, &key, &Object);

	if (key > Best_Intersection->Depth)
	  break;

	Tree = (PROJECT_TREE_NODE *)Object;

	if (Tree->is_leaf)
	{
	  /* Leaf --> test object */

	  if (Intersection(&New_Intersection, ((PROJECT_TREE_LEAF *)Tree)->Object, Ray))
	  {
	    if (New_Intersection.Depth < Best_Intersection->Depth)
	    {
	      *Best_Intersection = New_Intersection;
	      *Best_Object = ((PROJECT_TREE_LEAF *)Tree)->Object;
	      Found = TRUE;
	    }
	  }
	}
	else
	{
	  /* Check siblings of the unpruned node in 2d */

	  for (i = 0; i < Tree->Entries; i++)
	  {
	    Node = Tree->Entry[i];

#ifdef COUNT_PROJECT_TESTS
	    Project_Tests++;
#endif

	    if ((x >= Node->Project.x1) && (x <= Node->Project.x2) &&
		(y >= Node->Project.y1) && (y <= Node->Project.y2))
	    {
#ifdef COUNT_PROJECT_TESTS
	      Project_Tests_Succeeded++;
#endif

	      CheckAndEnqueue(Tree_Leaf_Queue, &Qsize, (OBJECT *)Node,
		&Node->Object->Bounds, &rayinfo);
	    }
	  }
	}
      }

      break;



    default :
      Fatal_Error("Unkown algorithm.\n");
  }

  return(Found);
}



/*****************************************************************************
*
* FUNCTION      : Project_triangle
*
* ARGUMENTS     : Project    - Triangle's projection
*                 P1, P2, P3 - Triangle's edges
*                 visible    - Flag if triangle is visible
*
* MODIFIED ARGS : Project, visible
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Project a triangle onto the screen.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Project_triangle (Project, P1, P2, P3, visible)
PROJECT *Project;
VECTOR P1, P2, P3;
int *visible;
{
  VECTOR Points[MAX_CLIP_POINTS];
  int i, number;
  int x, y;

  Points[0] = P1;
  Points[1] = P2;
  Points[2] = P3;

  number = 3;

  /* Clip triangle only if some quick tests say it's necessary.
     Assuming that only a few triangles need clipping this saves some time.
     (I don't need to write fabs(1+P?.z) since the tests succeed anyway if
      P?.z < -1. Hope the compiler doesn't change the tests' order!) */

  if ((P1.z < -1.0) || (P2.z < -1.0) || (P3.z < -1.0) ||
      (fabs(P1.x) > 0.5*(1.0+P1.z)) || (fabs(P1.y) > 0.5*(1.0+P1.z)) ||
      (fabs(P2.x) > 0.5*(1.0+P2.z)) || (fabs(P2.y) > 0.5*(1.0+P2.z)) ||
      (fabs(P3.x) > 0.5*(1.0+P3.z)) || (fabs(P3.y) > 0.5*(1.0+P3.z)))
  {
    Clip_Polygon(Points, &number, &VIEW_VX1, &VIEW_VX2, &VIEW_VY1, &VIEW_VY2,
				   VIEW_DX1,  VIEW_DX2,  VIEW_DY1,  VIEW_DY2);
  }

  if (!number)
    return;

  for (i = 0; i < number; i++)
  {
    if (Points[i].z == -1.0)
    {
      Points[i].x = Points[i].y = 0.0;
    }
    else
    {
      Points[i].x = Points[i].x / (1.0 + Points[i].z);
      Points[i].y = Points[i].y / (1.0 + Points[i].z);
    }

    x = Frame.Screen_Width/2  + (int)(Frame.Screen_Width  * Points[i].x);
    y = Frame.Screen_Height/2 - (int)(Frame.Screen_Height * Points[i].y);

    if (x < Project->x1) Project->x1 = x;
    if (x > Project->x2) Project->x2 = x;
    if (y < Project->y1) Project->y1 = y;
    if (y > Project->y2) Project->y2 = y;
  }

  *visible = TRUE;
}



/*****************************************************************************
*
* FUNCTION      : Project_rectangle
*
* ARGUMENTS     : Project        - Rectangle's projection
*                 P1, P2, P3, P4 - Rectangle's edges
*                 visible        - Flag if rectangle is visible
*
* MODIFIED ARGS : Project, visible
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Project a rectangle onto the screen.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Project_rectangle(Project, P1, P2, P3, P4, visible)
PROJECT *Project;
VECTOR P1, P2, P3, P4;
int *visible;
{
  VECTOR Points[MAX_CLIP_POINTS];
  int i, number;
  int x, y;

  Points[0] = P1;
  Points[1] = P2;
  Points[2] = P3;
  Points[3] = P4;

  number = 4;

  /* Clip square only if some quick tests say it's necessary.
     Assuming that only a few squares need clipping this saves some time.
     (I don't need to write fabs(1+P?.z) since the tests succeed anyway if
      P?.z < -1. Hope the compiler doesn't change the tests' order!) */

  if ((P1.z < -1.0) || (P2.z < -1.0) || (P3.z < -1.0) || (P4.z < -1.0) ||
      (fabs(P1.x) > 0.5*(1.0+P1.z)) || (fabs(P1.y) > 0.5*(1.0+P1.z)) ||
      (fabs(P2.x) > 0.5*(1.0+P2.z)) || (fabs(P2.y) > 0.5*(1.0+P2.z)) ||
      (fabs(P3.x) > 0.5*(1.0+P3.z)) || (fabs(P3.y) > 0.5*(1.0+P3.z)) ||
      (fabs(P4.x) > 0.5*(1.0+P4.z)) || (fabs(P4.y) > 0.5*(1.0+P4.z)))
  {
    Clip_Polygon(Points, &number, &VIEW_VX1, &VIEW_VX2, &VIEW_VY1, &VIEW_VY2,
				   VIEW_DX1,  VIEW_DX2,  VIEW_DY1,  VIEW_DY2);
  }

  if (!number)
    return;

  for (i = 0; i < number; i++)
  {
    if (Points[i].z == -1.0)
    {
      Points[i].x = Points[i].y = 0.0;
    }
    else
    {
      Points[i].x = Points[i].x / (1.0 + Points[i].z);
      Points[i].y = Points[i].y / (1.0 + Points[i].z);
    }

    x = Frame.Screen_Width/2  + (int)(Frame.Screen_Width  * Points[i].x);
    y = Frame.Screen_Height/2 - (int)(Frame.Screen_Height * Points[i].y);

    if (x < Project->x1) Project->x1 = x;
    if (x > Project->x2) Project->x2 = x;
    if (y < Project->y1) Project->y1 = y;
    if (y > Project->y2) Project->y2 = y;
  }

  *visible = TRUE;
}




/*****************************************************************************
*
* FUNCTION      : Project_BBox
*
* ARGUMENTS     : Project - Box's projection
*                 Points  - Box's edges
*                 visible - Flag if box is visible
*
* MODIFIED ARGS : Project, visible
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Project a box onto the screen.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Project_Bbox(Project, Points, visible)
PROJECT *Project;
VECTOR *Points;
int *visible;
{
  int vis;
  PROJECT New;

  New.x1 = MAX_VB_ENTRY;
  New.x2 = MIN_VB_ENTRY;
  New.y1 = MAX_VB_ENTRY;
  New.y2 = MIN_VB_ENTRY;

  vis = FALSE;

  Project_rectangle(&New, Points[0], Points[1], Points[3], Points[2], &vis);
  Project_rectangle(&New, Points[4], Points[5], Points[7], Points[6], &vis);
  Project_rectangle(&New, Points[0], Points[1], Points[5], Points[4], &vis);
  Project_rectangle(&New, Points[2], Points[3], Points[7], Points[6], &vis);
  Project_rectangle(&New, Points[1], Points[3], Points[7], Points[5], &vis);
  Project_rectangle(&New, Points[0], Points[2], Points[6], Points[4], &vis);

  if (vis)
  {
    if (New.x1 > Project->x1) Project->x1 = New.x1;
    if (New.x2 < Project->x2) Project->x2 = New.x2;
    if (New.y1 > Project->y1) Project->y1 = New.y1;
    if (New.y2 < Project->y2) Project->y2 = New.y2;
    *visible = TRUE;
  }
}



/*****************************************************************************
*
* FUNCTION      : Project_Bounds
*
* ARGUMENTS     : Project - Bounding box's projection
*                 Bounds  - Bounding box
*                 visible - Flag if bounding box is visible
*
* MODIFIED ARGS : Project, visible
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Project a bounding box onto the screen.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Project_Bounds(Project, Bounds, visible)
PROJECT *Project;
BBOX *Bounds;
int *visible;
{
  int i;
  VECTOR P[8];

  for (i = 0; i<8; i++)
  {
    P[i] = Bounds->Lower_Left;

    P[i].x += ((i & 1) ? Bounds->Lengths.x : 0.0);
    P[i].y += ((i & 2) ? Bounds->Lengths.y : 0.0);
    P[i].z += ((i & 4) ? Bounds->Lengths.z : 0.0);

    Transform_Point (&P[i]);
  }

  Project_Bbox(Project, &P[0], visible);
}



/*****************************************************************************
*
* FUNCTION      : Project_Bicubic_Patch
*
* ARGUMENTS     : Project - Projection
*                 Object  - Object
*
* MODIFIED ARGS : Project
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Project the bounding sphere of a bicubic patch onto the screen.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Project_Bicubic_Patch(Project, Object)
PROJECT *Project;
OBJECT *Object;
{
  BICUBIC_PATCH *patch;
  DBL x, y, z, r, a20, a02, a11, a10, a01, a00;
  MATRIX A;

  patch = (BICUBIC_PATCH *)Object;

  /* Set up 4x4 quadric matrix A */

  x = patch->Bounding_Sphere_Center.x;
  y = patch->Bounding_Sphere_Center.y;
  z = patch->Bounding_Sphere_Center.z;

  r = 1.0 / (patch->Bounding_Sphere_Radius * patch->Bounding_Sphere_Radius);

  A[0][0] = r;    A[0][1] = 0.0;  A[0][2] = 0.0;  A[0][3] = -x*r;
  A[1][0] = 0.0;  A[1][1] = r;    A[1][2] = 0.0;  A[1][3] = -y*r;
  A[2][0] = 0.0;  A[2][1] = 0.0;  A[2][2] = r;    A[2][3] = -z*r;
  A[3][0] = -x*r; A[3][1] = -y*r; A[3][2] = -z*r; A[3][3] = r*(x*x+y*y+z*z) - 1.0;

  Project_Vista(&A, NULL, &a20, &a02, &a11, &a10, &a01, &a00);

  /* Get vista's bounding rectangle */

  Get_Ellipse_Projection(Project, a20, a02, a11, a10, a01, a00);
}



/*****************************************************************************
*
* FUNCTION      : Project_Box
*
* ARGUMENTS     : Project - Projection
*                 Object  - Object
*                 visible - Flag if object is visible
*
* MODIFIED ARGS : Project, visible
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Project a box onto the screen.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Project_Box(Project, Object, visible)
PROJECT *Project;
OBJECT *Object;
int *visible;
{
  int i;
  VECTOR P[8];
  BOX *box;

  box = (BOX *)Object;

  for (i = 0; i<8; i++)
  {
    P[i] = box->bounds[0];

    if (i & 1) P[i].x = box->bounds[1].x;
    if (i & 2) P[i].y = box->bounds[1].y;
    if (i & 4) P[i].z = box->bounds[1].z;

    if (box->Trans != NULL)
      MTransPoint(&P[i], &P[i], box->Trans);

    Transform_Point (&P[i]);
  }

  Project_Bbox(Project, &P[0], visible);
}



/*****************************************************************************
*
* FUNCTION      : Project_Height_Field
*
* ARGUMENTS     : Project - Projection
*                 Object  - Object
*                 visible - Flag if object is visible
*
* MODIFIED ARGS : Project, visible
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Project the bounding box of a height field onto the screen.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Project_Height_Field(Project, Object, visible)
PROJECT *Project;
OBJECT *Object;
int *visible;
{
  int i;
  VECTOR P[8];
  HEIGHT_FIELD *hfield;

  hfield = (HEIGHT_FIELD *)Object;

  for (i = 0; i<8; i++)
  {
    P[i] = hfield->bounding_box->bounds[0];

    if (i & 1) P[i].x = hfield->bounding_box->bounds[1].x;
    if (i & 2) P[i].y = hfield->bounding_box->bounds[1].y;
    if (i & 4) P[i].z = hfield->bounding_box->bounds[1].z;

    if (hfield->Trans != NULL)
      MTransPoint(&P[i], &P[i], hfield->Trans);

    Transform_Point (&P[i]);
  }

  Project_Bbox(Project, &P[0], visible);
}



/*****************************************************************************
*
* FUNCTION      : Project_Triangle
*
* ARGUMENTS     : Project - Projection
*                 Object  - Object
*                 visible - Flag if object is visible
*
* MODIFIED ARGS : Project, visible
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Project a triangle onto the screen.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Project_Triangle(Project, Object, visible)
PROJECT *Project;
OBJECT *Object;
int *visible;
{
  int i, vis;
  VECTOR P[3];
  PROJECT New;

  New.x1 = MAX_VB_ENTRY;
  New.x2 = MIN_VB_ENTRY;
  New.y1 = MAX_VB_ENTRY;
  New.y2 = MIN_VB_ENTRY;

  P[0] = ((TRIANGLE *)Object)->P1;
  P[1] = ((TRIANGLE *)Object)->P2;
  P[2] = ((TRIANGLE *)Object)->P3;

  for (i = 0; i < 3; i++)
  {
    Transform_Point(&P[i]);
  }

  vis = FALSE;

  Project_triangle(&New, P[0], P[1], P[2], &vis);

  if (vis)
  {
    if (New.x1 > Project->x1) Project->x1 = New.x1;
    if (New.x2 < Project->x2) Project->x2 = New.x2;
    if (New.y1 > Project->y1) Project->y1 = New.y1;
    if (New.y2 < Project->y2) Project->y2 = New.y2;
    *visible = TRUE;
  }
}



/*****************************************************************************
*
* FUNCTION      : Project_Smooth_Triangle
*
* ARGUMENTS     : Project - Projection
*                 Object  - Object
*                 visible - Flag if object is visible
*
* MODIFIED ARGS : Project, visible
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Project a smooth triangle onto the screen.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Project_Smooth_Triangle(Project, Object, visible)
PROJECT *Project;
OBJECT *Object;
int *visible;
{
  int i, vis;
  VECTOR P[3];
  PROJECT New;

  New.x1 = MAX_VB_ENTRY;
  New.x2 = MIN_VB_ENTRY;
  New.y1 = MAX_VB_ENTRY;
  New.y2 = MIN_VB_ENTRY;

  P[0] = ((SMOOTH_TRIANGLE *)Object)->P1;
  P[1] = ((SMOOTH_TRIANGLE *)Object)->P2;
  P[2] = ((SMOOTH_TRIANGLE *)Object)->P3;

  for (i = 0; i < 3; i++)
  {
    Transform_Point(&P[i]);
  }

  vis = FALSE;

  Project_triangle(&New, P[0], P[1], P[2], &vis);

  if (vis)
  {
    if (New.x1 > Project->x1) Project->x1 = New.x1;
    if (New.x2 < Project->x2) Project->x2 = New.x2;
    if (New.y1 > Project->y1) Project->y1 = New.y1;
    if (New.y2 < Project->y2) Project->y2 = New.y2;
    *visible = TRUE;
  }
}



/*****************************************************************************
*
* FUNCTION      : Project_Sphere
*
* ARGUMENTS     : Project - Projection
*                 Object  - Oject
*
* MODIFIED ARGS : Project
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Project a sphere onto the screen.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Project_Sphere(Project, Object)
PROJECT *Project;
OBJECT *Object;
{
  SPHERE *sphere;
  DBL x, y, z, r, a20, a02, a11, a10, a01, a00;
  MATRIX A;

  sphere = (SPHERE *)Object;

  /* Set up 4x4 quadric matrix A0 */

  x = sphere->Center.x;
  y = sphere->Center.y;
  z = sphere->Center.z;

  r = 1.0/sphere->Radius_Squared;

  A[0][0] = r;    A[0][1] = 0.0;  A[0][2] = 0.0;  A[0][3] = -x*r;
  A[1][0] = 0.0;  A[1][1] = r;    A[1][2] = 0.0;  A[1][3] = -y*r;
  A[2][0] = 0.0;  A[2][1] = 0.0;  A[2][2] = r;    A[2][3] = -z*r;
  A[3][0] = -x*r; A[3][1] = -y*r; A[3][2] = -z*r; A[3][3] = r*(x*x+y*y+z*z) - 1.0;

  Project_Vista(&A, sphere->Trans, &a20, &a02, &a11, &a10, &a01, &a00);

  /* get vista's bounding rectangle */

  Get_Ellipse_Projection(Project, a20, a02, a11, a10, a01, a00);
}



/*****************************************************************************
*
* FUNCTION      : Get_Ellipse_Projection
*
* ARGUMENTS     : Project        - Projection
*                 a20, a02, a11,
*                 a10, a01, a00  - Parameters of the quadratic curve
*
* MODIFIED ARGS : Project
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Get the bounding rectangle around an ellipse, i.e. the
*   quadratic curve: a20*x*x + a02*y*y + a11*x*y + a10*x + a01*y + a00 = 0
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Get_Ellipse_Projection(Project, a20, a02, a11, a10, a01, a00)
PROJECT *Project;
DBL a20, a02, a11, a10, a01, a00;
{
  int x1i, x2i, y1i, y2i;
  DBL x1, y1, x2, y2;
  DBL a, b, c, d, k, k1, k2, k3, k4;

  x1 = y1 = -0.5;
  x2 = y2 = +0.5;

  k1 = a11/(-2.0*a20);
  k2 = a10/(-2.0*a20);

  k3 = a11/(-2.0*a02);
  k4 = a01/(-2.0*a02);

  a = a20 + a02*k3*k3 + a11*k3;

  b = a10 + 2.0*a02*k3*k4 + a01*k3 + a11*k4;

  c = a02*k4*k4 + a01*k4 + a00;

  d = b*b - 4.0*a*c;

  if (d > 0.0)
  {
    a = 2.0*a;

    d = sqrt(d);

    x1 = (-b+d)/a;
    x2 = (-b-d)/a;

    if (x1>x2)
    {
      k = x1; x1 = x2; x2 = k;
    }

    a = a02 + a20*k1*k1 + a11*k1;

    b = a01 + 2.0*a20*k1*k2 + a10*k1 + a11*k2;

    c = a20*k2*k2 + a10*k2 + a00;

    d = b*b - 4.0*a*c;

    if (d > 0.0)
    {
      a = 2.0*a;

      d = sqrt(d);

      y1 = (-b+d)/a;
      y2 = (-b-d)/a;

      if (y1>y2)
      {
	k = y1; y1 = y2; y2 = k;
      }
    }
  }

  if (x1 < -0.5) x1 = -0.5;
  if (x2 > +0.5) x2 = +0.5;
  if (y1 < -0.5) y1 = -0.5;
  if (y2 > +0.5) y2 = +0.5;

  x1i = Frame.Screen_Width/2 + (int)(Frame.Screen_Width * x1) - 1;
  x2i = Frame.Screen_Width/2 + (int)(Frame.Screen_Width * x2) + 1;

  y1i = Frame.Screen_Height/2 - (int)(Frame.Screen_Height * y2) - 1;
  y2i = Frame.Screen_Height/2 - (int)(Frame.Screen_Height * y1) + 1;

  if (x1i > Project->x1) Project->x1 = x1i;
  if (x2i < Project->x2) Project->x2 = x2i;

  if (y1i > Project->y1) Project->y1 = y1i;
  if (y2i < Project->y2) Project->y2 = y2i;
}



/*****************************************************************************
*
* FUNCTION      : Project_Vista
*
* ARGUMENTS     : A0             - Matrix defining the quadric surface
*                 Trans          - Transformation matrices
*                 a20, a02, a11,
*                 a10, a01, a00  - Parameters of the projected quadratic curve
*
* MODIFIED ARGS : a20, a02, a11, a10, a01, a00
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Project a quadric surface onto the screen (perspective projection)
*   and get the paramters of the resulting quadratic curve.
*
* CHANGES
*
*   -
*
******************************************************************************/

void Project_Vista(A0, Trans, a20, a02, a11, a10, a01, a00)
MATRIX *A0;
TRANSFORM *Trans;
DBL *a20, *a02, *a11, *a10, *a01, *a00;
{
  int i, j;
  MATRIX Tinv, A1, A2, B, Transposed_Matrix;
  DBL k1, k2, k3;

  /* Apply object transformations */

  if (Trans == NULL)
  {
    for (i=0;i<4;i++)
    {
      for (j=0;j<4;j++)
      {
	A1[i][j] = (*A0)[i][j];
      }
    }
  }
  else
  {
    MTranspose(&Transposed_Matrix, &Trans->inverse);
    MTimes(&A1, &Trans->inverse, A0);
    MTimes(&A1, &A1, &Transposed_Matrix);
  }

  /* Transform into viewing system */

  MTranspose(&Transposed_Matrix, &WC2VCinv);
  MTimes(&A2, &Transposed_Matrix, &A1);
  MTimes(&A2, &A2, &WC2VCinv);

  /* Calculate quadrics vista */

  MIdentity(&Tinv);

  Tinv[3][2] = -1.0;

  MTranspose(&Transposed_Matrix, &Tinv);

  MTimes(&B, &Transposed_Matrix, &A2);
  MTimes(&B, &B, &Tinv);

  k1 = (B[0][2]+B[2][0])/(-2.0*B[2][2]);

  k2 = (B[1][2]+B[2][1])/(-2.0*B[2][2]);

  k3 = (B[2][3]+B[3][2])/(-2.0*B[2][2]);

  *a20 = B[0][0] + k1*(B[0][2]+B[2][0]) + k1*k1*B[2][2];

  *a02 = B[1][1] + k2*(B[1][2]+B[2][1]) + k2*k2*B[2][2];

  *a10 = B[0][3]+B[3][0] + k1*(B[2][3]+B[3][2]) + k3*(B[0][2]+B[2][0]) + 2.0*k1*k3*B[2][2];

  *a01 = B[1][3]+B[3][1] + k2*(B[2][3]+B[3][2]) + k3*(B[1][2]+B[2][1]) + 2.0*k2*k3*B[2][2];

  *a11 = B[0][1]+B[1][0] + k1*(B[1][2]+B[2][1]) + k2*(B[0][2]+B[2][0]) + 2.0*k1*k2*B[2][2];

  *a00 = B[3][3] + k3*(B[2][3]+B[3][2]) + k3*k3*B[2][2];
}



/*****************************************************************************
*
* FUNCTION      : Transform_Point
*
* ARGUMENTS     : P - Point to transform
*
* MODIFIED ARGS : P
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Transform a point from the world coordinate system to the viewer's
*   coordinate system.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Transform_Point(P)
VECTOR *P;
{
  DBL x,y,z;

  x = P->x - O.x;
  y = P->y - O.y;
  z = P->z - O.z;

  P->x = U.x * x + U.y * y + U.z * z;
  P->y = V.x * x + V.y * y + V.z * z;
  P->z = W.x * x + W.y * y + W.z * z;
}



/*****************************************************************************
*
* FUNCTION      : Init_View_Coordinates
*
* ARGUMENTS     : none
*
* MODIFIED ARGS : none
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Init the matrices and vectors used to transform a point from
*   the world coordinate system to the viewer's coordinate system.
*
* CHANGES
*
*   -
*
******************************************************************************/

void Init_View_Coordinates()
{
  DBL k, up_length, right_length;
  MATRIX A, B;

  U = Frame.Camera->Right;
  V = Frame.Camera->Up;
  W = Frame.Camera->Direction;

  VAdd (O, Frame.Camera->Location, Frame.Camera->Direction);

  VNormalize(U,U);
  VNormalize(V,V);
  VNormalize(W,W);

  VDot(k, U,V);
  if (fabs(k) > EPSILON)
    Fatal_Error("Camera: Right-vector not perpendicular to Up-vector.\n");

  VDot(k, U,W);
  if (fabs(k) > EPSILON)
    Fatal_Error("Camera: Right-vector not perpendicular to Direction-vector.\n");

  VDot(k, V,W);
  if (fabs(k) > EPSILON)
    Fatal_Error("Camera: Up-vector not perpendicular to Direction-vector.\n");

  VLength (distance, Frame.Camera->Direction);

  VLength (up_length, Frame.Camera->Up);
  VLength (right_length, Frame.Camera->Right);

  VScaleEq (U, 1.0/right_length);
  VScaleEq (V, 1.0/up_length);
  VScaleEq (W, 1.0/distance);

  A[0][0] = U.x; A[0][1] = U.y; A[0][2] = U.z; A[0][3] = 0.0;
  A[1][0] = V.x; A[1][1] = V.y; A[1][2] = V.z; A[1][3] = 0.0;
  A[2][0] = W.x; A[2][1] = W.y; A[2][2] = W.z; A[2][3] = 0.0;
  A[3][0] = 0.0; A[3][1] = 0.0; A[3][2] = 0.0; A[3][3] = 1.0;

  B[0][0] = 1.0; B[0][1] = 0.0; B[0][2] = 0.0; B[0][3] = -O.x;
  B[1][0] = 0.0; B[1][1] = 1.0; B[1][2] = 0.0; B[1][3] = -O.y;
  B[2][0] = 0.0; B[2][1] = 0.0; B[2][2] = 1.0; B[2][3] = -O.z;
  B[3][0] = 0.0; B[3][1] = 0.0; B[3][2] = 0.0; B[3][3] = 1.0;

  MTimes(&WC2VC, &A, &B);
  MInvers(&WC2VCinv, &WC2VC);
}



/*****************************************************************************
*
* FUNCTION      : Get_Perspective_Projection
*
* ARGUMENTS     : Object   - Object to project
*                 Project  - Projection
*                 infinite - Flag if object is infinite
*
* MODIFIED ARGS : Project
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Get the perspective projection of a single object, i.e.
*   the smallest rectangle enclosing the object's image on the screen.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Get_Perspective_Projection(Object, Project, infinite)
OBJECT *Object;
PROJECT *Project;
int infinite;
{
  int visible;
  METHODS *Methods;

  visible = FALSE;

  Methods = Object->Methods;

  /* If the object is infinite, there's no sense of projecting */

  if (!infinite)
  {
    if ((Methods == &Box_Methods) ||
	(Methods == &Smooth_Triangle_Methods) ||
	(Methods == &Triangle_Methods) ||
	(Methods == &Csg_Height_Field_Methods) ||
	(Methods == &Height_Field_Methods))
    {
      if (Methods == &Box_Methods)
	Project_Box(Project, Object, &visible);

      if (Methods == &Csg_Height_Field_Methods)
	Project_Height_Field(Project, Object, &visible);

      if (Methods == &Height_Field_Methods)
	Project_Height_Field(Project, Object, &visible);

      if (Methods == &Smooth_Triangle_Methods)
	Project_Smooth_Triangle(Project, Object, &visible);

      if (Methods == &Triangle_Methods)
	Project_Triangle(Project, Object, &visible);
    }
    else
    {
      Project_Bounds(Project, &Object->Bounds, &visible);

/*
      if (Object->Methods == &Bicubic_Patch_Methods)
	Project_Bicubic_Patch(Project, Object);

      if (Object->Methods == &Ellipsoid_Methods)
	Project_Sphere(Project, Object);

      if (Object->Methods == &Sphere_Methods)
	Project_Sphere(Project, Object);
*/
    }
  }

  /* We don't want invisible objects to become visible due to
     an oversized bouding object, so we project the bounding
     object only if the object itself is visible! */

  if (visible & (Object->Bound != NULL))
  {
    if (Object->Bound->Methods == &Box_Methods)
      Project_Box(Project, Object->Bound, &visible);

    if (Object->Bound->Methods == &Ellipsoid_Methods)
      Project_Sphere(Project, Object->Bound);

    if (Object->Bound->Methods == &Sphere_Methods)
      Project_Sphere(Project, Object->Bound);
  }

  if (visible)
  {
    if (Options & ANTIALIAS)
    {
      /* Increase the rectangle to make sure that nothing will be missed.
	 For anti-aliased images increase by a larger amount. */

      Project->x1 = max (0,                     Project->x1 - 2);
      Project->x2 = min (Frame.Screen_Width-1,  Project->x2 + 2);
      Project->y1 = max (-1,                    Project->y1 - 2);
      Project->y2 = min (Frame.Screen_Height-1, Project->y2 + 2);
    }
    else
    {
      /* Increase the rectangle to make sure that nothing will be missed. */

      Project->x1 = max (0,                     Project->x1 - 1);
      Project->x2 = min (Frame.Screen_Width-1,  Project->x2 + 1);
      Project->y1 = max (0,                     Project->y1 - 1);
      Project->y2 = min (Frame.Screen_Height-1, Project->y2 + 1);
    }
  }
  else
  {
    if (!infinite)
    {
      /* Object is invisible (the camera can't see it) */

      Project->x1 = Project->y1 = MAX_VB_ENTRY;
      Project->x2 = Project->y2 = MIN_VB_ENTRY;
    }
  }
}



/*****************************************************************************
*
* FUNCTION      : Project_Object
*
* ARGUMENTS     : Object   - Object to project
*                 Project  - Projection
*
* MODIFIED ARGS : Project
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Get the projection of a single object depending on the camera
*   used (perspective/orthographic).
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Project_Object(Object, Project)
OBJECT *Object;
PROJECT *Project;
{
  int infinite;
  DBL Volume;

  /* Init project fields, assuming the object is visible! */

  Project->x1 = Project->y1 = MIN_VB_ENTRY;
  Project->x2 = Project->y2 = MAX_VB_ENTRY;

  BOUNDS_VOLUME(Volume, Object->Bounds);

  infinite = (Volume > INFINITE_VOLUME);

  Get_Perspective_Projection(Object, Project, infinite);

  Print_Point(POINT_MOD);
}



/*****************************************************************************
*
* FUNCTION      : Project_Bounding_Slab
*
* ARGUMENTS     : Project  - Projection
*                 Tree     - Current node/leaf
*                 Object   - Node/leaf in bounding slab hierarchy
*
* MODIFIED ARGS : Project, Tree
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Project the bounding slab hierarchy onto the screen and thus create
*   the vista buffer hierarchy.
*
* CHANGES
*
*   -
*
******************************************************************************/

static void Project_Bounding_Slab(Project, Tree, Object)
PROJECT *Project;
PROJECT_TREE_NODE **Tree;
OBJECT *Object;
{
  unsigned short int i;
  COMPOSITE *Comp;
  PROJECT Temp;
  PROJECT_TREE_LEAF *Leaf;
  PROJECT_TREE_NODE New;

  if (Object->Type & BOUNDING_OBJECT)
  {
    /* Current object is a bounding object, i.e. a node in the slab tree */

    Comp = (COMPOSITE *)Object;

    /* First, init new entry */

    New.Entries = 0;

    New.Object = Object;

    New.Project.x1 = New.Project.y1 = MAX_VB_ENTRY;
    New.Project.x2 = New.Project.y2 = MIN_VB_ENTRY;

    for (i = 0; i < BUNCHING_FACTOR; i++)
    {
      New.Entry[i] = NULL;
    }

    /* This is no leaf */

    New.is_leaf = FALSE;

    /* Second, get new entry, i.e. project node's entries */

    for (i = 0; i < Comp->Entries; i++)
    {
      Project_Bounding_Slab(&Temp, &New.Entry[New.Entries], Comp->Objects[i]);

      /* Use only visible entries */

      if (New.Entry[New.Entries] != NULL)
      {
	New.Project.x1 = min(New.Project.x1, Temp.x1);
	New.Project.x2 = max(New.Project.x2, Temp.x2);
	New.Project.y1 = min(New.Project.y1, Temp.y1);
	New.Project.y2 = max(New.Project.y2, Temp.y2);

	New.Entries++;
      }
    }

    /* If the new entry is visible we'll use it */

    if (New.Entries > 0)
    {
      /* If there's only one entry, we won't need a new node. */

      if (New.Entries == 1)
      {
	*Tree    = New.Entry[0];
	*Project = New.Project;
      }
      else
      {
	/* Allocate memory for new node in the vista tree (never freed!)  */

	*Tree = (PROJECT_TREE_NODE *)VB_malloc(sizeof(PROJECT_TREE_NODE));

	if (*Tree == NULL)
	{
	  Fatal_MAError("vista tree node");
	}

	**Tree = New;

	*Project = New.Project;
      }
    }
  }
  else
  {
    /* Current object is a normal object, i.e. a leaf in the slab tree */

    /* Get object's projection */

    Project_Object(Object, Project);

    /* Is the object visible? */

    if ((Project->x1 <= Project->x2) && (Project->y1 <= Project->y2))
    {
      /* Allocate memory for new leaf in the vista tree (never freed!)  */

      *Tree = (PROJECT_TREE_NODE *)VB_malloc(sizeof(PROJECT_TREE_LEAF));

      if (*Tree == NULL)
      {
	Fatal_MAError("vista tree leaf");
      }

      /* Init new leaf */

      Leaf = (PROJECT_TREE_LEAF *)(*Tree);

      Leaf->Object = Object;

      Leaf->Project = *Project;

      /* Yes, this is a leaf */

      Leaf->is_leaf = TRUE;
    }
  }
}



/*****************************************************************************
*
* FUNCTION      : Build_Vista_Tree
*
* ARGUMENTS     : none
*
* MODIFIED ARGS : none
*
* RETURN VALUE  : none
*
* AUTHOR        : Dieter Bayer, May 1994
*
* DESCRIPTION
*
*   Build the vista tree, i.e. the 2d representation of the bounding slab
*   hierarchy in image space.
*
*   This only works for perspective and orthographic cameras.
*
* CHANGES
*
*   -
*
******************************************************************************/

void Build_Vista_Tree()
{
  PROJECT Project;

  Root_Vista = NULL;

  if (Use_Slabs)
  {
    fprintf(stderr, "Building vista buffer");

    Begin_Point();

    Project_Bounding_Slab(&Project, &Root_Vista, Root_Object);

    End_Point();
  }
}



#endif