/*======================================================================*/
/* 3DLib								*/
/*									*/
/* Rendering-related functions						*/
/*									*/
/* AUTHOR:	Gabor Nagy						*/
/* DATE:	1996-Jan-09 23:52:26					*/
/*									*/
/* 3DLib(TM) Copyright (C) 1995 by Gabor Nagy. All rights reserved.	*/
/*======================================================================*/
#include <assert.h>
#include <float.h>
#include <math.h>
#include <stdarg.h>		// For varargs
#include <stdio.h>
#include <stdlib.h>

#include <ELists.h>
#include <EMalloc.h>
#include <EStrings.h>

#ifndef _E3D_h
#include <E3D/E3D.h>
#endif

#ifndef _E3DMath_h
#include <E3D/Math.h>
#endif

#ifndef _E3DMatrix_h
#include <E3D/Matrix.h>
#endif

#ifndef _E3DRay_h
#include <E3D/Ray.h>
#endif

#ifndef _E3DScene_h
#include <E3D/Scene.h>
#endif


#define	E3dRTRIANGLE_ALLOC_STEP	32

#undef HUGE
#ifndef HUGE
#define HUGE 1e30
#endif


E3dRenderer**	E3d_Renderers=NULL;
int		E3d_NumOfRenderers=0;


//================================================
// Find a Renderer by its name
//================================================
E3dRenderer* E3d_RendererFindByName(char* LName)
{
 E3dRenderer**	LRenderers=E3d_Renderers;

 if(LRenderers)
 {
  E3dRenderer*	LRenderer;
  unsigned int	LC, LN;

  for(LC=0, LN=E3d_NumOfRenderers;LC<LN;LC++)
  {
   LRenderer=LRenderers[LC];
   if(EStr_StringsEqual(LRenderer->Name, LName)) return(LRenderer);
  }
 }
 return(NULL);
}


//================================================
// Register a new LRenderer
//================================================
E3dRenderer* E3d_RendererRegister(E3dRenderer* LRendererTemplate)
{
 E3dRenderer**	LRenderers=E3d_Renderers;
 E3dRenderer*	LRenderer;
 unsigned int	LC, LN;

 if(LRenderers)
 {
  for(LC=0, LN=E3d_NumOfRenderers;LC<LN;LC++)
  {
   LRenderer=LRenderers[LC];
   if(EStr_StringsEqual(LRenderer->Name, LRendererTemplate->Name))
   {
    LRenderer->StructSize=LRendererTemplate->StructSize;
    LRenderer->RenderProc=LRendererTemplate->RenderProc;

    return(LRenderer);
   }
  }
 }

 if(LRenderers==NULL) E3d_Renderers=LRenderers=(E3dRenderer**)EMalloc(sizeof(E3dRenderer*));
 else
 {
  if((LRenderers=(E3dRenderer**)ERealloc(E3d_Renderers, sizeof(E3dRenderer*)*(E3d_NumOfRenderers+1)))!=NULL) E3d_Renderers=LRenderers;
 }

 if(LRenderers)
 {
  if((LRenderer=(E3dRenderer*)EMalloc(sizeof(E3dRenderer)))!=NULL)
  {
   LRenderers[E3d_NumOfRenderers]=LRenderer;

   LRenderer->Name=EStrDup(LRendererTemplate->Name);
   LRenderer->StructSize=LRendererTemplate->StructSize;
   LRenderer->RenderProc=LRendererTemplate->RenderProc;
   E3d_NumOfRenderers++;
  }
  return(LRenderer);
 }
 else return(NULL);
}


//================================================
// Deactivate a Renderer
//================================================
void E3d_RendererDeactivate(E3dRenderer* LRenderer)
{
 E3dRenderer**	LRenderers=E3d_Renderers;

 if(LRenderers)
 {
  unsigned int	LC, LN;

  for(LC=0, LN=E3d_NumOfRenderers;LC<LN;LC++)
  {
   if(LRenderers[LC]==LRenderer)
   {
    LRenderer->RenderProc=NULL;
   }
  }
 }
}


//================================================
// Remove a Renderer
//================================================
void E3d_RendererRemove(E3dRenderer* LRenderer)
{
 E3dRenderer**	LRenderers=E3d_Renderers;

 if(LRenderers)
 {
  unsigned int	LC, LN;

  for(LC=0, LN=E3d_NumOfRenderers;LC<LN;LC++)
  {
   if(LRenderers[LC]==LRenderer)
   {
    if(LRenderer->Name) EFree(LRenderer->Name);
    EFree(LRenderer);
    if(LC<(LN-1)) memmove(LRenderer, LRenderer+1, sizeof(E3dRenderer*)*(LN-LC-1));
    E3d_NumOfRenderers-=1;
    if(E3d_NumOfRenderers==0) { EFree(E3d_Renderers);E3d_Renderers=NULL; }
    else E3d_Renderers=(E3dRenderer**)ERealloc(E3d_Renderers, sizeof(E3dRenderer*)*E3d_NumOfRenderers);
    return;
   }
  }
 }
}



#define _E3dM_SetRTriangleT(mRTriangle, mITriangle, mVertices)\
{\
 mRTriangle->Normal.X=LPolygon->Normal.X;mRTriangle->Normal.Y=LPolygon->Normal.Y;mRTriangle->Normal.Z=LPolygon->Normal.Z;\
 LVertexNodes=mITriangle->VertexNodes;\
 LVertexNode=LVertexNodes[0];LVertex=mVertices+LVertexNode->VertexID;\
 mRTriangle->V0[E3dX]=LVertex->X;mRTriangle->V0[E3dY]=LVertex->Y;mRTriangle->V0[E3dZ]=LVertex->Z;\
 mRTriangle->Normal0.X=LVertexNode->Normal.X;mRTriangle->Normal0.Y=LVertexNode->Normal.Y;mRTriangle->Normal0.Z=LVertexNode->Normal.Z;\
 mRTriangle->S[0]=LVertexNode->S;mRTriangle->T[0]=LVertexNode->T;\
\
 LVertexNode=LVertexNodes[1];LVertex=mVertices+LVertexNode->VertexID;\
 mRTriangle->V1[E3dX]=LVertex->X;mRTriangle->V1[E3dY]=LVertex->Y;mRTriangle->V1[E3dZ]=LVertex->Z;\
 mRTriangle->Normal1.X=LVertexNode->Normal.X;mRTriangle->Normal1.Y=LVertexNode->Normal.Y;mRTriangle->Normal1.Z=LVertexNode->Normal.Z;\
 mRTriangle->S[1]=LVertexNode->S;mRTriangle->T[1]=LVertexNode->T;\
\
 LVertexNode=LVertexNodes[2];LVertex=mVertices+LVertexNode->VertexID;\
 mRTriangle->V2[E3dX]=LVertex->X;mRTriangle->V2[E3dY]=LVertex->Y;mRTriangle->V2[E3dZ]=LVertex->Z;\
 mRTriangle->Normal2.X=LVertexNode->Normal.X;mRTriangle->Normal2.Y=LVertexNode->Normal.Y;mRTriangle->Normal2.Z=LVertexNode->Normal.Z;\
 mRTriangle->S[2]=LVertexNode->S;mRTriangle->T[2]=LVertexNode->T;\
 mRTriangle->PolyGroup=LPolyGroup;\
}


#define _E3dM_SetRTriangle(mRTriangle, mVertexNode, mVertices, mNormalX, mNormalY, mNormalZ)\
{\
 mRTriangle->Normal.X=mNormalX;mRTriangle->Normal.Y=mNormalY;mRTriangle->Normal.Z=mNormalZ;\
 LVertex=mVertices+mVertexNode->VertexID;\
 mRTriangle->V0[E3dX]=LVertex->X;mRTriangle->V0[E3dY]=LVertex->Y;mRTriangle->V0[E3dZ]=LVertex->Z;\
 mRTriangle->Normal0.X=mVertexNode->Normal.X;mRTriangle->Normal0.Y=mVertexNode->Normal.Y;mRTriangle->Normal0.Z=mVertexNode->Normal.Z;\
 mRTriangle->S[0]=mVertexNode->S;mRTriangle->T[0]=mVertexNode->T;\
 mVertexNode++;\
\
 LVertex=mVertices+mVertexNode->VertexID;\
 mRTriangle->V1[E3dX]=LVertex->X;mRTriangle->V1[E3dY]=LVertex->Y;mRTriangle->V1[E3dZ]=LVertex->Z;\
 mRTriangle->Normal1.X=mVertexNode->Normal.X;mRTriangle->Normal1.Y=mVertexNode->Normal.Y;mRTriangle->Normal1.Z=mVertexNode->Normal.Z;\
 mRTriangle->S[1]=mVertexNode->S;mRTriangle->T[1]=mVertexNode->T;\
 mVertexNode++;\
\
 LVertex=mVertices+mVertexNode->VertexID;\
 mRTriangle->V2[E3dX]=LVertex->X;mRTriangle->V2[E3dY]=LVertex->Y;mRTriangle->V2[E3dZ]=LVertex->Z;\
 mRTriangle->Normal2.X=mVertexNode->Normal.X;mRTriangle->Normal2.Y=mVertexNode->Normal.Y;mRTriangle->Normal2.Z=mVertexNode->Normal.Z;\
 mRTriangle->S[2]=mVertexNode->S;mRTriangle->T[2]=mVertexNode->T;\
 mRTriangle->PolyGroup=LPolyGroup;\
}




//================================================================
// Get renderable triangles for a PolyGroup
//================================================================
int E3d_PolyGroupCountTriangles(E3dPolyGroup* LPolyGroup)
{
 E3dPolygon*		LPolygon;
 unsigned int		LNumOfRTriangles=0,
			LPolyCnt=LPolyGroup->NumOfPolygons;


 LPolygon=LPolyGroup->Polygons;
 do
 {
// Is this Polygon a triangle?
//
  if(LPolygon->NumOfVertices==3) LNumOfRTriangles++;
  else
  {
   if(LPolygon->NumOfTriangles) LNumOfRTriangles+=LPolygon->NumOfTriangles;
   else
   {
// If the Polygon is not triangulated for OpenGL (it was convex with no holes),
// we have to triangulate it here
//
//    LNumOfRTriangles+=E3d_TriangulatePolygonForRendering(LVertices, LNewVertices, LPolygon, &LRTriangles, LNumOfRTriangles, &LMaxNumOfRTriangles);
   }
  }

  LPolygon++;
 } while(--LPolyCnt);

 return(LNumOfRTriangles);
}


//================================================================
// Get renderable triangles for a PolyGroup
//================================================================
int E3d_PolyGroupTriangulateForRendering(E3dPolyGroup* LPolyGroup, E3dVertex* LVertices, E3dVertex* LNewVertices, E3dRTriangle** LRTrianglesPtr)
{
 E3dVertex*		LVertex;
 E3dPolygon*		LPolygon;
 E3dITriangle*		LITriangle;
 E3dRTriangle*		LRTriangles=NULL;
 E3dRTriangle*		LRTrianglesT;
 E3dRTriangle*		LRTriangle;
 E3dVertexNode*		LVertexNode;
 E3dVertexNode**	LVertexNodes;
 unsigned int		LC, LN, LNumOfRTriangles=0, LMaxNumOfRTriangles=0,
			LAllocStep,
			LPolyCnt=LPolyGroup->NumOfPolygons;


 LPolygon=LPolyGroup->Polygons;
 do
 {
  LVertexNode=LPolygon->VertexNodes;

  if(LVertexNode)
  {
// Is this Polygon a triangle?
//
   if(LPolygon->NumOfVertices==3)
   {
// Allocate new RTriangle if necessary
//
    if(LRTriangles==NULL)
    {
     if((LRTriangles=(E3dRTriangle*)EMalloc(sizeof(E3dRTriangle)*E3dRTRIANGLE_ALLOC_STEP))!=NULL)
     {
      LMaxNumOfRTriangles=E3dRTRIANGLE_ALLOC_STEP;
      LRTriangle=LRTriangles;
     }
    }
    else
    {
     if(LNumOfRTriangles>=LMaxNumOfRTriangles)
     {
      if((LRTrianglesT=(E3dRTriangle*)ERealloc(LRTriangles, sizeof(E3dRTriangle)*(LNumOfRTriangles+E3dRTRIANGLE_ALLOC_STEP)))!=NULL)
      {
       LRTriangles=LRTrianglesT;LMaxNumOfRTriangles+=E3dRTRIANGLE_ALLOC_STEP;
      }
     }
    }

    LRTriangle=LRTriangles+LNumOfRTriangles;
    LNumOfRTriangles++;

    _E3dM_SetRTriangle(LRTriangle, LVertexNode, LNewVertices, LPolygon->Normal.X, LPolygon->Normal.Y, LPolygon->Normal.Z);
    LRTriangle->Face=LPolygon;
   }
   else
   {
// If the Polygon had been triangulated, we can use those triangles
//
    if((LITriangle=LPolygon->ITriangles)!=NULL)
    {
     LN=LPolygon->NumOfTriangles;
     if(LN > E3dRTRIANGLE_ALLOC_STEP) LAllocStep = LN;
     else LAllocStep = E3dRTRIANGLE_ALLOC_STEP;

// Allocate new RTriangles if necessary
//
     if(LRTriangles == NULL)
     {
      if((LRTriangles = (E3dRTriangle*)EMalloc(sizeof(E3dRTriangle) * LAllocStep))!=NULL)
      {
       LMaxNumOfRTriangles = LAllocStep;
       LRTriangle = LRTriangles;
      }
     }
     else
     {
// Make sure we have enough triangles allocated
//
      if((LMaxNumOfRTriangles-LNumOfRTriangles) < LN)
      {
       if((LRTrianglesT = (E3dRTriangle*)ERealloc(LRTriangles, sizeof(E3dRTriangle)*(LNumOfRTriangles + LN)))!=NULL)
       {
	LRTriangles = LRTrianglesT;LMaxNumOfRTriangles = LNumOfRTriangles + LN;
       }
      }
     }

     LRTriangle = LRTriangles + LNumOfRTriangles;
     LNumOfRTriangles += LN;

// Go through the triangles of this Polygon
//
     for(LC=0;LC<LN;LC++, LITriangle++, LRTriangle++)
     {

      LVertexNodes=LITriangle->VertexNodes;
      _E3dM_SetRTriangleT(LRTriangle, LITriangle, LNewVertices);
      LRTriangle->Face=LPolygon;
     }

    }
    else
    {
// If the Polygon is not triangulated for OpenGL (it was convex with no holes),
// we have to triangulate it here
//
     LNumOfRTriangles+=E3d_TriangulatePolygonForRendering(LVertices, LNewVertices, LPolygon, &LRTriangles, LNumOfRTriangles, &LMaxNumOfRTriangles);
    }

   }
  }

  LPolygon++;
 } while(--LPolyCnt);

 if(LRTrianglesPtr) *LRTrianglesPtr=LRTriangles;
 return(LNumOfRTriangles);
}



//========================================================================
// Get the local bounding box of an RTriangleGroup used for rendering
//========================================================================
EBool E3d_RenderPolyGroupGetBoundingBox(E3dPolyGroup* LRenderPolyGroup, E3d3DPosition* LBBMin, E3d3DPosition* LBBMax)
{
 E3dRTriangle*	LRTriangle=LRenderPolyGroup->Triangles;
 E3dCoordinate*	LRVertex;
 unsigned int	LC, LN=LRenderPolyGroup->NumOfTriangles;


 if(LRTriangle)
 {
// Must initialize BBox on first vertex
//
  LRVertex=LRTriangle->V0;
  LBBMin->X = LBBMax->X = LRVertex[E3dX];
  LBBMin->Y = LBBMax->Y = LRVertex[E3dY];
  LBBMin->Z = LBBMax->Z = LRVertex[E3dZ];

  LRVertex=LRTriangle->V1;
  if(LRVertex[E3dX] < LBBMin->X) LBBMin->X = LRVertex[E3dX];
  if(LRVertex[E3dX] > LBBMax->X) LBBMax->X = LRVertex[E3dX];
  if(LRVertex[E3dY] < LBBMin->Y) LBBMin->Y = LRVertex[E3dY];
  if(LRVertex[E3dY] > LBBMax->Y) LBBMax->Y = LRVertex[E3dY];
  if(LRVertex[E3dZ] < LBBMin->Z) LBBMin->Z = LRVertex[E3dZ];
  if(LRVertex[E3dZ] > LBBMax->Z) LBBMax->Z = LRVertex[E3dZ];

  LRVertex=LRTriangle->V2;
  if(LRVertex[E3dX] < LBBMin->X) LBBMin->X = LRVertex[E3dX];
  if(LRVertex[E3dX] > LBBMax->X) LBBMax->X = LRVertex[E3dX];
  if(LRVertex[E3dY] < LBBMin->Y) LBBMin->Y = LRVertex[E3dY];
  if(LRVertex[E3dY] > LBBMax->Y) LBBMax->Y = LRVertex[E3dY];
  if(LRVertex[E3dZ] < LBBMin->Z) LBBMin->Z = LRVertex[E3dZ];
  if(LRVertex[E3dZ] > LBBMax->Z) LBBMax->Z = LRVertex[E3dZ];
  LRTriangle++;


  for(LC=1;LC<LN;LC++)
  {
   LRVertex=LRTriangle->V0;
   if(LRVertex[E3dX] < LBBMin->X) LBBMin->X = LRVertex[E3dX];
   if(LRVertex[E3dX] > LBBMax->X) LBBMax->X = LRVertex[E3dX];
   if(LRVertex[E3dY] < LBBMin->Y) LBBMin->Y = LRVertex[E3dY];
   if(LRVertex[E3dY] > LBBMax->Y) LBBMax->Y = LRVertex[E3dY];
   if(LRVertex[E3dZ] < LBBMin->Z) LBBMin->Z = LRVertex[E3dZ];
   if(LRVertex[E3dZ] > LBBMax->Z) LBBMax->Z = LRVertex[E3dZ];

   LRVertex=LRTriangle->V1;
   if(LRVertex[E3dX] < LBBMin->X) LBBMin->X = LRVertex[E3dX];
   if(LRVertex[E3dX] > LBBMax->X) LBBMax->X = LRVertex[E3dX];
   if(LRVertex[E3dY] < LBBMin->Y) LBBMin->Y = LRVertex[E3dY];
   if(LRVertex[E3dY] > LBBMax->Y) LBBMax->Y = LRVertex[E3dY];
   if(LRVertex[E3dZ] < LBBMin->Z) LBBMin->Z = LRVertex[E3dZ];
   if(LRVertex[E3dZ] > LBBMax->Z) LBBMax->Z = LRVertex[E3dZ];

   LRVertex=LRTriangle->V2;
   if(LRVertex[E3dX] < LBBMin->X) LBBMin->X = LRVertex[E3dX];
   if(LRVertex[E3dX] > LBBMax->X) LBBMax->X = LRVertex[E3dX];
   if(LRVertex[E3dY] < LBBMin->Y) LBBMin->Y = LRVertex[E3dY];
   if(LRVertex[E3dY] > LBBMax->Y) LBBMax->Y = LRVertex[E3dY];
   if(LRVertex[E3dZ] < LBBMin->Z) LBBMin->Z = LRVertex[E3dZ];
   if(LRVertex[E3dZ] > LBBMax->Z) LBBMax->Z = LRVertex[E3dZ];
   LRTriangle++;
  }

  LRenderPolyGroup->BBoxMin.X=LBBMin->X;
  LRenderPolyGroup->BBoxMin.Y=LBBMin->Y;
  LRenderPolyGroup->BBoxMin.Z=LBBMin->Z;
  LRenderPolyGroup->BBoxMax.X=LBBMax->X;
  LRenderPolyGroup->BBoxMax.Y=LBBMax->Y;
  LRenderPolyGroup->BBoxMax.Z=LBBMax->Z;

  return(TRUE);
 }

 return(FALSE);
}



//================================================================================
// Get the bounding box of a view-transformed scene prepared for rendering
//================================================================================
EBool E3d_SceneGetRenderBoundingBox(E3dScene* LScene)
{
 E3dPolyGroup**		LRenderPolyGroups=LScene->RenderPolyGroups;
 E3dPolyGroup*		LRenderPolyGroup;
 E3d3DPosition		LBBMin, LBBMax;
 E3d3DPosition		LTBBMin, LTBBMax;
 E3dCoordinate		LCenterX, LCenterY, LCenterZ;
 unsigned int		LC, LN=LScene->NumOfRenderPolyGroups;
 EBool			LFirst=TRUE, LResult=FALSE;


 for(LC=0;LC<LN;LC++)
 {
  LRenderPolyGroup=LRenderPolyGroups[LC];
  if(LFirst)
  {
   E3d_RenderPolyGroupGetBoundingBox(LRenderPolyGroup, &LBBMin, &LBBMax);
   LResult=TRUE;
   LFirst=FALSE;
  }
  else
  {
   if(E3d_RenderPolyGroupGetBoundingBox(LRenderPolyGroup, &LTBBMin, &LTBBMax))
   {
    if(LTBBMin.X < LBBMin.X) LBBMin.X = LTBBMin.X;
    if(LTBBMax.X > LBBMax.X) LBBMax.X = LTBBMax.X;

    if(LTBBMin.Y < LBBMin.Y) LBBMin.Y = LTBBMin.Y;
    if(LTBBMax.Y > LBBMax.Y) LBBMax.Y = LTBBMax.Y;

    if(LTBBMin.Z < LBBMin.Z) LBBMin.Z = LTBBMin.Z;
    if(LTBBMax.Z > LBBMax.Z) LBBMax.Z = LTBBMax.Z;
    LResult=TRUE;
   }
  }
 }

// Expand BBox by 1% and add a margin (in case one of the dimensions is 0)
//
 {
  E3dCoordinate	LScale=1.01, LMargin=0.1;

// Must scale around centroid of the BBox
//
  LCenterX=(LBBMax.X-LBBMin.X)*0.5+LBBMin.X;
  LCenterY=(LBBMax.Y-LBBMin.Y)*0.5+LBBMin.Y;
  LCenterZ=(LBBMax.Z-LBBMin.Z)*0.5+LBBMin.Z;

  LBBMin.X=(LBBMin.X - LCenterX)*LScale + LCenterX - LMargin;
  LBBMin.Y=(LBBMin.Y - LCenterY)*LScale + LCenterY - LMargin;
  LBBMin.Z=(LBBMin.Z - LCenterZ)*LScale + LCenterZ - LMargin;
  LBBMax.X=(LBBMax.X - LCenterX)*LScale + LCenterX + LMargin;
  LBBMax.Y=(LBBMax.Y - LCenterY)*LScale + LCenterY + LMargin;
  LBBMax.Z=(LBBMax.Z - LCenterZ)*LScale + LCenterZ + LMargin;
 }

 LScene->BBoxMin.X=LBBMin.X;
 LScene->BBoxMin.Y=LBBMin.Y;
 LScene->BBoxMin.Z=LBBMin.Z;
 LScene->BBoxMax.X=LBBMax.X;
 LScene->BBoxMax.Y=LBBMax.Y;
 LScene->BBoxMax.Z=LBBMax.Z;

 return(LResult);
}


void E3d_SetupTriangleGroup(E3dPolyGroup* LRenderPolyGroup, E3dMatrix LNormalMatrix)
{
 E3dRTriangle*	LRTriangle=LRenderPolyGroup->Triangles;
 E3dCoordinate*	LV0;
 E3dCoordinate*	LV1;
 E3dCoordinate*	LV2;
 E3dCoordinate	mX, mY, mZ;
 unsigned int	LPolyCnt=LRenderPolyGroup->NumOfTriangles;


// Build bounding box for this TriangleGroup
//
 E3d_RenderPolyGroupGetBoundingBox(LRenderPolyGroup, &(LRenderPolyGroup->BBoxMin), &(LRenderPolyGroup->BBoxMax));


// Transform Triangle and Vertex normals
//
 do
 {
// Precompute vectors: V0->V1 and V0->V2 to speed up the ray/triangle intersection tests
//
  LV0=LRTriangle->V0;
  LV1=LRTriangle->V1;
  LV2=LRTriangle->V2;
  LRTriangle->Edge1[E3dX]=LV1[E3dX]-LV0[E3dX];
  LRTriangle->Edge1[E3dY]=LV1[E3dY]-LV0[E3dY];
  LRTriangle->Edge1[E3dZ]=LV1[E3dZ]-LV0[E3dZ];
  LRTriangle->Edge2[E3dX]=LV2[E3dX]-LV0[E3dX];
  LRTriangle->Edge2[E3dY]=LV2[E3dY]-LV0[E3dY];
  LRTriangle->Edge2[E3dZ]=LV2[E3dZ]-LV0[E3dZ];

  mX=LRTriangle->Normal.X;mY=LRTriangle->Normal.Y;mZ=LRTriangle->Normal.Z;
  E3dM_MatrixTransform3x3(LNormalMatrix, LRTriangle->Normal.X, LRTriangle->Normal.Y, LRTriangle->Normal.Z);

  mX=LRTriangle->Normal0.X;mY=LRTriangle->Normal0.Y;mZ=LRTriangle->Normal0.Z;
  E3dM_MatrixTransform3x3(LNormalMatrix, LRTriangle->Normal0.X, LRTriangle->Normal0.Y, LRTriangle->Normal0.Z);

  mX=LRTriangle->Normal1.X;mY=LRTriangle->Normal1.Y;mZ=LRTriangle->Normal1.Z;
  E3dM_MatrixTransform3x3(LNormalMatrix, LRTriangle->Normal1.X, LRTriangle->Normal1.Y, LRTriangle->Normal1.Z);

  mX=LRTriangle->Normal2.X;mY=LRTriangle->Normal2.Y;mZ=LRTriangle->Normal2.Z;
  E3dM_MatrixTransform3x3(LNormalMatrix, LRTriangle->Normal2.X, LRTriangle->Normal2.Y, LRTriangle->Normal2.Z);

  LRTriangle->LastRayNumber=-1;
  LRTriangle++;
 } while(--LPolyCnt);
}


//========================================================
// http://www.cbloom.com/3d/techdocs/fast_ray_tri.txt
//========================================================
void FastTriPreproc(E3dRTriangle* LRTriangle)
{
/*
 edge1 = tri.vert1 - tri.vert0;
 edge2 = tri.vert2 - tri.vert1;
 edgePerp1 = tri.normal CROSS edge1
 edgePerp2 = tri.normal CROSS edge2
 edgePlane1 = { edgePerp1, edgePerp1*tri.vert0 }
 edgePlane2 = { edgePerp2, edgePerp2*tri.vert1 }

// A plane is just a vector and a constant; we can scale this
// normal to change the distances we want, so we make
//
 baryPlane1 = edgePlane1 / (edgePlane1.Distance(tri.vert2))
 baryPlane2 = edgePlane2 / (edgePlane1.Distance(tri.vert0))
*/


}


/*
//========================================================
// http://www.cbloom.com/3d/techdocs/fast_ray_tri.txt
//========================================================
EBool RayTriIntersect(const Triangle & tri,const Vec3 & LFrom, const Vec3 & LTo)
{
 float dTo = (plane.normal * LTo - plane.offset);	// Distance between "LTo" and the triangle's plane


 if(dTo <= 0) return(FALSE);		// Segment doesn't reach back side of triangle

 {
  float dFm = (plane.normal * LFm - plane.offset);	// Distance between "LFm" and the triangle's plane

  if(dFm > 0) return(FALSE);		// Segment starts already on the back side

  assert( dTo > 0 && dFm <= 0 );

  float denom = dTo - dFm;

  assert( denom > 0 );

  Vec3 temp = dTo * from + dFm * to;

  float uTimesDenom = (tri.bary1.vector DOT temp) - tri.bary1.offset * denom;

  if(uTimesDenom < 0 || uTimesDenom > denom) return(FALSE);	// Off triangle

  float vTimesDenom = (tri.bary2.vector DOT temp) - tri.bary2.offset * denom;

  if(vTimesDenom < 0 || (uTimesDenom+vTimesDenom) > denom) return(FALSE);	// Off triangle


// It's on the triangle, compute hit info :
//
  {
   float inv = 1 .0/ denom;
   *pT = dFm * inv;
   *pU = uTimesDenom * inv;
   *pV = vTimesDenom * inv;
  }
 }
 return(TRUE);
}
*/


//================================================================
// Get camera-transformed triangles from Meshes and Faces
// If LCamera is NULL, the result will be world-transformed.
//================================================================
unsigned int E3d_SceneGetRenderTriangles(E3dScene* LScene, E3dRenderInfo* LRenderInfo, E3dCamera* LCamera)
{
 E3dModel**		LRootModels;
 E3dModel*		LModel;
 E3dPolyGroup**		LPolyGroups;
 E3dPolyGroup*		LPolyGroup;
 E3dRTriangle*		LRTriangles;

 E3dMatrix		LMatrixVtx, LNormalMatrix;

 E3dGeometry**		LGeometries;
 E3dGeometry*		LGeometry;
 E3dVertex*		LNewVertices;
 E3dMesh*		LMesh;
 E3dVertexNode*		LVertexNodes;
 E3dCoordinate		mX, mY, mZ;
 E3dCoordinate		LSizeFactor=LRenderInfo->SizeFactor/200.0;	// For LODs
 unsigned int		LNumOfVertices;
 unsigned int		LNumOfRootModels, LGmCnt, LGmNum, LGCnt;
 unsigned int		LHC, LC,
			LNumOfRTriangles,
			LTotalNumOfPolygons=0, LTotalNumOfTriangles=0,
			LNumOfRenderPolyGroups;
 int			LIndex;
 EBool			LUseLODs=FALSE;


// If the Scene already has RenderPolyGroups cached, free them
//
 E3d_SceneFreeRenderTriangles(LScene);


// Transformed and triangulated Geometries
//
 LNumOfRenderPolyGroups=0;

 LTotalNumOfPolygons=0;LTotalNumOfTriangles=0;

 LNumOfRootModels=LScene->NumOfRootModels;LRootModels=LScene->RootModels;
 for(LHC=0;LHC<LNumOfRootModels;LHC++)
 {
  LModel=LRootModels[LHC];

  for(;LModel;LModel=LModel->Next)
  {
//    if(LModel->Selection != E3dSEL_NONE)
   if(LModel->Visible)
   {
    switch(LModel->Type)
    {
     case E3dMDL_NORMAL:
      if(LCamera)
      {
// Create local-to-viewer matrices to transform vertices and normals
//
       E3d_MatrixCopy(LCamera->WorldToViewerMatrix, LMatrixVtx);
       E3d_MatrixCopy(LCamera->WorldToViewerMatrix, LNormalMatrix);
       E3d_MatrixMult(LModel->LocalToWorldMatrix, LMatrixVtx);
       E3d_MatrixMult(LModel->LocalToWorldMatrix, LNormalMatrix);
      }
      else
      {
// If LCamera is NULL, transform to World coordinates
//
       E3d_MatrixCopy(LModel->LocalToWorldMatrix, LMatrixVtx);
       E3d_MatrixCopy(LModel->LocalToWorldMatrix, LNormalMatrix);
      }

      LNormalMatrix[M30]=0.0;LNormalMatrix[M31]=0.0;LNormalMatrix[M32]=0.0;


      LGmNum=LModel->NumOfGeometries;LGeometries=LModel->Geometries;
      for(LGmCnt=0;LGmCnt<LGmNum;LGmCnt++)
      {
       LGeometry=LGeometries[LGmCnt];

       if(LGeometry->Visible)
       {
        if(LUseLODs) { E3dM_GeometryGetLODForDrawing(LGeometry); }

	switch(LGeometry->GeoType)
	{
	 caseE3dMESH():
	  LMesh=(E3dMesh*)LGeometry;

	  LNumOfVertices=LMesh->NumOfVertices;

	  if((LNewVertices=E3d_VerticesAllocate(LNumOfVertices, TRUE))!=NULL)
	  {
	   E3dVertex*		LVertex;
	   E3dVertex*		LVertices;
	   E3dVertex*		LNewVertex;


	   LVertices=LMesh->Vertices;

// Transform Mesh Vertices into the selected coordinate system
//
	   LVertex=LMesh->Vertices;
	   LNewVertex=LNewVertices;
	   for(LC=0;LC<LNumOfVertices;LC++, LVertex++, LNewVertex++)
	   {
	    mX=LVertex->X;mY=LVertex->Y;mZ=LVertex->Z;
	    E3dM_MatrixTransform3x4(LMatrixVtx, LNewVertex->X, LNewVertex->Y, LNewVertex->Z);
	   }

	   LPolyGroups=LMesh->PolyGroups;

	   for(LGCnt=0;LGCnt<LMesh->NumOfPolyGroups;LGCnt++)
	   {
	    LPolyGroup=LPolyGroups[LGCnt];

	    if((LPolyGroup->Visible)&&(LPolyGroup->NumOfPolygons>0))
	    {
	     if(ELst_AddPointerChk((void***)(&(LScene->RenderPolyGroups)), &(LScene->NumOfRenderPolyGroups), LPolyGroup))
	     {
	      LPolyGroup->RefCnt+=1;
	     }
	     else
	     {
// The PolyGroup is already in the list, so it must belong to an instanced Mesh
//
	      E3dPolyGroup*	LNewPolyGroup=E3d_PolyGroupClone(LPolyGroup, 0);

	      LNewPolyGroup->DrawingMaterial=LPolyGroup->DrawingMaterial;
	      LNewPolyGroup->DrawingMaterial->RefCnt+=1;


	      if(ELst_AddPointer((void***)(&(LScene->RenderPolyGroups)), &(LScene->NumOfRenderPolyGroups), LNewPolyGroup))
	      {
	       LNewPolyGroup->RefCnt+=1;
	      }
	      LPolyGroup=LNewPolyGroup;
	     }


	     if(LPolyGroup)
	     {
// TrianglesBase tells where the RTriangles of this PolyGroup start in a global RTriangles array.
// These global RTriangles arrays might be used by tools that need to store more information on
// triangles than the internal structures do.
// TrianglesBase should be computed by these tools, so we'll initialize them to "NOT_SET".
//
//
	      LPolyGroup->TrianglesBase=-1;
	      LPolyGroup->PolygonsBase=-1;

	      LNumOfRTriangles=E3d_PolyGroupTriangulateForRendering(LPolyGroup, LVertices, LNewVertices, &LRTriangles);

	      LTotalNumOfTriangles+=LNumOfRTriangles;
	      LTotalNumOfPolygons+=LPolyGroup->NumOfPolygons;
	      LPolyGroup->NumOfTriangles=LNumOfRTriangles;
	      LPolyGroup->Triangles=LRTriangles;

	      E3d_SetupTriangleGroup(LPolyGroup, LNormalMatrix);
	     }
	    }

	   }


	   EFree(LNewVertices);
	  }
	 break;

	 case E3dGEO_FACE:
	  {
	   E3dFace*	LFace=(E3dFace*)LGeometry;

	   if(LFace->NumOfTriangles==0) E3d_FaceTesselate(LFace);
	   LNumOfVertices=LFace->NumOfVertices;

	   if((LNewVertices=E3d_VerticesAllocate(LNumOfVertices, TRUE))!=NULL)
	   {
	    E3dVertex*		LVertex;
	    E3dVertex*		LVertices;
	    E3dVertex*		LNewVertex;

	    LVertices=LFace->Vertices;

// Transform Face vertices into the selected coordinate system
//
	    LVertex=LFace->Vertices;
	    LNewVertex=LNewVertices;
	    for(LC=0;LC<LNumOfVertices;LC++, LVertex++, LNewVertex++)
	    {
	     mX=LVertex->X;mY=LVertex->Y;mZ=LVertex->Z;
	     E3dM_MatrixTransform3x4(LMatrixVtx, LNewVertex->X, LNewVertex->Y, LNewVertex->Z);
	    }

	    if((LNumOfRTriangles=LFace->NumOfTriangles)>0)
	    {
	     E3dTriangle*	LTriangle=LFace->Triangles;
	     E3dRTriangle*	LRTriangle;
	     unsigned int	LTC;

// Allocate a fake PolyGroup for this Face
//
	     LPolyGroup=E3d_PolyGroupAllocate();

	     if(ELst_AddPointer((void***)(&(LScene->RenderPolyGroups)), &(LScene->NumOfRenderPolyGroups), LPolyGroup))
	     {
	      LPolyGroup->RefCnt+=1;
	     }

	     LRTriangles=LRTriangle=(E3dRTriangle*)EMalloc(sizeof(E3dRTriangle)*LNumOfRTriangles);

// Go through the triangles of this Face
//
	     for(LTC=0;LTC<LNumOfRTriangles;LTC++, LTriangle++, LRTriangle++)
	     {
	      LVertexNodes=LTriangle->VertexNodes;
	      _E3dM_SetRTriangle(LRTriangle, LVertexNodes, LNewVertices, LFace->Normal.X, LFace->Normal.Y, LFace->Normal.Z);
	      LRTriangle->Face=LFace;
	     }

	     LTotalNumOfTriangles+=LNumOfRTriangles;
	     LTotalNumOfPolygons++;			// A Face is just a Polygon with optional holes...

	     LPolyGroup->Flags=E3dPOLY_NORMALS_UPTODATE;
	     LPolyGroup->NumOfTriangles=LNumOfRTriangles;
	     LPolyGroup->Triangles=LRTriangles;

// TrianglesBase tells where the RTriangles of this PolyGroup start in a global RTriangles array.
// These global RTriangles arrays might be used by tools that need to store more information on
// triangles than the internal structures do.
// TrianglesBase should be computed by these tools, so we'll initialize them to "NOT_SET".
//
//
	     LPolyGroup->TrianglesBase=-1;
	     LPolyGroup->PolygonsBase=-1;

	     LPolyGroup->DrawingMaterial=LFace->DrawingMaterial;if(LFace->DrawingMaterial) LFace->DrawingMaterial->RefCnt+=1;


	     E3d_SetupTriangleGroup(LPolyGroup, LNormalMatrix);

//Printf("Face %p LPolyGroup %p  Tris %d %p\n", LFace, LPolyGroup, LPolyGroup->NumOfTriangles, LPolyGroup->Triangles);fflush(stdout);
	    }

	    EFree(LNewVertices);
	   }
	  }
	 break;
	}

       }
      }

     break;

/*
     E3dMDL_LIGHT:
//	printf("Light\n");
     break;
*/

     default:
     break;
    }


   }
  }
 }

printf("Polygons: %d  Triangles: %d\n", LTotalNumOfPolygons, LTotalNumOfTriangles);fflush(stdout);

 return(LScene->NumOfRenderPolyGroups);
}


//================================================
// Free renderable triangles of a Scene
//================================================
EBool E3d_SceneFreeRenderTriangles(E3dScene* LScene)
{
 if(LScene->RenderPolyGroups)
 {
  E3dPolyGroup**	LRenderPolyGroups=LScene->RenderPolyGroups;
  E3dPolyGroup*		LRenderPolyGroup;
  unsigned int		LC, LNumOfRenderPolyGroups=LScene->NumOfRenderPolyGroups;

  for(LC=0;LC<LNumOfRenderPolyGroups;LC++)
  {
   LRenderPolyGroup=LRenderPolyGroups[LC];
   E3d_PolyGroupFreeRenderTriangles(LRenderPolyGroup);

   E3d_PolyGroupFree(LRenderPolyGroup);
  }

  EFree(LRenderPolyGroups);
  LScene->RenderPolyGroups=NULL;
  LScene->NumOfRenderPolyGroups=0;

  return(TRUE);
 }

 return(FALSE);
}


//========================================================
// Add an E3dRTriangle to the Voxels it occupies
//========================================================
void E3d_GridAddTriangleNN(E3dVoxel* LVoxels, E3dRTriangle* LRTriangle, E3dPolyGroup* LGroup, E3d3DPosition* LSceneBBMin, E3dCoordinate LVoxelXSizeInv, E3dCoordinate LVoxelYSizeInv, E3dCoordinate LVoxelZSizeInv, int LGridXCount, int LGridYCount, int LGridZCount)
{
/*
 E3dVoxel*	LVoxel;
 E3dCoordinate*	LV0=LRTriangle->V0;	// The triangle's Vertices
 E3dCoordinate*	LV1=LRTriangle->V1;
 E3dCoordinate*	LV2=LRTriangle->V2;
 E3dCoordinate	LX0, LY0, LZ0, LX1, LY1, LZ1;
 E3dCoordinate	LEdgeLen, LDX, LDY, LDZ, LF;

 unsigned int	LC, LNumOfChords;
 int		LPageSize = LGridXCount*LGridYCount;
 int		LRowSize = LGridXCount;
 int		LX, LY, LZ;


// Get the length of the V1-V2 edge
//
 LX1=LV1[E3dX];
 LY1=LV1[E3dY];
 LZ1=LV1[E3dZ];
 LDX=LV2[E3dX]-LX1;
 LDY=LV2[E3dY]-LY1;
 LDZ=LV2[E3dZ]-LZ1;
 LEdgeLen=sqrt(LDX*LDX + LDY*LDY + LDZ*LDZ);
 LNumOfChords=((int)LEdgeLen)*8;
 LF=1.0/((E3dCoordinate)LNumOfChords);
 LDX=LDX*LF*LVoxelXSizeInv;
 LDY=LDY*LF*LVoxelYSizeInv;
 LDZ=LDZ*LF*LVoxelZSizeInv;


// Do a DDA from V0 to the points on the V1-V2 edge
// to walk the Voxels this Triangle occupies
//
 LX0=(LV0[E3dX]-LSceneBBMin->X)*LVoxelXSizeInv;
 LY0=(LV0[E3dY]-LSceneBBMin->Y)*LVoxelYSizeInv;
 LZ0=(LV0[E3dZ]-LSceneBBMin->Z)*LVoxelZSizeInv;
 LX1=(LX1-LSceneBBMin->X)*LVoxelXSizeInv;
 LY1=(LY1-LSceneBBMin->Y)*LVoxelYSizeInv;
 LZ1=(LZ1-LSceneBBMin->Z)*LVoxelZSizeInv;



 for(LC=0;LC<LNumOfChords;LC++)
 {
  E3d_RayDDAInitRb(LX0, LY0, LZ0, LX1, LY1, LZ1);

  while(E3d_DDANextVoxelRb(&LX, &LY, &LZ))
  {
// Make sure we're inside the Grid
//
   if((LX>=0)&&(LX<LGridXCount)&&(LY>=0)&&(LY<LGridYCount)&&(LZ>=0)&&(LZ<LGridZCount))
   {
    LVoxel=LVoxels + LZ*LPageSize + LY*LRowSize + LX;
    if(ELst_AddPointerChk((void***)(&(LVoxel->Triangles)), &(LVoxel->NumOfTriangles), LRTriangle))
    {
     if(LVoxel->Groups==NULL) LVoxel->Groups=(E3dPolyGroup**)EMalloc(sizeof(E3dPolyGroup*));
     else LVoxel->Groups=(E3dPolyGroup**)ERealloc(LVoxel->Groups, sizeof(E3dPolyGroup*)*LVoxel->NumOfTriangles);
     LVoxel->Groups[LVoxel->NumOfTriangles-1]=LGroup;
    }
   }
  }

// Next chord
//
  LX1+=LDX;
  LY1+=LDY;
  LZ1+=LDZ;
 }
*/
}




EBool E3d_RayIntersectBoundingBox(E3dCoordinate* origin, E3dCoordinate* dir, E3dCoordinate* minB, E3dCoordinate* maxB, E3dCoordinate* LTRet)
{
#define NUMDIM  3
#define RIGHT   0
#define LEFT    1
#define MIDDLE  2

 char		inside = TRUE;
 char		quadrant[NUMDIM];
 register int	i;
 int		whichPlane;
 E3dCoordinate	maxT[NUMDIM];
 E3dCoordinate	candidatePlane[NUMDIM];

// Find candidate planes; this loop can be avoided if
// rays cast all from the eye(assume perpsective view)
//
 for (i=0; i<NUMDIM; i++)
 {
   if(origin[i] < minB[i])
   {
     quadrant[i] = LEFT;
     candidatePlane[i] = minB[i];
     inside = FALSE;
   } else if(origin[i] > maxB[i])
	  {
	    quadrant[i] = RIGHT;
	    candidatePlane[i] = maxB[i];
	    inside = FALSE;
   } else quadrant[i] = MIDDLE;
 }


// Ray origin inside bounding box
//
 if(inside) return(TRUE);


// Calculate T distances to candidate planes
//
 for(i = 0; i < NUMDIM; i++)
 {
  if(quadrant[i] != MIDDLE && dir[i] !=0.) maxT[i] = (candidatePlane[i]-origin[i]) / dir[i];
  else maxT[i] = -1.0;
 }

// Get largest of the maxT's for final choice of intersection
//
 whichPlane = 0;
 for (i = 1; i < NUMDIM; i++) if(maxT[whichPlane] < maxT[i]) whichPlane = i;

// Check if final candidate is actually inside the box
//
 if(maxT[whichPlane] < 0.) return(FALSE);

 *LTRet=maxT[whichPlane];

 return(TRUE);	// Ray hits box
}



//========================================================
// Intersect an edge with an axis-aligned box
//========================================================
EBool E3d_EdgeIntersectAABox(E3dCoordinate* LBoxMin, E3dCoordinate* LBoxMax, E3dCoordinate* LV0, E3dCoordinate* LV1)
{
 E3dCoordinate		LX0=LV0[E3dX], LY0=LV0[E3dY], LZ0=LV0[E3dZ],
			LX1=LV1[E3dX], LY1=LV1[E3dY], LZ1=LV1[E3dZ],
			LD, LDX=LX1-LX0, LDY=LY1-LY0, LDZ=LZ1-LZ0,
			LT, LBBVar, LBBUMin, LBBVMin, LBBUMax, LBBVMax,
			LU, LV, LMin, LMax;

// Left YZ face
//
 E3dM_MINMAX(LX0, LX1, LMin, LMax);
 LD=LX1-LX0;
 LBBUMin=LBoxMin[E3dY];LBBUMax=LBoxMax[E3dY];
 LBBVMin=LBoxMin[E3dZ];LBBVMax=LBoxMax[E3dZ];

 LBBVar=LBoxMin[E3dX];
 if( (LMin <= LBBVar) && (LMax >= LBBVar) )
 {
  LT=(LBBVar-LX0)/LD;
  LU=LY0+LDY*LT;
  LV=LZ0+LDZ*LT;
  if((LU>=LBBUMin)&&(LU<LBBUMax)&&(LV>=LBBVMin)&&(LV<LBBVMax)) return(TRUE);
 }

// Right YZ face
//
 LBBVar=LBoxMax[E3dX];
 if( (LMin < LBBVar) && (LMax >= LBBVar) )
 {
  LT=(LBBVar-LX0)/LD;
  LU=LY0+LDY*LT;
  LV=LZ0+LDZ*LT;
  if((LU>=LBBUMin)&&(LU<LBBUMax)&&(LV>=LBBVMin)&&(LV<LBBVMax)) return(TRUE);
 }


// Top XZ face
//
 E3dM_MINMAX(LY0, LY1, LMin, LMax);
 LD=LY1-LY0;
 LBBUMin=LBoxMin[E3dX];LBBUMax=LBoxMax[E3dX];
 LBBVMin=LBoxMin[E3dZ];LBBVMax=LBoxMax[E3dZ];

 LBBVar=LBoxMin[E3dY];
 if( (LMin <= LBBVar) && (LMax >= LBBVar) )
 {
  LT=(LBBVar-LY0)/LD;
  LU=LX0+LDX*LT;
  LV=LZ0+LDZ*LT;
  if((LU>=LBBUMin)&&(LU<LBBUMax)&&(LV>=LBBVMin)&&(LV<LBBVMax)) return(TRUE);
 }

// Bottom XZ face
//
 LBBVar=LBoxMax[E3dY];
 if( (LMin < LBBVar) && (LMax >= LBBVar) )
 {
  LT=(LBBVar-LY0)/LD;
  LU=LX0+LDX*LT;
  LV=LZ0+LDZ*LT;
  if((LU>=LBBUMin)&&(LU<LBBUMax)&&(LV>=LBBVMin)&&(LV<LBBVMax)) return(TRUE);
 }


// Front XY face
//
 E3dM_MINMAX(LZ0, LZ1, LMin, LMax);
 LD=LZ1-LZ0;
 LBBUMin=LBoxMin[E3dX];LBBUMax=LBoxMax[E3dX];
 LBBVMin=LBoxMin[E3dY];LBBVMax=LBoxMax[E3dY];

 LBBVar=LBoxMin[E3dZ];
 if( (LMin <= LBBVar) && (LMax >= LBBVar) )
 {
  LT=(LBBVar-LZ0)/LD;
  LU=LX0+LDX*LT;
  LV=LY0+LDY*LT;
  if((LU>=LBBUMin)&&(LU<LBBUMax)&&(LV>=LBBVMin)&&(LV<LBBVMax)) return(TRUE);
 }

// Back XY face
//
 LBBVar=LBoxMax[E3dZ];
 if( (LMin < LBBVar) && (LMax >= LBBVar) )
 {
  LT=(LBBVar-LZ0)/LD;
  LU=LX0+LDX*LT;
  LV=LY0+LDY*LT;
  if((LU>=LBBUMin)&&(LU<LBBUMax)&&(LV>=LBBVMin)&&(LV<LBBVMax)) return(TRUE);
 }



 return(FALSE);
}


// Axis-aligned ray directions
//
E3d3DPosition	_PXDir= {  1.0,  0.0,  0.0 },
		_PYDir= {  0.0,  1.0,  0.0 },
		_PZDir= {  0.0,  0.0,  1.0 },
		_MXDir= { -1.0,  0.0,  0.0 },
		_MYDir= {  0.0, -1.0,  0.0 },
		_MZDir= {  0.0,  0.0, -1.0 };


//========================================================
// Test if the given Triangle is the given Voxel
//========================================================
EBool E3d_TriangleInVoxel(E3dScene* LScene, E3dRTriangle* LRTriangle, int LX, int LY, int LZ)
{
 E3dCoordinate*	LV0=LRTriangle->V0;
 E3dCoordinate*	LV1=LRTriangle->V1;
 E3dCoordinate*	LV2=LRTriangle->V2;
 E3dCoordinate*	LEdge1=LRTriangle->Edge1;
 E3dCoordinate*	LEdge2=LRTriangle->Edge2;
 E3d3DPosition	LOrigin;
 E3dCoordinate	LVoxelMin[3], LVoxelMax[3];
 E3dCoordinate	LGX, LGY, LGZ,
		LT;


 LGX=(E3dCoordinate)LX*LScene->VoxelXSize+LScene->BBoxMin.X;
 LGY=(E3dCoordinate)LY*LScene->VoxelYSize+LScene->BBoxMin.Y;
 LGZ=(E3dCoordinate)LZ*LScene->VoxelZSize+LScene->BBoxMin.Z;



// Check if any edges of the Voxel intersect the Triangle
//

// Originate 3 edges (rays) at a time, from the same corner of the Voxel
//
 LOrigin.X=LGX;
 LOrigin.Y=LGY;
 LOrigin.Z=LGZ;
 if(E3d_RayIntersectWithTriangleEO(&LOrigin, &_PXDir, LV0, LEdge1, LEdge2, &LT))
 {
  if(LT<LScene->VoxelXSize) return(TRUE);
 }
 if(E3d_RayIntersectWithTriangleEO(&LOrigin, &_PYDir, LV0, LEdge1, LEdge2, &LT))
 {
  if(LT<LScene->VoxelYSize) return(TRUE);
 }
 if(E3d_RayIntersectWithTriangleEO(&LOrigin, &_PZDir, LV0, LEdge1, LEdge2, &LT))
 {
  if(LT<LScene->VoxelZSize) return(TRUE);
 }

 LOrigin.X=LGX;
 LOrigin.Y=LGY+LScene->VoxelYSize;
 LOrigin.Z=LGZ+LScene->VoxelZSize;
 if(E3d_RayIntersectWithTriangleEO(&LOrigin, &_PXDir, LV0, LEdge1, LEdge2, &LT))
 {
  if(LT<LScene->VoxelXSize) return(TRUE);
 }
 if(E3d_RayIntersectWithTriangleEO(&LOrigin, &_MYDir, LV0, LEdge1, LEdge2, &LT))
 {
  if(LT<LScene->VoxelYSize) return(TRUE);
 }
 if(E3d_RayIntersectWithTriangleEO(&LOrigin, &_MZDir, LV0, LEdge1, LEdge2, &LT))
 {
  if(LT<LScene->VoxelZSize) return(TRUE);
 }

 LOrigin.X=LGX+LScene->VoxelXSize;
 LOrigin.Y=LGY;
 LOrigin.Z=LGZ+LScene->VoxelZSize;
 if(E3d_RayIntersectWithTriangleEO(&LOrigin, &_MXDir, LV0, LEdge1, LEdge2, &LT))
 {
  if(LT<LScene->VoxelXSize) return(TRUE);
 }
 if(E3d_RayIntersectWithTriangleEO(&LOrigin, &_PYDir, LV0, LEdge1, LEdge2, &LT))
 {
  if(LT<LScene->VoxelYSize) return(TRUE);
 }
 if(E3d_RayIntersectWithTriangleEO(&LOrigin, &_MZDir, LV0, LEdge1, LEdge2, &LT))
 {
  if(LT<LScene->VoxelZSize) return(TRUE);
 }

 LOrigin.X=LGX+LScene->VoxelXSize;
 LOrigin.Y=LGY+LScene->VoxelYSize;
 LOrigin.Z=LGZ;
 if(E3d_RayIntersectWithTriangleEO(&LOrigin, &_MXDir, LV0, LEdge1, LEdge2, &LT))
 {
  if(LT<LScene->VoxelXSize) return(TRUE);
 }
 if(E3d_RayIntersectWithTriangleEO(&LOrigin, &_MYDir, LV0, LEdge1, LEdge2, &LT))
 {
  if(LT<LScene->VoxelYSize) return(TRUE);
 }
 if(E3d_RayIntersectWithTriangleEO(&LOrigin, &_PZDir, LV0, LEdge1, LEdge2, &LT))
 {
  if(LT<LScene->VoxelZSize) return(TRUE);
 }





// Now check if any of the edges of the Triangle intersect a face of the Voxel
//
 LVoxelMin[E3dX]=LGX;
 LVoxelMin[E3dY]=LGY;
 LVoxelMin[E3dZ]=LGZ;
 LVoxelMax[E3dX]=LGX+LScene->VoxelXSize;
 LVoxelMax[E3dY]=LGY+LScene->VoxelYSize;
 LVoxelMax[E3dZ]=LGZ+LScene->VoxelZSize;


// V0-V1 edge
//
 if(E3d_EdgeIntersectAABox(LVoxelMin, LVoxelMax, LV0, LV1)) return(TRUE);

// V0-V2 edge
//
 if(E3d_EdgeIntersectAABox(LVoxelMin, LVoxelMax, LV0, LV2)) return(TRUE);

// V1-V2 edge
//
 if(E3d_EdgeIntersectAABox(LVoxelMin, LVoxelMax, LV1, LV2)) return(TRUE);


 return(FALSE);
}


//========================================================
// Add an E3dRTriangle to the Voxels it occupies
//========================================================
void E3d_GridAddTriangle(E3dScene* LScene, E3dRTriangle* LRTriangle, E3dPolyGroup* LPolyGroup)
{
 E3dVoxel*	LVoxel;
 E3dCoordinate*	LV0=LRTriangle->V0;	// The Triangle's Vertices
 E3dCoordinate*	LV1=LRTriangle->V1;
 E3dCoordinate*	LV2=LRTriangle->V2;
 E3dCoordinate	LBBMinX=LScene->BBoxMin.X, LBBMinY=LScene->BBoxMin.Y, LBBMinZ=LScene->BBoxMin.Z;

 int		LPageSize = LScene->GridXCount*LScene->GridYCount;
 int		LRowSize = LScene->GridXCount;
 int		LXMin, LYMin, LZMin,
		LXMax, LYMax, LZMax,
		LX, LY, LZ,
		LVoxelIndex, LV0VoxelIndex, LV1VoxelIndex, LV2VoxelIndex;


 LXMin=LXMax=(int)((LV0[E3dX]-LBBMinX)*LScene->VoxelXSizeInv);
 LYMin=LYMax=(int)((LV0[E3dY]-LBBMinY)*LScene->VoxelYSizeInv);
 LZMin=LZMax=(int)((LV0[E3dZ]-LBBMinZ)*LScene->VoxelZSizeInv);
 LV0VoxelIndex=LZMin*LPageSize + LYMin*LRowSize + LXMin;

 LX=(int)((LV1[E3dX]-LBBMinX)*LScene->VoxelXSizeInv);
 LY=(int)((LV1[E3dY]-LBBMinY)*LScene->VoxelYSizeInv);
 LZ=(int)((LV1[E3dZ]-LBBMinZ)*LScene->VoxelZSizeInv);
 LV1VoxelIndex=LZ*LPageSize + LY*LRowSize + LX;
 if(LX<LXMin) LXMin=LX;else if(LX>LXMax) LXMax=LX;
 if(LY<LYMin) LYMin=LY;else if(LY>LYMax) LYMax=LY;
 if(LZ<LZMin) LZMin=LZ;else if(LZ>LZMax) LZMax=LZ;

 LX=(int)((LV2[E3dX]-LBBMinX)*LScene->VoxelXSizeInv);
 LY=(int)((LV2[E3dY]-LBBMinY)*LScene->VoxelYSizeInv);
 LZ=(int)((LV2[E3dZ]-LBBMinZ)*LScene->VoxelZSizeInv);
 LV2VoxelIndex=LZ*LPageSize + LY*LRowSize + LX;
 if(LX<LXMin) LXMin=LX;else if(LX>LXMax) LXMax=LX;
 if(LY<LYMin) LYMin=LY;else if(LY>LYMax) LYMax=LY;
 if(LZ<LZMin) LZMin=LZ;else if(LZ>LZMax) LZMax=LZ;


 for(LZ=LZMin;LZ<=LZMax;LZ++)
 {
  for(LY=LYMin;LY<=LYMax;LY++)
  {
   for(LX=LXMin;LX<=LXMax;LX++)
   {
    LVoxelIndex=LZ*LPageSize + LY*LRowSize + LX;
    if((LVoxelIndex==LV0VoxelIndex) || (LVoxelIndex==LV1VoxelIndex) || (LVoxelIndex==LV2VoxelIndex) || E3d_TriangleInVoxel(LScene, LRTriangle, LX, LY, LZ))
    {
     LVoxel=LScene->Voxels + LVoxelIndex;
     ELst_AddPointerChk((void***)(&(LVoxel->Triangles)), &(LVoxel->NumOfTriangles), LRTriangle);
/*
     if(ELst_AddPointerChk((void***)(&(LVoxel->Triangles)), &(LVoxel->NumOfTriangles), LRTriangle))
     {
      if(LVoxel->PolyGroups==NULL) LVoxel->PolyGroups=(E3dPolyGroup**)EMalloc(sizeof(E3dPolyGroup*));
      else LVoxel->PolyGroups=(E3dPolyGroup**)ERealloc(LVoxel->PolyGroups, sizeof(E3dPolyGroup*)*LVoxel->NumOfTriangles);
      LVoxel->PolyGroups[LVoxel->NumOfTriangles-1]=LPolyGroup;
     }
*/
    }


   }
  }
 }
}


//========================================================================================
// Subdivide Scene into a uniform grid
// Must call E3d_GetRenderTriangles() first, to fill out LScene->RenderPolyGroups!
//========================================================================================
E3dVoxel* E3d_SceneVoxelize(E3dScene* LScene)
{
 E3dVoxel*	LVoxels;
 E3dRenderInfo*	LRenderInfo=&(LScene->RenderInfo);
 unsigned int	LGridXCount=LRenderInfo->GridXCount, LGridYCount=LRenderInfo->GridYCount, LGridZCount=LRenderInfo->GridZCount;


 if(LScene->Voxels) E3d_SceneFreeVoxels(LScene);

 if(E3d_SceneGetRenderBoundingBox(LScene));
 else return(NULL);

 LVoxels=(E3dVoxel*)EMalloc(sizeof(E3dVoxel)*LGridXCount*LGridYCount*LGridZCount);
 if(LVoxels)
 {
  E3dPolyGroup**	LRenderPolyGroups;
  E3dPolyGroup*		LRenderPolyGroup;
  E3dRTriangle*		LRTriangle;
  E3dCoordinate		LXSpan=LScene->BBoxMax.X-LScene->BBoxMin.X,
			LYSpan=LScene->BBoxMax.Y-LScene->BBoxMin.Y,
			LZSpan=LScene->BBoxMax.Z-LScene->BBoxMin.Z;
  long			LStartTime;
  unsigned int		LC, LN, LTC, LTN;


  LScene->Voxels=LVoxels;

// We have GridAllocatedCount to avoid the application changing the grid resolution while a grid is allocated.
// This assures that E3d_SceneFreeVoxels() won't free the wrong number of Voxels.
//
  LScene->GridXCount=LGridXCount;
  LScene->GridYCount=LGridYCount;
  LScene->GridZCount=LGridZCount;

// Initialize global static variables for Voxel traversal
//
  LScene->RayCounter=0;

  LScene->VoxelXSize=LXSpan/(E3dCoordinate)LGridXCount;
  LScene->VoxelYSize=LYSpan/(E3dCoordinate)LGridYCount;
  LScene->VoxelZSize=LZSpan/(E3dCoordinate)LGridZCount;

  LScene->VoxelXSizeInv=1.0 / LScene->VoxelXSize;
  LScene->VoxelYSizeInv=1.0 / LScene->VoxelYSize;
  LScene->VoxelZSizeInv=1.0 / LScene->VoxelZSize;


// Maximum distance: the largest of:
// - the diagonal of the Scene BBox
// - the eye-to-far-corner distance
//
// FIXME
//
//  _MaxDistance = sqrt(LXSpan*LXSpan + LYSpan*LYSpan + LZSpan*LZSpan);
//  _MaxDistance = HUGE;


// Clear the grid
//
  memset(LVoxels, 0, sizeof(E3dVoxel)*LGridXCount*LGridYCount*LGridZCount);

  LStartTime=E3d_GetSeconds();

printf("Loading Geometries into grid");fflush(stdout);

// Go through the Triangles and put them in the appropriate voxels
//
  LN=LScene->NumOfRenderPolyGroups;
  LRenderPolyGroups=LScene->RenderPolyGroups;
  for(LC=0;LC<LN;LC++)
  {
   LRenderPolyGroup=LRenderPolyGroups[LC];

   LRTriangle=LRenderPolyGroup->Triangles;
   LTN=LRenderPolyGroup->NumOfTriangles;

   for(LTC=0;LTC<LTN;LTC++, LRTriangle++)
   {
    E3d_GridAddTriangle(LScene, LRTriangle, LRenderPolyGroup);
   }
  }
printf(" OK %d second(s)\n", (int)(E3d_GetSeconds()-LStartTime));fflush(stdout);
 }
 return(LScene->Voxels);
}


//================================================
// Free the uniform-grid of a Scene
//================================================
EBool E3d_SceneFreeVoxels(E3dScene* LScene)
{
 E3dVoxel*	LVoxels=LScene->Voxels;

 if(LVoxels)
 {
  E3dVoxel*	LVoxel=LVoxels;
  unsigned int	LC, LN;


  LN=LScene->GridXCount*LScene->GridYCount*LScene->GridZCount;
  for(LC=0;LC<LN;LC++, LVoxel++)
  {
   if(LVoxel->Triangles) EFree(LVoxel->Triangles);
//   if(LVoxel->PolyGroups) EFree(LVoxel->PolyGroups);
  }
  EFree(LVoxels);
  LScene->Voxels=NULL;

  return(TRUE);
 }

 return(FALSE);
}


//================================================================
// Initialize a RenderInfo structure
//================================================================
void E3d_RenderInfoDefault(E3dRenderInfo* LRenderInfo)
{
 LRenderInfo->Camera=NULL;

 LRenderInfo->Cameras=NULL;
 LRenderInfo->NumOfCameras=0;

 LRenderInfo->Renderer=NULL;

// Image
//
 LRenderInfo->ImageXSize=320;
 LRenderInfo->ImageYSize=240;
 LRenderInfo->WorkingImageXSize=0;
 LRenderInfo->WorkingImageYSize=0;

 LRenderInfo->WindowXSize=320;
 LRenderInfo->WindowYSize=240;
 LRenderInfo->StartPixelSize=8;

 LRenderInfo->AntiAliasing=E3dAA_NONE;
 LRenderInfo->CreateHDRImage=FALSE;

 LRenderInfo->Brightness=1.0;
 LRenderInfo->Gamma=2.2;


 LRenderInfo->DoubleBuffer=FALSE;
 LRenderInfo->Interlace=FALSE;
 LRenderInfo->ClearFrameBuffer=FALSE;
 LRenderInfo->ClearZBuffer=FALSE;

 Ec_SetRGBAiColor(&(LRenderInfo->BackgroundRGBAiColor), 0x0000, 0x0000, 0x0000, 0x0000);
 LRenderInfo->BackfaceCulling=TRUE;

 LRenderInfo->AlphaBlending=FALSE;

 LRenderInfo->Algorithm=E3dRENDER_RAYTRACING;

// Ray-tracer
//
 LRenderInfo->Shadows=TRUE;
 LRenderInfo->Reflections=TRUE;
 LRenderInfo->Transparency=TRUE;
 LRenderInfo->Refractions=TRUE;
 LRenderInfo->MaxRayDepth=7;

 LRenderInfo->GridXCount=100;
 LRenderInfo->GridYCount=100;
 LRenderInfo->GridZCount=100;



 LRenderInfo->FogColor.R=0.5;LRenderInfo->FogColor.G=0.5;LRenderInfo->FogColor.B=0.5;
 LRenderInfo->SceneAmbientColor.R=E3dDEF_AMBIENCE;LRenderInfo->SceneAmbientColor.G=E3dDEF_AMBIENCE;LRenderInfo->SceneAmbientColor.B=E3dDEF_AMBIENCE;
 LRenderInfo->FogNear=0.1;LRenderInfo->FogFar=100.0;
 LRenderInfo->Fog=FALSE;
 LRenderInfo->Textures=TRUE;

 LRenderInfo->SizeFactor=0.0;

 LRenderInfo->CanvasImage=NULL;

 LRenderInfo->UpdateFlags=E3dUF_ALL;
 LRenderInfo->LinesDone=0;
 LRenderInfo->CurrentPixelSize=0;


 LRenderInfo->Rendering=FALSE;
 LRenderInfo->GoRender=FALSE;
 LRenderInfo->FirstIteration=FALSE;
 LRenderInfo->RestartRender=FALSE;
}


#define M_CopyField(mField)	LRenderInfo->mField=LSRenderInfo->mField

//================================================================
// Copy the changeable fields of a RenderInfo structure
//================================================================
void E3d_RenderInfoCopy(E3dRenderInfo* LSRenderInfo, E3dRenderInfo* LRenderInfo)
{
 M_CopyField(ImageXSize);
 M_CopyField(ImageYSize);
 M_CopyField(WindowXSize);
 M_CopyField(WindowYSize);
 M_CopyField(StartPixelSize);
 M_CopyField(AntiAliasing);
 M_CopyField(CreateHDRImage);
 M_CopyField(SizeFactor);
 M_CopyField(Brightness);
 M_CopyField(Gamma);
 M_CopyField(MaxRayDepth);
 M_CopyField(FogNear);
 M_CopyField(FogFar);
 M_CopyField(BackgroundRGBAiColor);
 M_CopyField(Algorithm);
 M_CopyField(GridXCount);
 M_CopyField(GridYCount);
 M_CopyField(GridZCount);
 M_CopyField(DoubleBuffer);
 M_CopyField(Interlace);
 M_CopyField(ClearFrameBuffer);
 M_CopyField(ClearZBuffer);
 M_CopyField(BackfaceCulling);
 M_CopyField(AlphaBlending);
 M_CopyField(Textures);
 M_CopyField(Shadows);
 M_CopyField(Reflections);
 M_CopyField(Transparency);
 M_CopyField(Refractions);
 M_CopyField(Fog);
}


//================================================================
// Clear Ray numbers on the renderable Triangles of a Scene
//================================================================
void E3d_SceneResetRayCounters(E3dScene* LScene)
{
 if(LScene->RenderPolyGroups)
 {
  E3dPolyGroup**	LRenderPolyGroups=LScene->RenderPolyGroups;
  E3dPolyGroup*		LRenderPolyGroup;
  E3dRTriangle*		LRTriangle;
  unsigned int		LC, LTC, LTN, LNumOfRenderPolyGroups=LScene->NumOfRenderPolyGroups;


  for(LC=0;LC<LNumOfRenderPolyGroups;LC++, LRenderPolyGroup++)
  {
   LRenderPolyGroup=LRenderPolyGroups[LC];

   LTN=LRenderPolyGroup->NumOfTriangles;LRTriangle=LRenderPolyGroup->Triangles;
   for(LTC=0;LTC<LTN;LTC++, LRTriangle++)
   {
    LRTriangle->LastRayNumber=-1;
   }
  }
 }
}


void E3d_SceneResetForRendering(E3dScene* LScene)
{
 LScene->ShadowRayLastHit=NULL;
 LScene->TestRayLastHit=NULL;
 LScene->RayCounter=0;
 E3d_SceneResetRayCounters(LScene);
}


//========================================================================================================
// Voxel stepping 3D-DDA for uniform-grid spatial subdivision
//========================================================================================================

#define EPSILON	0.000001


/*
#define x2voxel(px) 	((LScene->VoxelXSize<EPSILON) ? 0 : ((px) - LScene->BBoxMin.X) * LScene->VoxelXSizeInv)
#define y2voxel(py) 	((LScene->VoxelYSize<EPSILON) ? 0 : ((py) - LScene->BBoxMin.Y) * LScene->VoxelYSizeInv)
#define z2voxel(pz) 	((LScene->VoxelZSize<EPSILON) ? 0 : ((pz) - LScene->BBoxMin.Z) * LScene->VoxelZSizeInv)
*/

#define x2voxel(px) 	( ((px) - LScene->BBoxMin.X) * LScene->VoxelXSizeInv )
#define y2voxel(py) 	( ((py) - LScene->BBoxMin.Y) * LScene->VoxelYSizeInv )
#define z2voxel(pz) 	( ((pz) - LScene->BBoxMin.Z) * LScene->VoxelZSizeInv )

#define VoxelToX(px)	((px) * LScene->VoxelXSize + LScene->BBoxMin.X)
#define VoxelToY(py)	((py) * LScene->VoxelYSize + LScene->BBoxMin.Y)
#define VoxelToZ(pz)	((pz) * LScene->VoxelZSize + LScene->BBoxMin.Z)




int E3d_BoundsIntersectingSegment(E3dScene* LScene, E3dRay* LRay, float* tmin, float* tmax)
{
 float	t, mindist, maxdist;
 float	dir, pos;

 mindist = *tmin;
 maxdist = *tmax;

 pos = LRay->X;
 dir = LRay->DirX;

 if(dir < 0)
 {
   t = (LScene->BBoxMin.X - pos) / dir;
   if(t < *tmin) return(FALSE);
   if(t <= *tmax) *tmax = t;
   t = (LScene->BBoxMax.X - pos) / dir;
   if(t >= *tmin)
   {
     if(t > *tmax * (1.+EPSILON)) return(FALSE);
     *tmin = t;
   }
 }
 else if(dir > 0.)
      {
       t = (LScene->BBoxMax.X - pos) / dir;
       if(t < *tmin) return(FALSE);
       if(t <= *tmax) *tmax = t;
       t = (LScene->BBoxMin.X - pos) / dir;
       if(t >= *tmin)
       {
        if(t > *tmax * (1.+EPSILON)) return(FALSE);
        *tmin = t;
       }
      } else if(pos < LScene->BBoxMin.X || pos > LScene->BBoxMax.X) return(FALSE);


 pos = LRay->Y;
 dir = LRay->DirY;

 if(dir < 0)
 {
   t = (LScene->BBoxMin.Y - pos) / dir;
   if(t < *tmin) return FALSE;
   if(t <= *tmax) *tmax = t;
   t = (LScene->BBoxMax.Y - pos) / dir;
   if(t >= *tmin)
   {
    if(t > *tmax * (1.+EPSILON)) return FALSE;
    *tmin = t;
   }
 }
 else if(dir > 0.)
      {
       t = (LScene->BBoxMax.Y - pos) / dir;
       if(t < *tmin) return(FALSE);
       if(t <= *tmax) *tmax = t;
       t = (LScene->BBoxMin.Y - pos) / dir;
       if(t >= *tmin)
       {
        if(t > *tmax * (1.+EPSILON)) return FALSE;
        *tmin = t;
       }
      } else if(pos < LScene->BBoxMin.Y || pos > LScene->BBoxMax.Y) return FALSE;
 

 pos = LRay->Z;
 dir = LRay->DirZ;
 
 if(dir < 0)
 {
   t = (LScene->BBoxMin.Z - pos) / dir;
   if(t < *tmin) return FALSE;
   if(t <= *tmax) *tmax = t;
   t = (LScene->BBoxMax.Z - pos) / dir; 
   if(t >= *tmin)
   {
     if(t > *tmax * (1.+EPSILON)) return FALSE;
     *tmin = t;
   }
 }
 else if(dir > 0.)
      {
       t = (LScene->BBoxMax.Z - pos) / dir;
       if(t < *tmin) return FALSE;
       if(t <= *tmax) *tmax = t;
       t = (LScene->BBoxMin.Z - pos) / dir;
       if(t >= *tmin)
       {
        if(t > *tmax * (1.+EPSILON)) return FALSE;
        *tmin = t;
       }
      }
      else if(pos < LScene->BBoxMin.Z || pos > LScene->BBoxMax.Z) return FALSE;


// If *tmin == mindist, then there was no "near"
// intersection farther than EPSILON away.
//
 if(*tmin == mindist)
 {
  if(*tmax < maxdist) return(TRUE);
 }
 else
 {
  if(*tmin < maxdist) return(TRUE);
 }
 return(FALSE);		// Hit, but not closer than maxdist
}


//========================================================================
// Compute t0, ray's minimal intersection with the whole grid and
// position P of this intersection. Returns TRUE if the grid bounds are
// intersected and FALSE if the ray passes along the voxel grid.
//========================================================================
int E3d_RayGridBoundsIntersect(E3dScene* LScene, E3dRay* LRay)
{
 E3dCoordinate	LMinDist=0.0, LMaxDist=HUGE;
 E3dCoordinate	LTMin=LMinDist, LTMax=LMaxDist;
 int		LHit = E3d_BoundsIntersectingSegment(LScene, LRay, &LTMin, &LTMax);

 if(LHit)
 {
// Reduce maxdist
//
  if(LTMin == LMinDist)
  {
   if(LTMax < LMaxDist) LMaxDist = LTMax;
  }
  else
  {
   if(LTMin < LMaxDist) LMaxDist = LTMin;
  }

  LRay->X = LRay->X + LTMin * LRay->DirX;
  LRay->Y = LRay->Y + LTMin * LRay->DirY;
  LRay->Z = LRay->Z + LTMin * LRay->DirZ;
  LRay->LastT=LTMin;
 }
 else return(FALSE);

 return(TRUE);
}


//================================================
// Initialize Voxel-traversal for a Ray
//================================================
void E3d_RayDDAInitStartVoxel(E3dScene* LScene, E3dRay* LRay)
{
// Compute the Voxel where this intersection occurs
//
 LRay->GridX = x2voxel(LRay->X); if(LRay->GridX >= LScene->GridXCount) LRay->GridX = LScene->GridXCount-1;
 LRay->GridY = y2voxel(LRay->Y); if(LRay->GridY >= LScene->GridYCount) LRay->GridY = LScene->GridYCount-1;
 LRay->GridZ = z2voxel(LRay->Z); if(LRay->GridZ >= LScene->GridZCount) LRay->GridZ = LScene->GridZCount-1;
}



//================================================
// Initialize Voxel-traversal for a Ray
//================================================
void E3d_RayDDAInit(E3dScene* LScene, E3dRay* LRay)
{
 E3dCoordinate	LInv;


// Compute the Voxel where this intersection occurs
//
 LRay->GridX = x2voxel(LRay->X); if(LRay->GridX >= LScene->GridXCount) LRay->GridX = LScene->GridXCount-1;
 LRay->GridY = y2voxel(LRay->Y); if(LRay->GridY >= LScene->GridYCount) LRay->GridY = LScene->GridYCount-1;
 LRay->GridZ = z2voxel(LRay->Z); if(LRay->GridZ >= LScene->GridZCount) LRay->GridZ = LScene->GridZCount-1;


// LRay->GridDeltaX   - distance increment along the ray to the adjacent Voxel in the X direction.
// LRay->GridNextX    - the total distance from the ray origin to the next Voxel in the X direction.
// LRay->GridStepX    - either +1 or -1 accroding to the ray X direction. 
// LRay->GridOutX     - -1 or _GridXCount: the first X Voxel index outside the grid
//

// X
//
 if(fabs(LRay->DirX) > EPSILON)
 {
   LInv=1.0 / LRay->DirX;
   if(LRay->DirX > 0.0)
   {
     LRay->GridDeltaX = LScene->VoxelXSize * LInv;
     LRay->GridNextX = (VoxelToX(LRay->GridX+1) - LRay->X) * LInv;
     LRay->GridStepX = 1; LRay->GridOutX = LScene->GridXCount;
   }
   else
   {
     LRay->GridDeltaX = LScene->VoxelXSize * -LInv;
     LRay->GridNextX = (VoxelToX(LRay->GridX) - LRay->X) * LInv;
     LRay->GridStepX = LRay->GridOutX = -1;
   }
 }
 else
 {
   LRay->GridDeltaX = 0.0;
   LRay->GridNextX = HUGE;
 }

// Y
//
 if(fabs(LRay->DirY) > EPSILON)
 {
   LInv=1.0 / LRay->DirY;
   if(LRay->DirY > 0.0)
   {
     LRay->GridDeltaY = LScene->VoxelYSize * LInv;
     LRay->GridNextY = (VoxelToY(LRay->GridY+1) - LRay->Y) * LInv;
     LRay->GridStepY = 1; LRay->GridOutY = LScene->GridYCount;
   }
   else
   {
     LRay->GridDeltaY = LScene->VoxelYSize * -LInv;
     LRay->GridNextY = (VoxelToY(LRay->GridY) - LRay->Y) * LInv;
     LRay->GridStepY = LRay->GridOutY = -1;
   }
 }
 else
 {
   LRay->GridDeltaY = 0.0;
   LRay->GridNextY = HUGE;
 }

// Z
//
 if(fabs(LRay->DirZ) > EPSILON)
 {
   LInv=1.0 / LRay->DirZ;
   if(LRay->DirZ > 0.0)
   {
     LRay->GridDeltaZ = LScene->VoxelZSize * LInv;
     LRay->GridNextZ = (VoxelToZ(LRay->GridZ+1) - LRay->Z) * LInv;
     LRay->GridStepZ = 1; LRay->GridOutZ = LScene->GridZCount;
   }
   else
   {
     LRay->GridDeltaZ = LScene->VoxelZSize * -LInv;
     LRay->GridNextZ = (VoxelToZ(LRay->GridZ) - LRay->Z) * LInv;
     LRay->GridStepZ = LRay->GridOutZ = -1;
   }
 }
 else
 {
   LRay->GridDeltaZ = 0.0;
   LRay->GridNextZ = HUGE;
 }
}


//================================================================
// Initialize Voxel-traversal for a Ray
// This function assumes that the starting Voxel is already set
//================================================================
void E3d_RayDDAInitStartSet(E3dScene* LScene, E3dRay* LRay)
{
 E3dCoordinate	LInv;


// LRay->GridDeltaX   - distance increment along the ray to the adjacent Voxel in the X direction.
// LRay->GridNextX    - the total distance from the ray origin to the next Voxel in the X direction.
// LRay->GridStepX    - either +1 or -1 accroding to the ray X direction. 
// LRay->GridOutX     - -1 or _GridXCount: the first X Voxel index outside the grid
//

// X
//
 if(fabs(LRay->DirX) > EPSILON)
 {
   LInv=1.0 / LRay->DirX;
   if(LRay->DirX > 0.0)
   {
     LRay->GridDeltaX = LScene->VoxelXSize * LInv;
     LRay->GridNextX = (VoxelToX(LRay->GridX+1) - LRay->X) * LInv;
     LRay->GridStepX = 1; LRay->GridOutX = LScene->GridXCount;
   }
   else
   {
     LRay->GridDeltaX = LScene->VoxelXSize * -LInv;
     LRay->GridNextX = (VoxelToX(LRay->GridX) - LRay->X) * LInv;
     LRay->GridStepX = LRay->GridOutX = -1;
   }
 }
 else
 {
   LRay->GridDeltaX = 0.0;
   LRay->GridNextX = HUGE;
 }

// Y
//
 if(fabs(LRay->DirY) > EPSILON)
 {
   LInv=1.0 / LRay->DirY;
   if(LRay->DirY > 0.0)
   {
     LRay->GridDeltaY = LScene->VoxelYSize * LInv;
     LRay->GridNextY = (VoxelToY(LRay->GridY+1) - LRay->Y) * LInv;
     LRay->GridStepY = 1; LRay->GridOutY = LScene->GridYCount;
   }
   else
   {
     LRay->GridDeltaY = LScene->VoxelYSize * -LInv;
     LRay->GridNextY = (VoxelToY(LRay->GridY) - LRay->Y) * LInv;
     LRay->GridStepY = LRay->GridOutY = -1;
   }
 }
 else
 {
   LRay->GridDeltaY = 0.0;
   LRay->GridNextY = HUGE;
 }

// Z
//
 if(fabs(LRay->DirZ) > EPSILON)
 {
   LInv=1.0 / LRay->DirZ;
   if(LRay->DirZ > 0.0)
   {
     LRay->GridDeltaZ = LScene->VoxelZSize * LInv;
     LRay->GridNextZ = (VoxelToZ(LRay->GridZ+1) - LRay->Z) * LInv;
     LRay->GridStepZ = 1; LRay->GridOutZ = LScene->GridZCount;
   }
   else
   {
     LRay->GridDeltaZ = LScene->VoxelZSize * -LInv;
     LRay->GridNextZ = (VoxelToZ(LRay->GridZ) - LRay->Z) * LInv;
     LRay->GridStepZ = LRay->GridOutZ = -1;
   }
 }
 else
 {
   LRay->GridDeltaZ = 0.0;
   LRay->GridNextZ = HUGE;
 }
}




#define E3dM_RayDDANextVoxel()\
 if(LRay->GridNextX <= LRay->GridNextY && LRay->GridNextX <= LRay->GridNextZ)\
 {							/* LRay->GridNextX is smallest	*/\
   LRay->GridX += LRay->GridStepX;\
   LRay->GridNextX += LRay->GridDeltaX;\
   LStillInGrid = LRay->GridX - LRay->GridOutX;		/* FALSE if LRay->GridX==LRay->GridOutX	*/\
 }\
 else if(LRay->GridNextY <= LRay->GridNextZ)\
      {							/* LRay->GridNextY is smallest	*/\
	LRay->GridY += LRay->GridStepY;\
	LRay->GridNextY += LRay->GridDeltaY;\
	LStillInGrid = LRay->GridY - LRay->GridOutY;\
      }\
      else\
      {							/* LRay->GridNextZ is smallest	*/\
	LRay->GridZ += LRay->GridStepZ;\
	LRay->GridNextZ += LRay->GridDeltaZ;\
	LStillInGrid = LRay->GridZ - LRay->GridOutZ;\
      }



//========================================================================================================
// Advance ray to the next Voxel. Needs setup with E3d_RayDDAInit().
//
// Return value
//  FALSE (0) if the current Voxel was the last in the grid intersected by the ray, otherwise non-zero
//========================================================================================================
static int E3d_RayDDANextVoxel(E3dRay* LRay)
{
 int	LStillInGrid;


 if(LRay->GridNextX <= LRay->GridNextY && LRay->GridNextX <= LRay->GridNextZ)
 {							// LRay->GridNextX is smallest
   LRay->GridX += LRay->GridStepX;
   LRay->GridNextX += LRay->GridDeltaX;
   LStillInGrid = LRay->GridX - LRay->GridOutX;			// FALSE if LRay->GridX==LRay->GridOutX
 }
 else if(LRay->GridNextY <= LRay->GridNextZ)
      {							// LRay->GridNextY is smallest
	LRay->GridY += LRay->GridStepY;
	LRay->GridNextY += LRay->GridDeltaY;
	LStillInGrid = LRay->GridY - LRay->GridOutY;
      }
      else
      {							// LRay->GridNextZ is smallest
	LRay->GridZ += LRay->GridStepZ;
	LRay->GridNextZ += LRay->GridDeltaZ;
	LStillInGrid = LRay->GridZ - LRay->GridOutZ;
      }

 return(LStillInGrid);
}


//================================================================================
// Test a ray against objects in the Voxel-grid
// Only front-facing triangles will be accepted.
//================================================================================
E3dRTriangle* E3d_RayIntersectVoxelsBackFaceCull(E3dScene* LScene, E3dRay* LRay, E3dCoordinate* LUVRet, void* LExcludeFace)
{
 E3dVoxel*	LVoxel;
 E3dRTriangle*	LRTriangle;
 E3dRTriangle*	LClosestRTriangle=NULL;
 E3dCoordinate	LX, LY, LZ,
		LU, LV,
		LClosestT=0.0;
 int		LX1, LY1, LZ1,
		LVoxelIndex, LHitVoxelIndex=-1,
		LPageSize = LScene->GridXCount*LScene->GridYCount;
 int		LRowSize = LScene->GridXCount;
 unsigned int	LC, LN, LInVoxelHits=0;
// unsigned int	LHits=0;
// int		LClosestHitVoxelIndex=-1;


 LScene->RayCounter++;

 E3d_RayDDAInit(LScene, LRay);

 do
 {
// Check the Triangles in this Voxel
//
  LVoxelIndex=LRay->GridZ*LPageSize + LRay->GridY*LRowSize + LRay->GridX;
  LVoxel = LScene->Voxels + LVoxelIndex;

  LN=LVoxel->NumOfTriangles;
  if(LN)
  {
   for(LC=0;LC<LN;LC++)
   {
    LRTriangle=LVoxel->Triangles[LC];

    if(LRTriangle->Face!=LExcludeFace)
    {
/*
     if(LRTriangle->LastRayNumber < _RayCounter)
     {
*/

      if(E3d_RayIntersectWithTriangleBackCullE(LRay, LRTriangle->V0, LRTriangle->Edge1, LRTriangle->Edge2, &LU, &LV))
      {
       E3dCoordinate	LT=LRay->LastT;

// Make sure the hit was really in this Voxel
//
       LX=LRay->X+LT*LRay->DirX;
       LY=LRay->Y+LT*LRay->DirY;
       LZ=LRay->Z+LT*LRay->DirZ;
       LX1=(int)((LX-LScene->BBoxMin.X) * LScene->VoxelXSizeInv),
       LY1=(int)((LY-LScene->BBoxMin.Y) * LScene->VoxelYSizeInv),
       LZ1=(int)((LZ-LScene->BBoxMin.Z) * LScene->VoxelZSizeInv);

       LHitVoxelIndex=LZ1*LPageSize + LY1*LRowSize + LX1;

       if(LVoxelIndex==LHitVoxelIndex)
       {
	 if((LInVoxelHits==0)||(LRay->LastT<LClosestT))
	 {
	   LInVoxelHits++;
	   LRay->HitX=LX;
	   LRay->HitY=LY;
	   LRay->HitZ=LZ;
	   LUVRet[0]=LU;
	   LUVRet[1]=LV;
	   LClosestRTriangle=LRTriangle;
	   LClosestT=LRay->LastT;
	 }
       }
/*
       else
       {
	 if((LHits==0)||(LRay->LastT<LClosestT))
	 {
	   LRay->HitX=LX;
	   LRay->HitY=LY;
	   LRay->HitZ=LZ;
	   LUVRet[0]=LU;
	   LUVRet[1]=LV;
	   LClosestRTriangle=LRTriangle;
	   LClosestT=LRay->LastT;
	   LClosestHitVoxelIndex=LHitVoxelIndex;
	 }
       }

       LHits++;
*/
      }

// Mark this Triangle as already having been tested by this ray
//
/*
      LRTriangle->LastRayNumber = _RayCounter;
     }
     else
     {
// This Triangle was already tested against this Ray
//
      if((LClosestRTriangle==LRTriangle)&&(LInVoxelHits==0))
      {
       if(LVoxelIndex==LClosestHitVoxelIndex) LInVoxelHits++;
      }
     }

//    else { Printf("DoubleTestBFC\n");fflush(stdout); }
*/

    }

   }
  }
  if(LInVoxelHits)
  {
   LRay->LastT=LClosestT;
   return(LClosestRTriangle);
  }
 } while(E3d_RayDDANextVoxel(LRay));

 return(NULL);
}


//================================================================================
// Test a ray against objects in the uniform-grid
// Both sides of triangles are considered a hit
//================================================================================
E3dRTriangle* E3d_RayIntersectVoxels(E3dScene* LScene, E3dRay* LRay, E3dCoordinate* LUVRet, void* LExcludeFace)
{
 E3dVoxel*	LVoxel;
 E3dRTriangle**	LRTriangles;
 E3dRTriangle*	LRTriangle;
 E3dRTriangle*	LClosestRTriangle=NULL;
 E3dCoordinate	LX, LY, LZ, LU, LV,
		LClosestT=0.0;
 int		LX1, LY1, LZ1,
		LVoxelIndex, LHitVoxelIndex, LPageSize = LScene->GridXCount*LScene->GridYCount;
 int		LRowSize = LScene->GridXCount;
 unsigned int	LC, LN, LHits=0;
 EBool		LFrontFaceHit=TRUE;


 LScene->RayCounter++;

 E3d_RayDDAInit(LScene, LRay);

 do
 {
// Check the Triangles in this Voxel
//
  LVoxelIndex=LRay->GridZ*LPageSize + LRay->GridY*LRowSize + LRay->GridX;
  LVoxel = LScene->Voxels + LVoxelIndex;

  if((LN=LVoxel->NumOfTriangles)>0)
  {
   LRTriangles=LVoxel->Triangles;
   for(LC=0;LC<LN;LC++)
   {
    LRTriangle=LVoxel->Triangles[LC];

    if(LRTriangle->Face!=LExcludeFace)
    {
     if(E3d_RayIntersectWithTriangleFrontBackE(LRay, LRTriangle->V0, LRTriangle->Edge1, LRTriangle->Edge2, &LU, &LV))
     {
      {
// Make sure the hit was really in this Voxel
//
       LX=LRay->X+LRay->LastT*LRay->DirX;
       LY=LRay->Y+LRay->LastT*LRay->DirY;
       LZ=LRay->Z+LRay->LastT*LRay->DirZ;
       LX1=(int)((LX-LScene->BBoxMin.X) * LScene->VoxelXSizeInv),
       LY1=(int)((LY-LScene->BBoxMin.Y) * LScene->VoxelYSizeInv),
       LZ1=(int)((LZ-LScene->BBoxMin.Z) * LScene->VoxelZSizeInv);

       LHitVoxelIndex=LZ1*LPageSize + LY1*LRowSize + LX1;

       if(LVoxelIndex==LHitVoxelIndex)
       {
	if((LHits==0)||(LRay->LastT<LClosestT))
	{
	  LHits++;
	  LRay->HitX=LX;
	  LRay->HitY=LY;
	  LRay->HitZ=LZ;
	  LUVRet[0]=LU;
	  LUVRet[1]=LV;
	  LClosestRTriangle=LRTriangle;
	  LFrontFaceHit=LRay->FrontFaceHit;
	  LClosestT=LRay->LastT;
	}
       }
      }

     }
    }
   }
  }
  if(LHits)
  {
   LRay->LastT=LClosestT;
   LRay->FrontFaceHit=LFrontFaceHit;
//printf("Closest: %f\n", LClosestT);fflush(stdout);

   return(LClosestRTriangle);
  }
 } while(E3d_RayDDANextVoxel(LRay));

 return(NULL);
}


//================================================================================
// Test a ray against objects in the uniform-grid.
// This version tests whether there is a hit within the given distance.
//================================================================================
EBool E3d_ShadowRayIntersectVoxels(E3dScene* LScene, E3dRay* LRay, void* LExcludeFace, E3dCoordinate LDistance)
{
 E3dRTriangle*	LRTriangle;


 LScene->RayCounter++;


 if(LScene->ShadowRayLastHit)
 {
  LRTriangle=LScene->ShadowRayLastHit;

  if(LRTriangle->Face!=LExcludeFace)
  {
// Mark this Triangle as already having been tested by this ray
//
   LRTriangle->LastRayNumber = LScene->RayCounter;
//    if(E3d_RayIntersectWithTriangleE(LRay, LRTriangle->V0, LRTriangle->Edge1, LRTriangle->Edge2, &LU, &LV))
   if(E3d_RayIntersectWithTriangleBackCullENoUV(LRay, LRTriangle->V0, LRTriangle->Edge1, LRTriangle->Edge2))
   {
//printf("ShadowLastHit %p\n", LRTriangle);fflush(stdout);

    if(LRay->LastT<LDistance) return(TRUE);
   }
  }
 }


 if(LScene->Voxels)
 {
  E3dVoxel*		LVoxel;
  int			LVoxelIndex, LPageSize = LScene->GridXCount*LScene->GridYCount;
  int			LRowSize = LScene->GridXCount;
  unsigned int		LC;


  E3d_RayDDAInit(LScene, LRay);

  do
  {
// Check the Triangles in this Voxel
//
   LVoxelIndex=LRay->GridZ*LPageSize + LRay->GridY*LRowSize + LRay->GridX;
   LVoxel = LScene->Voxels + LVoxelIndex;


   for(LC=0;LC<LVoxel->NumOfTriangles;LC++)
   {
    LRTriangle=LVoxel->Triangles[LC];

    if(LRTriangle->Face!=LExcludeFace)
    {
// Check if this triangle has been tested by this ray
//
     if(LRTriangle->LastRayNumber < LScene->RayCounter)
     {
// Mark this Triangle as already having been tested by this ray
//
      LRTriangle->LastRayNumber = LScene->RayCounter;

//      if(E3d_RayIntersectWithTriangleE(LRay, LRTriangle->V0, LRTriangle->Edge1, LRTriangle->Edge2, &LU, &LV))
      if(E3d_RayIntersectWithTriangleBackCullENoUV(LRay, LRTriangle->V0, LRTriangle->Edge1, LRTriangle->Edge2))
      {
       if(LRay->LastT<LDistance)
       {
	LScene->ShadowRayLastHit=LRTriangle;
	return(TRUE);
       }
      }
     }
    }
   }
  } while(E3d_RayDDANextVoxel(LRay));
 }

 return(FALSE);
}


//================================================================================
// Test a Ray against objects in the uniform-grid.
// This version just tests whether there is a hit.
//================================================================================
E3dRTriangle* E3d_RayIntersectVoxelsTest(E3dScene* LScene, E3dRay* LRay, void* LExcludeFace)
{
 E3dRTriangle*	LRTriangle;


 LScene->RayCounter++;




 if(LScene->TestRayLastHit)
 {
  LRTriangle=LScene->TestRayLastHit;

  if(LRTriangle->Face!=LExcludeFace)
  {
// Mark this Triangle as already having been tested by this ray
//
   LRTriangle->LastRayNumber = LScene->RayCounter;
   if(E3d_RayIntersectWithTriangleENoUV(LRay, LRTriangle->V0, LRTriangle->Edge1, LRTriangle->Edge2))
   {
//printf("TestLastHit %p\n", LRTriangle);fflush(stdout);
    return(LRTriangle);
   }
  }
 }


 E3d_RayDDAInit(LScene, LRay);

 {
  E3dVoxel*		LVoxel;
  int			LVoxelIndex, LPageSize = LScene->GridXCount*LScene->GridYCount;
  int			LRowSize = LScene->GridXCount;
  unsigned int		LC;


  do
  {
// Check the Triangles in this Voxel
//
   LVoxelIndex=LRay->GridZ*LPageSize + LRay->GridY*LRowSize + LRay->GridX;
   LVoxel = LScene->Voxels + LVoxelIndex;


   for(LC=0;LC<LVoxel->NumOfTriangles;LC++)
   {
    LRTriangle=LVoxel->Triangles[LC];

    if(LRTriangle->Face!=LExcludeFace)
    {
// Check if this triangle has been tested by this ray
//
     if(LRTriangle->LastRayNumber < LScene->RayCounter)
     {
// Mark this Triangle as already having been tested by this ray
//
      LRTriangle->LastRayNumber = LScene->RayCounter;

//      if(E3d_RayIntersectWithTriangleBackCullENoUV(LRay, LRTriangle->V0, LRTriangle->Edge1, LRTriangle->Edge2))
      if(E3d_RayIntersectWithTriangleENoUV(LRay, LRTriangle->V0, LRTriangle->Edge1, LRTriangle->Edge2))
      {
       LScene->TestRayLastHit=LRTriangle;
       return(LRTriangle);
      }
     }
    }
   }
  } while(E3d_RayDDANextVoxel(LRay));
 }

 return(NULL);
}


//================================================================================
// Test a Ray against objects in the uniform-grid.
// This version just tests whether there is a hit.
//================================================================================
E3dRTriangle* E3d_RayIntersectVoxelsTestStartSet(E3dScene* LScene, E3dRay* LRay, void* LExcludeFace)
{
 E3dRTriangle*	LRTriangle;


 LScene->RayCounter++;

 if(LScene->TestRayLastHit)
 {
  LRTriangle=LScene->TestRayLastHit;

  if(LRTriangle->Face!=LExcludeFace)
  {
// Mark this Triangle as already having been tested by this ray
//
   LRTriangle->LastRayNumber = LScene->RayCounter;
   if(E3d_RayIntersectWithTriangleENoUV(LRay, LRTriangle->V0, LRTriangle->Edge1, LRTriangle->Edge2))
   {
//printf("TestLastHit %p\n", LRTriangle);fflush(stdout);
    return(LRTriangle);
   }
  }
 }


 E3d_RayDDAInitStartSet(LScene, LRay);

 {
  E3dVoxel*		LVoxel;
  int			LVoxelIndex, LPageSize = LScene->GridXCount*LScene->GridYCount;
  int			LRowSize = LScene->GridXCount;
  unsigned int		LC;

  do
  {
// Check the Triangles in this Voxel
//
   LVoxelIndex=LRay->GridZ*LPageSize + LRay->GridY*LRowSize + LRay->GridX;
   LVoxel = LScene->Voxels + LVoxelIndex;

   for(LC=0;LC<LVoxel->NumOfTriangles;LC++)
   {
    LRTriangle=LVoxel->Triangles[LC];

    if(LRTriangle->Face!=LExcludeFace)
    {
// Check if this triangle has been tested by this ray
//
     if(LRTriangle->LastRayNumber < LScene->RayCounter)
     {
// Mark this Triangle as already having been tested by this ray
//
      LRTriangle->LastRayNumber = LScene->RayCounter;

//      if(E3d_RayIntersectWithTriangleBackCullENoUV(LRay, LRTriangle->V0, LRTriangle->Edge1, LRTriangle->Edge2))
      if(E3d_RayIntersectWithTriangleENoUV(LRay, LRTriangle->V0, LRTriangle->Edge1, LRTriangle->Edge2))
      {
       LScene->TestRayLastHit=LRTriangle;
       return(LRTriangle);
      }
     }
    }
   }
  } while(E3d_RayDDANextVoxel(LRay));
 }

 return(NULL);
}


//================================================================
// Compute illumination values (R,G,B) for a surface point
//================================================================
void E3d_IlluminatePoint(E3dScene* LScene, E3dSurfacePoint* LPoint, E3dMaterial* LMaterial, EcRGBAfColor* LIllumination)
{
 E3dFColorComponent	LR, LG, LB;
 int			LC;


 switch(LMaterial->Type)
 {
  default:	LR=LG=LB=0.0;break;

  case E3dMAT_CUSTOM:
   {
    E3dShader*	LShader=LMaterial->Shader;

    if(LShader)
    {
     E3dShaderSurfaceProc	LSurfaceProc=LShader->Class->SurfaceProc;

     if(LSurfaceProc)
     {
      LSurfaceProc(LScene, LPoint, LMaterial, LIllumination);
      return;
     }
    }
    LIllumination->R=0.0f;LIllumination->G=0.0f;LIllumination->B=0.0f;
    return;
   }

  case E3dMAT_CONSTANT:
   LR=LMaterial->Diffuse.R;
   LG=LMaterial->Diffuse.G;
   LB=LMaterial->Diffuse.B;
  break;

  case E3dMAT_LAMBERT:
   {
    E3dLight**		LLights=LScene->Lights;
    E3dFColorComponent	LDiffuseR, LDiffuseG, LDiffuseB;
    E3dCoordinate	LDX, LDY, LDZ, LNX, LNY, LNZ, LDistance, LDistInv;
    E3dCoordinate	LX, LY, LZ;
    E3dAngle		LDotProd;
    E3dLight*		LLight;
    double		LAttn;
    unsigned int	LNumOfLights=LScene->NumOfLights;


    LNX=LPoint->NormalX;LNY=LPoint->NormalY;LNZ=LPoint->NormalZ;
    LDiffuseR=LMaterial->Diffuse.R;
    LDiffuseG=LMaterial->Diffuse.G;
    LDiffuseB=LMaterial->Diffuse.B;

    LR=LMaterial->Ambient.R;
    LG=LMaterial->Ambient.G;
    LB=LMaterial->Ambient.B;

    LX=LPoint->X;LY=LPoint->Y;LZ=LPoint->Z;

    for(LC=0;LC<LNumOfLights;LC++)
    {
     LLight=LLights[LC];
     LDX=LLight->Position.X-LX;
     LDY=LLight->Position.Y-LY;
     LDZ=LLight->Position.Z-LZ;

     LDistance=LDX*LDX+LDY*LDY+LDZ*LDZ;
     if(LDistance>0.0)
     {
      LDistance=sqrt(LDistance);
      LDistInv=1.0/LDistance;
      LDX*=LDistInv;LDY*=LDistInv;LDZ*=LDistInv;	// LDX, LDY, LDZ: Normalized direction vector to the current light source

// If LD . LN <= 0.0 there is no diffuse or specular light at LPosition (no 'negative lighting')
//
      if((LDotProd=LDX*LNX+LDY*LNY+LDZ*LNZ)>0.0)
      {
       if(LLight->QuadraticAttenuation>0.0)
       {
	LAttn=1.0/(LLight->ConstAttenuation+(LLight->QuadraticAttenuation)*LDistance*LDistance);
	LDotProd*=LAttn;
       }
       LR+=LLight->Color.R*LDotProd*LDiffuseR;
       LG+=LLight->Color.G*LDotProd*LDiffuseG;
       LB+=LLight->Color.B*LDotProd*LDiffuseB;
      }
     }
    }
   }
  break;

  case E3dMAT_PHONG:
  case E3dMAT_BLINN:
   {
    E3dCamera*		LCamera=LScene->RenderInfo.Camera;
    E3dLight**		LLights=LScene->Lights;
    E3dFColorComponent	LDiffuseR, LDiffuseG, LDiffuseB, LSpecularR, LSpecularG, LSpecularB;
    E3dCoordinate	LDX, LDY, LDZ, LNX, LNY, LNZ, LVX, LVY, LVZ, LDistance, LDistInv, LL;
    E3dCoordinate	LX, LY, LZ;
    E3dAngle		LDotProd;
    E3dLight*		LLight;
    double		LAttn;
    unsigned int	LNumOfLights=LScene->NumOfLights;


    LNX=LPoint->NormalX;LNY=LPoint->NormalY;LNZ=LPoint->NormalZ;
    LDiffuseR=LMaterial->Diffuse.R;
    LDiffuseG=LMaterial->Diffuse.G;
    LDiffuseB=LMaterial->Diffuse.B;
    LSpecularR=LMaterial->Specular.R;
    LSpecularG=LMaterial->Specular.G;
    LSpecularB=LMaterial->Specular.B;

    LR=LMaterial->Ambient.R;
    LG=LMaterial->Ambient.G;
    LB=LMaterial->Ambient.B;

    LX=LPoint->X;LY=LPoint->Y;LZ=LPoint->Z;

// Calculate point-to-viewer unit vector
//
    LVX=LCamera->Position.X-LX;
    LVY=LCamera->Position.Y-LY;
    LVZ=LCamera->Position.Z-LZ;

// Normalize LV
//
    LL=sqrt(LVX*LVX+LVY*LVY+LVZ*LVZ);
    if(LL>0.0)
    {
     LL=1.0/LL;
     LVX*=LL;LVY*=LL;LVZ*=LL;
    }
    else	// LPosition is in the viewpoint
    {
     LIllumination->R=0.0;LIllumination->G=0.0;LIllumination->B=0.0;
     return;
    }

    for(LC=0;LC<LNumOfLights;LC++)
    {
     LLight=LLights[LC];
     LDX=LLight->Position.X-LX;
     LDY=LLight->Position.Y-LY;
     LDZ=LLight->Position.Z-LZ;

     LDistance=LDX*LDX+LDY*LDY+LDZ*LDZ;
     if(LDistance>0.0)
     {
      LDistance=sqrt(LDistance);
      LDistInv=1.0/LDistance;
      LDX*=LDistInv;LDY*=LDistInv;LDZ*=LDistInv;	// LDX, LDY, LDZ: Normalized direction vector to the current light source

// If LD . LN <= 0.0 there is no diffuse or specular light at LPosition (no 'negative lighting')
//
      if((LDotProd=LDX*LNX+LDY*LNY+LDZ*LNZ)>0.0)
      {
       E3dCoordinate	LSX, LSY, LSZ;

       if(LLight->QuadraticAttenuation>0.0)
       {
	LAttn=1.0/(LLight->ConstAttenuation+(LLight->QuadraticAttenuation)*LDistance*LDistance);
       }
       else LAttn=1.0;

// Diffuse component
//
       LDotProd*=LAttn;
       LR+=LLight->Color.R*LDotProd*LDiffuseR;
       LG+=LLight->Color.G*LDotProd*LDiffuseG;
       LB+=LLight->Color.B*LDotProd*LDiffuseB;

// Specular component
//

// Get sum of LPosition-to-viewpoint vector and Vertex to light source vector
//

// Add the LPosition-to-light source vector to LV
//
       LSX=LVX+LDX;
       LSY=LVY+LDY;
       LSZ=LVZ+LDZ;
     
// Normalize result
//
       LL=sqrt(LSX*LSX+LSY*LSY+LSZ*LSZ);
       if(LL>0.0)
       {
	LL=1.0/LL;
	LSX*=LL;LSY*=LL;LSZ*=LL;

	if((LDotProd=LSX*LNX+LSY*LNY+LSZ*LNZ)>0.0)
	{
	 LDotProd=pow(LDotProd, LMaterial->Specularity);
	 LR+=LLight->Color.R*LDotProd*LSpecularR;
	 LG+=LLight->Color.G*LDotProd*LSpecularG;
	 LB+=LLight->Color.B*LDotProd*LSpecularB;
	}
       }
      }
      else // The LPosition is in the viewpoint
      {
      }
     }
    }
   }
  break;
 }

 LIllumination->R=LR;LIllumination->G=LG;LIllumination->B=LB;
}



//========================================================
// Compute illumination values (R,G,B) for a SurfacePoint
// Simple shadow-ray testing
//========================================================
void E3d_IlluminatePointShadowing(E3dScene* LScene, E3dSurfacePoint* LPoint, E3dMaterial* LMaterial, EcRGBAfColor* LIllumination, E3dRTriangle* LExcludeRTriangle)
{
 E3dFColorComponent	LR, LG, LB;
 int			LC;


 switch(LMaterial->Type)
 {
  default:	LR=LG=LB=0.0;break;

  case E3dMAT_CUSTOM:
   {
    E3dShader*	LShader=LMaterial->Shader;

    if(LShader)
    {
     E3dShaderSurfaceProc	LSurfaceProc=LShader->Class->SurfaceProc;

     if(LSurfaceProc)
     {
      LSurfaceProc(LScene, LPoint, LMaterial, LIllumination);
      return;
     }
    }
    LIllumination->R=0.0f;LIllumination->G=0.0f;LIllumination->B=0.0f;
    return;
   }

  case E3dMAT_CONSTANT:
   LR=LMaterial->Diffuse.R;
   LG=LMaterial->Diffuse.G;
   LB=LMaterial->Diffuse.B;
  break;

  case E3dMAT_LAMBERT:
   {
    E3dLight**		LLights=LScene->Lights;
    E3dFColorComponent	LDiffuseR, LDiffuseG, LDiffuseB;
    E3dCoordinate	LDX, LDY, LDZ, LNX, LNY, LNZ, LDistance, LDistInv;
    E3dCoordinate	LX, LY, LZ;
    E3dAngle		LDotProd;
    E3dLight*		LLight;
    double		LAttn;

    E3dRay		LShadowRay;		// For shadows
    E3dFColorComponent	LDiffuseFactor;
    unsigned int	LNumOfLights=LScene->NumOfLights;


    LNX=LPoint->NormalX;LNY=LPoint->NormalY;LNZ=LPoint->NormalZ;
    LDiffuseR=LMaterial->Diffuse.R;
    LDiffuseG=LMaterial->Diffuse.G;
    LDiffuseB=LMaterial->Diffuse.B;

    LR=LMaterial->Ambient.R;
    LG=LMaterial->Ambient.G;
    LB=LMaterial->Ambient.B;

    LX=LPoint->X;LY=LPoint->Y;LZ=LPoint->Z;

    for(LC=0;LC<LNumOfLights;LC++)
    {
     LLight=LLights[LC];
     LDX=LLight->Position.X-LX;
     LDY=LLight->Position.Y-LY;
     LDZ=LLight->Position.Z-LZ;

     LDistance=LDX*LDX+LDY*LDY+LDZ*LDZ;
     if(LDistance>0.0)
     {
      LDistance=sqrt(LDistance);
      LDistInv=1.0/LDistance;
      LDX*=LDistInv;LDY*=LDistInv;LDZ*=LDistInv;	// LDX, LDY, LDZ: Normalized direction vector to the current light source


// Fire a shadow-test ray
//
      if(LLight->ShadowBrightness<1.0)
      {
       LShadowRay.X=LX;
       LShadowRay.Y=LY;
       LShadowRay.Z=LZ;
       LShadowRay.DirX=LDX;
       LShadowRay.DirY=LDY;
       LShadowRay.DirZ=LDZ;
       if(E3d_ShadowRayIntersectVoxels(LScene, &LShadowRay, LExcludeRTriangle->Face, LDistance))
       {
	LDiffuseFactor=LLight->ShadowBrightness;
       }
       else
       {
	LDiffuseFactor=1.0;
       }
      }
      else LDiffuseFactor=1.0;

// If LD . LN <= 0.0 there is no diffuse or specular light at LPosition (no 'negative lighting')
//
      if((LDotProd=LDX*LNX+LDY*LNY+LDZ*LNZ)>0.0)
      {
       if(LLight->QuadraticAttenuation>0.0)
       {
	LAttn=1.0/(LLight->ConstAttenuation+(LLight->QuadraticAttenuation)*LDistance*LDistance);
	LDotProd*=LAttn;
       }
       LDotProd*=LDiffuseFactor;
       LR+=LLight->Color.R*LDotProd*LDiffuseR;
       LG+=LLight->Color.G*LDotProd*LDiffuseG;
       LB+=LLight->Color.B*LDotProd*LDiffuseB;
      }
     }
    }
   }
  break;

  case E3dMAT_PHONG:
  case E3dMAT_BLINN:
   {
    E3dCamera*		LCamera=LScene->RenderInfo.Camera;
    E3dLight**		LLights=LScene->Lights;
    E3dLight*		LLight;
    E3dFColorComponent	LDiffuseR, LDiffuseG, LDiffuseB, LSpecularR, LSpecularG, LSpecularB;
    E3dCoordinate	LDX, LDY, LDZ, LNX, LNY, LNZ, LVX, LVY, LVZ, LDistance, LDistInv, LL;
    E3dCoordinate	LX, LY, LZ;
    E3dAngle		LDotProd;
    double		LAttn;

    E3dRay		LShadowRay;
    E3dFColorComponent	LDiffuseFactor, LSpecularFactor;	// For shadows
    unsigned int	LNumOfLights=LScene->NumOfLights;


    LNX=LPoint->NormalX;LNY=LPoint->NormalY;LNZ=LPoint->NormalZ;
    LDiffuseR=LMaterial->Diffuse.R;
    LDiffuseG=LMaterial->Diffuse.G;
    LDiffuseB=LMaterial->Diffuse.B;
    LSpecularR=LMaterial->Specular.R;
    LSpecularG=LMaterial->Specular.G;
    LSpecularB=LMaterial->Specular.B;

    LR=LMaterial->Ambient.R;
    LG=LMaterial->Ambient.G;
    LB=LMaterial->Ambient.B;

    LX=LPoint->X;LY=LPoint->Y;LZ=LPoint->Z;

// Calculate point-to-viewer unit vector
//
    LVX=LCamera->Position.X-LX;
    LVY=LCamera->Position.Y-LY;
    LVZ=LCamera->Position.Z-LZ;

// Normalize LV
//
    LL=sqrt(LVX*LVX+LVY*LVY+LVZ*LVZ);
    if(LL>0.0)
    {
     LL=1.0/LL;
     LVX*=LL;LVY*=LL;LVZ*=LL;
    }
    else	// LPosition is in the viewpoint
    {
     LIllumination->R=0.0;LIllumination->G=0.0;LIllumination->B=0.0;
     return;
    }

    for(LC=0;LC<LNumOfLights;LC++)
    {
     LLight=LLights[LC];
     LDX=LLight->Position.X-LX;
     LDY=LLight->Position.Y-LY;
     LDZ=LLight->Position.Z-LZ;

//Printf("Light %f,%f,%f\n", LLight->Position.X, LLight->Position.Y, LLight->Position.Z);fflush(stdout);

     LDistance=LDX*LDX+LDY*LDY+LDZ*LDZ;
     if(LDistance!=0.0)
     {
      LDistance=sqrt(LDistance);
      LDistInv=1.0/LDistance;
      LDX*=LDistInv;LDY*=LDistInv;LDZ*=LDistInv;	// LDX, LDY, LDZ: Normalized direction vector to the current light source


// Fire a shadow-test ray
//
      if(LLight->ShadowBrightness<1.0)
      {
       LShadowRay.X=LX;
       LShadowRay.Y=LY;
       LShadowRay.Z=LZ;
       LShadowRay.DirX=LDX;
       LShadowRay.DirY=LDY;
       LShadowRay.DirZ=LDZ;
       if(E3d_ShadowRayIntersectVoxels(LScene, &LShadowRay, LExcludeRTriangle->Face, LDistance))
       {
	LDiffuseFactor=LLight->ShadowBrightness;
	LSpecularFactor=0.0;
       }
       else
       {
	LDiffuseFactor=1.0;
	LSpecularFactor=1.0;
       }
      }
      else
      {
       LDiffuseFactor=1.0;
       LSpecularFactor=1.0;
      }


// If LD . LN <= 0.0 there is no diffuse or specular light at LPosition (no 'negative lighting')
//
      if((LDotProd=LDX*LNX+LDY*LNY+LDZ*LNZ)>0.0)
      {
       E3dCoordinate	LSX, LSY, LSZ;

       if(LLight->QuadraticAttenuation>0.0)
       {
	LAttn=1.0/(LLight->ConstAttenuation+(LLight->QuadraticAttenuation)*LDistance*LDistance);
       }
       else LAttn=1.0;

// Diffuse component
//
       LDotProd*=LAttn*LDiffuseFactor;
       LR+=LLight->Color.R*LDotProd*LDiffuseR;
       LG+=LLight->Color.G*LDotProd*LDiffuseG;
       LB+=LLight->Color.B*LDotProd*LDiffuseB;

// Specular component
//

// Get sum of LPosition-to-viewpoint vector and Vertex to light source vector
//

// Add the LPosition-to-light vector to LV
//
       LSX=LVX+LDX;
       LSY=LVY+LDY;
       LSZ=LVZ+LDZ;
     
// Normalize result
//
       LL=sqrt(LSX*LSX+LSY*LSY+LSZ*LSZ);
       if(LL>0.0)
       {
	LL=1.0/LL;
	LSX*=LL;LSY*=LL;LSZ*=LL;

	if((LDotProd=LSX*LNX+LSY*LNY+LSZ*LNZ)>0.0)
	{
	 LDotProd=pow(LDotProd, LMaterial->Specularity)*LSpecularFactor;
	 LR+=LLight->Color.R*LDotProd*LSpecularR;
	 LG+=LLight->Color.G*LDotProd*LSpecularG;
	 LB+=LLight->Color.B*LDotProd*LSpecularB;
	}
       }
      }
      else // LPosition is in the viewpoint
      {
      }
     }
    }
   }
  break;
 }

 LIllumination->R=LR;LIllumination->G=LG;LIllumination->B=LB;
}
