---

RT_Grid.cc

Go to the documentation of this file.

/*
    File:       RT_Grid.cc

    Function:   Defines a hierarchical grid class for ray-tracing 
                optimisation

    Author:     Andrew Willmott, ajw+@cs.cmu.edu
*/

#include "RT_Grid.h"
#include "RT_RayStats.h"
#include "FCompare.h"
#ifdef RT_VIS
#include "gcl/Renderer.h"
#endif

#define RT_SetLock(x) 
#define RT_UnsetLock(x)
#define RT_Lock Int

#define DBG_COUT if (0) cerr


// --- Globals ----------------------------------------------------------------


Int     RT_Grid::sStamp = 1;
Int     RT_Grid::sTag = 0x00000001;     
                                // By default, the first group is enabled.
Int     RT_Grid::sMem = 0;
Int     RT_Grid::sTriMem = 0;
Int     RT_Grid::sGenMem = 0;
Int     RT_Grid::sPointMem = 0;
Int     RT_Grid::sMemLimit = 0;

Int     RT_Grid::sMaxCells      = 40;   
                                // Maximum cells along one side of a side
Int     RT_Grid::sMaxCellPrims  = RT_Grid::sMaxCells;
                                // Need this many prims in a cell before
                                // considering them as subgrid candidates.

GCLReal RT_Grid::sPushSizeLimit = 0.05;
                                // A primitive must be smaller than this 
                                // fraction of a cell for it to be a subgrid
                                // candidate.
GCLReal RT_Grid::sMaxCellRatio  = 0.20;
                                // A subgrid's cell size must be this much
                                // smaller than its parent's.
GCLReal RT_Grid::sMinDensity    = 0.05;
                                // A grid must have at least this primitive
                                // density to be considered worthwhile creating
                                // it.
GCLReal RT_Grid::sMaxExpand     = 3.00;
                                // The maximum amount to expand a subgrid in
                                // order to merge it with another.
GCLReal RT_Grid::sEpsilon       = 1e-15;
                                // Nearness epsilon. For intersection 
                                // tests where a casting primitive isn't
                                // provided, this helps avoid self-shadowing /
                                // intersection.


// --- Internal data structures -----------------------------------------------


struct PrimEntry
{
    union
    {
    RT_Prim*    prim;
    RT_Grid*    grid;
    };
    PrimEntry*  next;
};

typedef PrimEntry* PrimEntryPtr;

// We assume no memory is malloc'd at $8.

const PrimEntryPtr      PRIM_GRID_TAG = (PrimEntry*) 0x8;   

struct PrimListEntry
{
    PrimEntry**     list;
    PrimListEntry*  next;
};


// --- GridStore methods ------------------------------------------------------


Void GridStore::Init(Int size)
{
    Int     spill;
    Int     i;
    
    spill = size & PAGE_PTR_MASK;   // size of last page
    if (spill == 0)
        spill = PAGE_PTR_SIZE;
        
    numPages = ((size - 1) >> PAGE_PTR_BITS) + 1;
    pages = RT_ArrayAlloc(PrimEntry**, numPages);
    mem = sizeof(PrimEntry**) * numPages;
    
    for (i = 0; i < numPages - 1; i++)
    {
        pages[i] = RT_ArrayAllocClr(PrimEntry*, PAGE_PTR_SIZE);
        Assert(pages[i] != 0, "Out of memory");
        mem += sizeof(PrimEntry*) * PAGE_PTR_SIZE;
    }
    
    pages[numPages - 1] = RT_ArrayAllocClr(PrimEntry*, spill);
    mem += sizeof(PrimEntry*) * spill;
}

Void GridStore::Free()
{
    Int     i;
    
    if (pages)
    {
        for (i = 0; i < numPages; i++)
            RT_ArrayFree(pages[i]);
        
        RT_ArrayFree(pages);
    }
}

// --- Useful #defines --------------------------------------------------------


#ifdef RT_USE_FC
#define FEQG(a,b) F_EQG10(a,b)
#define FSGN(x) ((F_EQ10((x), 0)) ? 0 : (((x) < 0) ? -1 : 1))
#else
#define FEQG(a, b) (b >= a)
// the epsilon is necessary here so that planes that just touch grid cells
// get added to them... old old bug.
#define FSGN(x) (((x) > sEpsilon ) ? 1 : (((x) < -sEpsilon) ? -1 : 0))
#endif

// --- RT_Grid methods --------------------------------------------------------


Void RT_Grid::MasterInit()
{
    sMem = 0;
    sTriMem = 0;
    sGenMem = 0;
    sPointMem = 0;
    Init();
}

Void RT_Grid::Init()
{
    totalPrims = 0;
    dirty = false;
    level = 0;

    max.MakeBlock(-vl_inf);
    min.MakeBlock(vl_inf);

    totalCells = 0;
    addPrims = 0;
    addObjs = 0;
    addedObjs = 0;
    subGrids = 0;
    next = 0;
    hitPrim = 0;
    hitT = 0.0;
    stamp = 0;
    parent = 0;

    sMem += sizeof(RT_Grid);
}

Void RT_Grid::Free()
// Delete subfields of the grid.
{
    RT_Grid         *subGrid, *nextGrid = 0;
    Int             i;
    PrimEntry       *entry, *nextEntry = 0;
    PrimListEntry   *tlEntry, *tlNextEntry = 0;
    ObjEntry        *objEntry, *objNextEntry = 0;
    
    // Delete subgrids
        
    for (subGrid = subGrids; subGrid; subGrid = nextGrid)
    {
        nextGrid = subGrid->next;
        subGrid->Free();
        RT_Free(subGrid);
    }

    // Delete the grid storage

    if (totalCells > 0)
    {
        for (i = 0; i < totalCells; i++)
            for (entry = grid[i]; entry && entry != PRIM_GRID_TAG; 
                 entry = nextEntry)
            {
                nextEntry = entry->next;
                RT_Free(entry);
            }
             
        grid.Free();
    }
    
    for (tlEntry = addPrims; tlEntry; tlEntry = tlNextEntry)
    {
        tlNextEntry = tlEntry->next;
        RT_Free(tlEntry);
    }
    
    for (objEntry = addObjs; objEntry; objEntry = objNextEntry)
    {
        objNextEntry = objEntry->next;
        RT_Free(objEntry);
    }   

    for (objEntry = addedObjs; objEntry; objEntry = objNextEntry)
    {
        objNextEntry = objEntry->next;
        RT_Free(objEntry);
    }   
}

Void RT_Grid::AddObject(RT_Object *object)
{
    ObjEntry*   newEntry;
        
    UpdateBounds(object);

    newEntry = RT_Alloc(ObjEntry);
    sMem += sizeof(ObjEntry);
    
    newEntry->object = object;
    newEntry->next = addObjs;
    addObjs = newEntry;
    
    dirty = true;
}

RT_Lock *gGridLock = 0;

Void RT_Grid::PrepareGrid()
// Prepare grid for raytracing.
{
    ObjEntry    *entry, *nextEntry;
    Int         j;
    
    if (gGridLock)
        RT_SetLock(*gGridLock);

    if (!dirty)
    {
        // Another process has Prepare()'d this grid while we were
        // waiting for the lock...
        
        if (gGridLock)
            RT_UnsetLock(*gGridLock);
        return; 
    }
    
    if (addObjs)
    {
        RT_Tri*     tri;
        RT_GenPtr*  gen;

        SetupGrid();

        // Add object triangles
        
        for (entry = addObjs; entry; entry = nextEntry)
        {
            nextEntry = entry->next;    

            tri = entry->object->tris;
            for (j = 0; j < entry->object->numTris; j++, tri++)
                AddTriangle(tri);
            totalPrims += entry->object->numTris;
        
            gen = entry->object->gens;
            for (j = 0; j < entry->object->numGens; j++, gen++)
                AddGeneric(*gen);
            totalPrims += entry->object->numGens;
        
            sTriMem += sizeof(RT_Tri) * entry->object->numTris;
            sGenMem += sizeof(RT_GenPtr) * entry->object->numGens;
            sPointMem += entry->object->numPoints * sizeof(Point);
            
            if (nextEntry == 0)
                entry->next = addedObjs;
        }

        // transfer objects to the added objects list

        addedObjs = addObjs;
        addObjs = 0;
    }
                
    // Now, create nested grids.
    
    dirty = false;

    if (sMemLimit && sMem >= sMemLimit)
    {
        cerr << "*** out of memory: aborting grid creation... >" 
             << level << endl;
        cerr << "*** memory = " << sMem << " bytes." << endl;
        if (gGridLock)
            RT_UnsetLock(*gGridLock);
        return;
    }
    
    CreateNestedGrids();

    if (gGridLock)
        RT_UnsetLock(*gGridLock);
}


Void RT_Grid::UpdateBounds(Point& cellMin, Point& cellMax)
{
    FindMinElts(cellMin, min, min);
    FindMaxElts(cellMax, max, max);
}

Void RT_Grid::UpdateBounds(RT_Object *object)
{
    Int     i;
    Point*  points = object->points;
            
    for (i = 0; i < object->numPoints; i++)
		::UpdateBounds(points[i], min, max);
    
    for (i = 0; i < object->numGens; i++)
        object->gens[i]->UpdateBounds(min, max);
}

Void RT_Grid::SetupGrid()
{
    GCLReal maxSpan;
    Vector  centre, span;
    Int     i;

    // Calculate the boundaries for the grid.  We already know
    // the world-space boundaries for the object.

    // XXX fix more elegantly
    // expand the box a little so we don't suffer from numerical
    // problems for surfaces lying along one of the sides of the grid
    const GCLReal   kGridFudge = 1e-15;
    Vector          fudge;

    max += fudge.MakeBlock(kGridFudge);
    min -= fudge;

    span = max - min;
    maxSpan = MaxElt(span);
    cellSize = maxSpan / sMaxCells;
    areaLimit = sPushSizeLimit * sqr(cellSize);     
    
    // Calculate how many grid cells the object will have in each 
    // direction.  We've already assigned the widest direction a 
    // preset number of grid cells and the narrower directions fewer 
    // cells.
    
    for (i = 0; i < 3; i++)
    {
        numCells[i] = (Int) ceil(sMaxCells * span[i] / maxSpan);

        // Ensure at least 1 cell in all directions
        if (numCells[i] == 0)       
            numCells[i] = 1;
        
        max[i] = min[i] + numCells[i] * cellSize;
    }

    // Allocate a one-dimensional array for the grid.

    totalCells = numCells[0] * numCells[1] * numCells[2];
    xSpan = numCells[1] * numCells[2];
    ySpan = numCells[2];

    grid.Init(totalCells);
    sMem += (Int) grid.MemoryUsage();
}


Void RT_Grid::AddGeneric(RT_GenPtr gen)
// Add a generic object to this grid.
{
    Vector      tMin, tMax;
    Vector      cellMin, vertex;
    Int         i;
    Int         x, y, z, cellNum, cellsAdded;
    PrimEntry   **last, *entry, *newEntry;
    Int         start[3], end[3];
    
    // We want to enter the polygon into all of the grid cells 
    // it intersects.  We first limit our test of grid cells to 
    // those in which the polygon's bounding box resides.

    gen->FindBounds(tMin, tMax);
            
    // Find the equivalent cell bounds

    for (i = 0; i < 3; i++)
    {
        start[i] = (Int) ((tMin[i] - min[i]) / cellSize);
        end[i]   = (Int) ((tMax[i] - min[i]) / cellSize);

        // Crop to the grid

        if (end[i] >= numCells[i])
            end[i] = numCells[i] - 1;
        else if (end[i] < 0)
            end[i] = 0;
            
        if (start[i] >= numCells[i])
            start[i] = numCells[i] - 1;
        else if (start[i] < 0)
            start[i] = 0;
    }
        
    // Finally, the generic is added to every grid cell in its bounding box

    cellsAdded = 0;
    for (x = start[0]; x <= end[0]; x++)
        for (y = start[1]; y <= end[1]; y++)
            for (z = start[2]; z <= end[2]; z++)
            {
                cellsAdded++;
                newEntry = RT_Alloc(PrimEntry);
                sMem += sizeof(PrimEntry);
                cellNum = x * xSpan + y * ySpan + z;
                
                newEntry->prim = gen;

                // Insert gen in order of (equivalent) area...
                
                last = &(grid[cellNum]);    // the last entry we looked at
                
                for (entry = *last; entry; entry = entry->next)
                {
                    // Only bother sorting prims < areaLimit
                    
                    if (entry->prim->area <= gen->area || 
                        entry->prim->area < areaLimit) 
                        break;
                        
                    last = &(entry->next);
                }
                
                *last = newEntry;
                newEntry->next = entry;
            }

#ifdef RT_OPAC
    gen->cells = cellsAdded;
#endif
}


Void RT_Grid::AddTriangle(RT_Tri *tri)
// Add a triangle to this grid.
{
    Vector      tMin, tMax;
    Vector      cellMin, vertex;
    Point       *p1, *p2, *p3;
    Int         i;
    Int         x, y, z, sign, cellNum, cellsAdded;
    PrimEntry   **last, *entry, *newEntry;
    Point*      points = tri->object->points;
    Int         start[3], end[3];
    GCLReal     temp;
    
    // We want to enter the polygon into all of the grid cells 
    // it intersects.  We first limit our test of grid cells to 
    // those in which the polygon's bounding box resides.

    dirty = true;
    p1 = points + tri->v[0];
    p2 = points + tri->v[1];
    p3 = points + tri->v[2];

    FindMaxElts(*p1, *p2, tMax);
    FindMaxElts(*p3, tMax, tMax);
    
    FindMinElts(*p1, *p2, tMin);
    FindMinElts(*p3, tMin, tMin);
        
    // Find the equivalent cell bounds

    for (i = 0; i < 3; i++)
    {
        start[i] = (Int) ((tMin[i] - min[i]) / cellSize);
        end[i]   = (Int) ((tMax[i] - min[i]) / cellSize);

        // Crop to the grid

        if (end[i] >= numCells[i])
            end[i] = numCells[i] - 1;
        else if (end[i] < 0)
            end[i] = 0;
            
        if (start[i] >= numCells[i])
            start[i] = numCells[i] - 1;
        else if (start[i] < 0)
            start[i] = 0;
    }
    
    // Finally, the triangle is added to every grid cell in the bounding box
    // which straddles the plane of the triangle. 

    cellsAdded = 0;
    for (x = start[0]; x <= end[0]; x++)
        for (y = start[1]; y <= end[1]; y++)
            for (z = start[2]; z <= end[2]; z++)
            {
                // We test each vertex of the cell by plugging it into the
                // triangle's plane equation. If a vertex has a different
                // sign from the preceeding vertices, we know the plane
                // straddles this cube, and we add the triangle to the
                // cell. A result of zero means the plane touches the vertex, 
                // which is also grounds for adding the triangle. 
                
                cellMin = min + Vector(x, y, z) * cellSize;
                cellNum = x * xSpan + y * ySpan + z;
                temp = dot(tri->normal, cellMin) + tri->d;
                sign = FSGN(temp);
                
                if (sign == 0)  // The first sign is zero: we can stop.
                    i = 1;
                else
                    for (i = 1; i < 8; i++) 
                    {
                        vertex = cellMin;
                        
                        if (i & 1) vertex[0] += cellSize;
                        if (i & 2) vertex[1] += cellSize;
                        if (i & 4) vertex[2] += cellSize;
                        
                        temp = dot(tri->normal, vertex) + tri->d;

                        if (sign != FSGN(temp))
                            break;
                    }

                // Add the triangle if we have a straddle

                if (i < 8)
                {
                    cellsAdded++;
                    newEntry = RT_Alloc(PrimEntry);
                    sMem += sizeof(PrimEntry);
                    
                    newEntry->prim = tri;

                    // Insert tri in order of area...
                    
                    last = &(grid[cellNum]);    // the last entry we looked at
                    
                    for (entry = *last; entry; entry = entry->next)
                    {
                        // Only bother sorting tris < areaLimit
                        
                        if (entry->prim->area <= tri->area || 
                            entry->prim->area < areaLimit) 
                            break;
                            
                        last = &(entry->next);
                    }
                    
                    *last = newEntry;
                    newEntry->next = entry;
                }
            }

#ifdef RT_OPAC
    tri->cells = cellsAdded;
#endif
}

Void RT_Grid::RayClampEntryPoint(Int index, Bool backEntry[3],
                                        Vector& gridPoint, Int cell[3])
{
    if (backEntry[index])
    {
        // If we're entering from the back end, we're entering the 
        // largest numbered cell in this direction. 

        cell[index] = numCells[index] - 1;
        gridPoint[index] = (GCLReal) numCells[index];
    }
    else
    {
        // From the front, it's the smallest, which is conveniently zero.

        cell[index] = 0;
        gridPoint[index] = 0.0;
    }
}  


Bool RT_Grid::FindGridStart(

    const Point&    start,
    const Vector&   direction,
    Int             cell[3],
    Int             increment[3],
    Int             limit[3], 
    GCLReal&        tRay,       
                    // start of ray in grid 
    Vector&         tDelta,     
                    // delta t for a step of size cellSize in each direction
    Vector&         tMax 
                    // delta t to the next grid wall in each direction
    )
{
    Int     i;
    Vector  gridPoint;
    Vector  in, out;
    Bool    backEntry[3];
    Point   point;
                
    gridPoint = (start - min) / cellSize;

    if (0.0 <= gridPoint[0] && gridPoint[0] <= numCells[0] &&
        0.0 <= gridPoint[1] && gridPoint[1] <= numCells[1] &&
        0.0 <= gridPoint[2] && gridPoint[2] <= numCells[2])
    {
        // We are starting inside the grid: find which cell
        
        tRay = 0.0;
        
        for (i = 0; i < 3; i++)
        {
            cell[i] = (Int) floor(gridPoint[i]);

            if (cell[i] == numCells[i])
                cell[i]--;
        }
    }
    else
    {
        // Otherwise, we have to compute the entry point
        // into the grid.

        for (i = 0; i < 3; i++)
        {
            // Avoid the divide by zero.  This would occur if the ray
            // has no component in one or two directions.
    
            if (abs(direction[i]) > (GCLReal) 0.0)
            {
                in[i]     = (min[i] - start[i]) / direction[i];
                out[i]    = (max[i] - start[i]) / direction[i];
            }
            else
            {
                in[i]     = -vl_inf;
                out[i]    =  vl_inf;
            }

            // We're really concerned only with the entry point of the ray 
            // into the grid.  Thus, if we've got our entry and exit points
            // backwards, then correct the entry point.
    
            if (backEntry[i] = (in[i] > out[i]))
                in[i] = out[i];
        }

        // The actual entry point into the grid occurs at the maximum of the
        // just-computed entry values (i.e., we enter the grid when we've
        // crossed the X _and_ the Y _and_ the Z planes of the Fujimoto 
        // grid.

        tRay = MaxElt(in);          
        point = start + tRay * direction;
    
        // Compute the entry point in terms of grid space.  (I.e., in
        // what cell do we enter the grid?)
        
        gridPoint = (point - min) / cellSize;
        
        for (i = 0; i < 3; i++)
            cell[i] = (Int) floor(gridPoint[i]);
        
        // the cell direction corresponding to the maximum component(s)
        // of 'in' will lie along the boundary of the grid. To 
        // prevent spilling over the boundary due to round-off error, we clamp
        // the direction(s) back to the grid boundary.
        
        if (in[0] > in[1])  
        {
            if (in[0] >= in[2]) 
            {
                RayClampEntryPoint(0, backEntry, gridPoint, cell);
    
                if (in[0] == in[2]) 
                    RayClampEntryPoint(2, backEntry, gridPoint, cell);
            }
            else
                RayClampEntryPoint(2, backEntry, gridPoint, cell);
        }
        else if (in[0] == in[1])
        {
            if (in[2] >= in[0])
            {
                RayClampEntryPoint(2, backEntry, gridPoint, cell);
    
                if (in[2] == in[0])
                {
                    RayClampEntryPoint(0, backEntry, gridPoint, cell);
                    RayClampEntryPoint(1, backEntry, gridPoint, cell);
                }
            }
            else
            {
                RayClampEntryPoint(0, backEntry, gridPoint, cell);
                RayClampEntryPoint(1, backEntry, gridPoint, cell);
            }
        }
        else
        {
            if (in[1] >= in[2])
            {
                RayClampEntryPoint(1, backEntry, gridPoint, cell);
    
                if (in[1] == in[2])
                    RayClampEntryPoint(2, backEntry, gridPoint, cell);
            }
            else
                RayClampEntryPoint(2, backEntry, gridPoint, cell);
        }
        
        // If the ray completely misses the Fujimoto grid, then we 
        // of course didn't hit anything.

        for (i = 0; i < 3; i++)
            if (cell[i] < 0 || numCells[i] <= cell[i])
                return(false);
    }
    
    // Calculate auxillary info used for traversing the grid
    
    for (i = 0; i < 3; i++)
    {
        GCLReal scale = cellSize;
    
        // The ray that goes through the Fujimoto grid is parameterized,
        // i.e., the ray is of the form start + t * direction, where t is
        // the parameter.  For any given direction (i.e. X, Y, or Z),
        // the amount that t changes as the ray crosses to the next cell
        // in the same direction is fixed.  Calculate what that change in
        // t is and store the information.  We're going to use this 
        // information to help us traverse the grid.

        tDelta[i] = (abs(direction[i]) == (GCLReal) 0.0) ?
            0.0 : scale / abs(direction[i]);
    
        // Note the sign of the direction vector of the ray.

        increment[i] = FSGN(direction[i]);
        limit[i] = (increment[i] < 0) ? -1 : numCells[i];

        // Make a note of what the value of the ray parameter t will
        // be when we cross over the next orthogonal plane in this 
        // direction.  (Of course, if the direction vector has no 
        // component in this direction, then we'll _never_ reach the 
        // next grid cell in this direction.)

        if (increment[i] == 0)
            tMax[i] = vl_inf;
        else if (increment[i] == 1)
            tMax[i] = tRay + tDelta[i] * (cell[i] + 1 - gridPoint[i]);
        else
        {
            // The the grid point is integral, i.e. if it lies exactly on
            // one of the orthogonal planes, then the next movement will
            // cross the entire cell.
        
            if (gridPoint[i] == floor(gridPoint[i]))
                tMax[i] = tRay + tDelta[i]; 
            else
                tMax[i] = tRay + tDelta[i] * 
                    (gridPoint[i] - floor(gridPoint[i]));
        }
    }
    
    return(tDelta[0] != 0.0 || tDelta[1] != 0.0 || tDelta[2] != 0);
}


// --- Intersection methods ---------------------------------------------------


Bool RT_Grid::IntersectionExists(
        const Point&    start,
        const Point&    end,
        RT_Prim*        origPrim,
        RayStats*       stats
    )
// Finds out if any intersection occurs between start and end. Ignores 
// intersections with origPrim, if it is not nil.
{
    Vector      direction = end - start;
    Vector      reflectPoint;
    GCLReal     t, tRay;
    Int         i, increment[3], limit[3];  
    Vector      tDelta, tMax;
    Int         nextCellDir;
    Int         cell[3];
    PrimEntry*  entry;
    
    if (stats) stats->rays++;

    if (sqrlen(direction) < 1e-15)
        return(false);

    if (dirty)
        PrepareGrid();

    if (level == 0)         
    {
        // just starting off: ignore all previously stamped subgrids.
        RT_Grid::IncTime();
        RT_Prim::IncTime();
    }
    else if (IsStamped())
    {
        // return cached result
        if (stats) stats->d2++;
        return(hitPrim != 0);
    }
    else
    {
#ifdef RT_GRID_CACHE
        // ensure result is cached
        Stamp();
#endif
    }
    
    if (stats) stats->grids++;
    hitPrim = 0;

    // Figure out grid parameters, and return if we completely miss the grid
    if (!FindGridStart(start, direction, cell, increment, limit, tRay,
                       tDelta, tMax))
        return(false);
            
    // Traverse the grid...

    for (i = 0; i < 256; i++)
    {
        // We are only after a hit, not the closest hit, so as soon as we get
        // one, we can return.
        if (stats) stats->voxels++;

        entry = grid[cell[0] * xSpan + cell[1] * ySpan + cell[2]];

        for (; entry; entry = entry->next)
        {
            if (stats) stats->polysTested++;

            if (entry->next == PRIM_GRID_TAG)
            {
                // If this is a subgrid entry, call the subgrid...

                if (entry->grid->IntersectionExists(start, end,
                                                    origPrim, stats))
                {
                    hitPrim = entry->grid->hitPrim;
                    hitT = entry->grid->hitT;
                    return(true);
                }
                
                break;  // a grid entry is always last.
            }
            
            if (entry->prim == origPrim)
                continue;
            if ((entry->prim->flags & TRI_INACTIVE) || !(sTag & entry->prim->object->tag))
                continue;
#ifdef RT_SUPPORT_GENS
            if (entry->prim->primType == rt_gen)
            {
                RT_GenPtr gen = *(RT_GenPtr*) entry->prim;
                                
                if (!gen->IsStamped())
                    gen->Intersect(start, direction, 0, 1);
                
                if (gen->IsTrue(GEN_HIT))
                {
                    hitPrim = gen;
                    hitT = gen->hitT;
                    DBG_COUT << "hit generic at " << hitT << endl;
                    return(true);
                }
            }
            else if (entry->prim->primType == rt_triangle)
#else
            if(1)
#endif
            {
                RT_Tri*     thisTri = (RT_Tri*) entry->prim;
                    
                // Find ray intersection with tri's plane. Bail if there 
                // isn't one.
                GCLReal normDotDirection = dot(thisTri->normal, direction);
    
                if (thisTri->flags & TRI_2SIDED_V)
                {
                    if (normDotDirection == 0.0)
                        continue;
                }
                else
                    if (normDotDirection >= 0.0)
                        continue;

                t = (thisTri->d + dot(thisTri->normal, start)) / 
                    -normDotDirection;

                // Bail if the intersection occurs before 'start' or after
                // 'end'. the epsilon here accounts for self-shadowing.
                if (t <= sEpsilon || t >= (1.0 - sEpsilon))
                    continue;

                // Okay, let's do the actual test. If we've hit the tri 
                // itself, we can return.

                if (stats) stats->intersections++;

                reflectPoint = start + t * direction;
                thisTri->PointIsInside(reflectPoint);
                        
                if (thisTri->IsTrue(TRI_HIT))
                {
                    if (stats) stats->goodHits++;
                    hitPrim = thisTri;
                    hitT = t;
                    return(true);
                }
            }
            else
                DBG_COUT << "unknown primitive type: " << 
                    entry->prim->primType << endl;
        }   

        nextCellDir = MinEltIndex(tMax);
        cell[nextCellDir] += increment[nextCellDir];

        // If we've exited the grid, then bail -- we didn't find any
        // intersections.
        if (cell[nextCellDir] == limit[nextCellDir])
            return(false);
        else
            tMax[nextCellDir] += tDelta[nextCellDir];
    }

#ifdef DEBUG
    cerr << "*** STUCK IN GRID: tMax: " << tMax << tDelta << " for " << start << end << endl;
#endif
}

RT_Prim *RT_Grid::ClosestIntersection(
        const Point&    start,              // Start of the ray
        const Vector&   direction,          // direction of the same
        RT_Prim*        origPrim,           // originating primitive
        RayStats*       stats               // ray statistics.
    )
// Finds the closest intersection to the start point. Ignores 
// intersections with origPrim, if it is not nil.
{   
    Vector      reflectPoint;
    Int         increment[3],  limit[3];    
    GCLReal     tRayNext, tRayLast, t, tMin;
    Vector      tDelta, tMax;
    Int         nextCellDir;
    Int         cell[3];
    PrimEntry*  entry;

    if (stats) stats->rays++;

    if (dirty)
        PrepareGrid();

    if (level == 0)         
    {
        // just starting off: ignore all previously stamped subgrids + prims
        DBG_COUT << "clearing stamps" << endl;
        RT_Prim::IncTime();
        RT_Grid::IncTime();
    }
    else if (IsStamped())
    {
        // return cached result
        DBG_COUT << "returning cached result" << endl;
        if (stats) stats->d2++;
        return(hitPrim);
    }
    else
    {
#ifdef RT_GRID_CACHE
        // ensure result is cached
        Stamp();
#endif
    }
                
    if (stats) stats->grids++;
    hitPrim = 0;
    
    // Figure out grid's parameters, and return if we completely miss its 
    // bounding box.
    if (!FindGridStart(start, direction, cell, increment, limit, tRayNext,
                       tDelta, tMax))
        return(hitPrim);
    
    // Traverse the grid & keep track of the closest intersection found so far.
    // tRayLast and tRayNext bracket the t values of the part of the ray that 
    // passes through the current cell. We're only interested in hits that
    // occur within these values.
    tRayLast = Max(sEpsilon, tRayNext);

    for (;;)
    {
        // We have to check all the triangles within this cell, since they
        // are in no particular order with respect to the ray.  
        if (stats) stats->voxels++;

        entry = grid[cell[0] * xSpan + cell[1] * ySpan + cell[2]];
        nextCellDir = MinEltIndex(tMax);    // the next cell after this...
        tRayNext = tMax[nextCellDir];       // t at which we leave this cell

        // Update record of delta t to the next cell in each direction.
        cell[nextCellDir] += increment[nextCellDir];
        tMax[nextCellDir] += tDelta[nextCellDir];
        
        // we're only interested in hits before tMin
        if (cell[nextCellDir] == limit[nextCellDir])
            tMin = vl_inf;
        else
            tMin = tRayNext;

        DBG_COUT << "searching " << tRayLast << " - " << tRayNext << endl;

        // Iterate through the primitives in this cell...
        for (; entry; entry = entry->next)
        {
            DBG_COUT << "checking cell prim" << endl;
            if (entry->next == PRIM_GRID_TAG)
            {
                // If this is a subgrid entry 
                if (entry->grid->ClosestIntersection(start, direction, 
                    origPrim, stats))
                {
                    t = entry->grid->hitT;

                    if (FEQG(tRayLast, t) && FEQG(t, tMin))
                    {
                        hitPrim = entry->grid->hitPrim;
                        hitT = t;
                        tMin = t;
                        DBG_COUT << "hit grid at " << t << endl;
                    }
                }
                break;  // RT_Grid is always last.
            }

            if (entry->prim == origPrim)
                continue;           
            if ((entry->prim->flags & TRI_INACTIVE) || !(sTag & entry->prim->object->tag))
                continue;

            DBG_COUT << "intersecting" << endl;
#ifdef RT_SUPPORT_GENS
            if (entry->prim->primType == rt_gen)
            {
                RT_GenPtr gen = *(RT_GenPtr*) entry->prim;
                                
                if (!gen->IsStamped())
                    gen->Intersect(start, direction, tRayLast, tMin);
                
                if (gen->IsTrue(GEN_HIT))
                {
                    tMin = gen->hitT;
                    hitT = tMin;
                    t = hitT;
                    hitPrim = gen;
                    if (stats) stats->hits++;
                    DBG_COUT << "hit gen at " << hitT << endl;
                }
                
            }
            else if (entry->prim->primType == rt_triangle)
#else
            if(1)
#endif
            {           
                RT_Tri*     thisTri = (RT_Tri*) entry->prim;
                    
                if (stats) stats->polysTested++;
                
                // Find ray intersection with tri's plane.
                // Bail if there isn't one.
                GCLReal normDotDirection = dot(thisTri->normal, direction);
    
                if (thisTri->flags & TRI_2SIDED_R)
                {
                    if (normDotDirection == 0.0)
                        continue;
                }
                else
                    if (normDotDirection >= 0.0)
                        continue;

                t = (thisTri->d + dot(thisTri->normal, start)) / 
                    -normDotDirection;
                    
                // Bail if it occurs before start of the cell, or after tMin
                if ((!FEQG(tRayLast, t) || !FEQG(t, tMin)))
                {
                    DBG_COUT << "hit outside cell" << endl;
                    continue;
                }
                // Okay, lets do the actual test. If we've hit the tri itself,
                // we can return a hit after we've finished this voxel.
                if (stats) stats->intersections++;
                reflectPoint = start + t * direction;
                
                if (!thisTri->IsStamped())
                    thisTri->PointIsInside(reflectPoint);
                
                if (thisTri->IsTrue(TRI_HIT))
                {
                    tMin = t;       
                    hitT = t;
#ifdef RT_GRID_SHADE_DEBUG
                    thisTri->id = level;
#endif
                    hitPrim = thisTri;  
                    if (stats) stats->hits++;
                    DBG_COUT << "hit tri at " << t << endl;
                }
            }
            else
                DBG_COUT << "unknown prim type: " << entry->prim->primType 
                         << endl;
        }

        if (hitPrim)
        {
            // We hit something inside this cell. Mission over, let's go home...
            if (stats) stats->goodHits++;
            return(hitPrim);
        }

        // If we've exited the grid, then bail--we didn't find any
        // intersections.

        if (cell[nextCellDir] == limit[nextCellDir])
            return(hitPrim);

        tRayLast = tRayNext;
    }
}


// --- Nested grid setup ------------------------------------------------------


Void RT_Grid::CreateNestedGrids()
// Finds cells with too many primitives in them, and creates new grids to hold
// them.
{
    Int             x, y, z, count;
    PrimEntry       *entry, *startEntry, **pushStart;
    RT_Grid         *subGrid, *bestGrid;
    PrimListEntry   *newPrimListEntry;
    Vector          cellMin, cellMax, bmin, bmax;
    GCLReal         newBVol, expRatio;
    GCLReal         minRatio;
    Vector          cellSpan, primsMin, primsMax, primsSpan;

    // Loop over the entire grid.
    
    for (x = 0; x < numCells[0]; x++)
        for (y = 0; y < numCells[1]; y++)
            for (z = 0; z < numCells[2]; z++)
            {
                count = 0;
                startEntry = grid[x * xSpan + y * ySpan + z];
                
                // + Don't push prims larger than areaLimit...
                // This limit exists for two reasons: if a primitive is very
                // large compared to its grid, it will create far too many 
                // prim entries, plus there's also not much gain in pushing
                // a primitive that's already on the order of the size 
                // of a voxel.
                
                // find start of pushable prims (the entries are sorted
                // by area)
                pushStart = &(grid[x * xSpan + y * ySpan + z]);
                for (entry = startEntry; entry; entry = entry->next)
                    if (entry->prim->area < areaLimit)
                        break;
                    else
                        pushStart = &(entry->next);

                // Count pushable prims
                for (entry = *pushStart; entry; entry = entry->next)
                    count++;
                                    
                // Add this cell's prims to a subgrid if there are too many of
                // them
                if (count > sMaxCellPrims * (level + 1))
                {
                    // Find the bounds of the prims

                    cellMin = min + Vector(x,  y,  z) * cellSize;
                    cellMax.MakeBlock(cellSize);
                    cellMax += cellMin;
                    primsMin.MakeBlock(+vl_inf);
                    primsMax.MakeBlock(-vl_inf);
                                        
                    for (entry = *pushStart; entry; entry = entry->next)
                    {
                        if (entry->prim->primType == rt_triangle)
                            ((RT_Tri*) entry->prim)->
                                UpdateBounds(primsMin, primsMax);
                        else if (entry->prim->primType == rt_gen)
                            (*(RT_GenPtr*) entry->prim)->
                                UpdateBounds(primsMin, primsMax);
                        else
                            DBG_COUT << "illegal prim type" << endl;
                    }
                    
                    // Trim the bounds to the boundaries of the cell

                    FindMaxElts(primsMin, cellMin, primsMin);
                    FindMinElts(primsMax, cellMax, primsMax);
                    
                    // We need a new grid. See whether we can merge with one of
                    // the already-created subgrids.
                    
                    minRatio = sMaxExpand;
                    bestGrid = 0;
                    primsSpan = primsMax - primsMin;
                    newBVol = primsSpan[0] * primsSpan[1] * primsSpan[2];
                    
                    for (subGrid = subGrids; subGrid; subGrid = subGrid->next)
                    {
                        GCLReal cellRatio;
                        
                        // Find the bounds of the putatative new grid.
                                                    
                        FindMinElts(subGrid->min, primsMin, bmin);
                        FindMaxElts(subGrid->max, primsMax, bmax);

                        // We want to ensure the cellSize of subgrids doesn't 
                        // get too large... so we find the ratio of the cell 
                        // size of the new grid to the parent (cellRatio)
                        // and the increase in volume over the old grid.
                        
                        cellRatio = MaxElt(bmax - bmin) / 
                            (sMaxCells * cellSize); 
                        expRatio = BoxVol(bmin,  bmax) /
                            (BoxVol(subGrid->min,  subGrid->max) + newBVol);

                        // If both are acceptable, this is our new best grid

                        if (expRatio < minRatio && cellRatio < sMaxCellRatio)
                        {
                            minRatio = expRatio;
                            bestGrid = subGrid;
                        }
                    }
                    
                    if (!bestGrid)
                    {
                        // We need a brand-spanking-new grid if we didn't 
                        // find one to merge with.
                    
                        bestGrid = RT_Alloc(RT_Grid);
                        
                        bestGrid->Init();
                        bestGrid->parent = this;
                        bestGrid->level = level + 1;
                        bestGrid->next = subGrids;
                        bestGrid->dirty = true; 
                        subGrids = bestGrid;
                    }
                    
                    // Add a pointer to the pushable prims to the prim list of
                    // the subgrid
                    
                    newPrimListEntry = RT_Alloc(PrimListEntry);
                    sMem += sizeof(PrimListEntry);
                    
                    newPrimListEntry->list = pushStart;
                    newPrimListEntry->next = bestGrid->addPrims;
                    bestGrid->addPrims = newPrimListEntry;
                                        
                    bestGrid->UpdateBounds(primsMin, primsMax);
                    bestGrid->totalPrims += count;
                }   
            }
            
    CullSubGrids();
    PushPrimitives();
}

Void RT_Grid::CullSubGrids()
// Cull subgrids that aren't worth the extra space.
{
    GCLReal     r;
    Int         estCellTotal;
    Vector      span;
    RT_Grid*    subGrid;
    RT_Grid**   last;
        
    last = &subGrids;
    
    for (subGrid = *last; subGrid; subGrid = *last)
    {
        span = subGrid->max - subGrid->min;
        r = sMaxCells / MaxElt(span);
        estCellTotal = Int(r * r * r * GCLReal(span[0] * span[1] * span[2]));
        
        if (subGrid->totalPrims / GCLReal(estCellTotal) < 
            sMinDensity * (1 - level / 2.0))
        {
            // Delete the subgrid from the list.
            
            *last = subGrid->next;
            subGrid->Free();
            RT_Free(subGrid);
        }
        else
            last = &(subGrid->next);
    }
}

Void RT_Grid::PushPrimitives()
// Shift primitives down to the grids created by CreateNestedGrids routine...
{
    PrimEntry       *entry, *nextEntry;
    RT_Grid         *subGrid;
    PrimListEntry   *primList;
    
    for (subGrid = subGrids; subGrid; subGrid = subGrid->next)
    {
        // Create the grids!
        subGrid->SetupGrid();
        subGrid->dirty = true;  
        RT_Prim::IncTime();
                        
        // Push the primitives!
        for (primList = subGrid->addPrims; primList; primList = primList->next)
        {
            // Each primList is a list of primitives pushed to this grid
            // from a particular cell.
            
            for (entry = *(primList->list); entry; entry = nextEntry)
            {
                nextEntry = entry->next;
                
                // Add primitives, avoiding duplicates
                
                if (!entry->prim->IsStamped())
                {
                    if (entry->prim->primType == rt_triangle)
                        subGrid->AddTriangle((RT_Tri*) entry->prim);
                    else if (entry->prim->primType == rt_gen)
                        subGrid->AddGeneric(*(RT_GenPtr*) entry->prim);
                    else
                        DBG_COUT << "illegal prim type" << endl;
                        
                    entry->prim->Stamp();
                }
                else
                    totalPrims--;

                RT_Free(entry);
                sMem -= sizeof(PrimEntry);
            }
                        
            // Remove the pushed primitives from the cell list, and replace
            // them with a grid entry.

            *primList->list = RT_Alloc(PrimEntry);
            sMem += sizeof(PrimEntry);
            (*primList->list)->next = PRIM_GRID_TAG;
            (*primList->list)->grid = subGrid;
        }
        
        addPrims = 0;
    }
}


// --- Support routines -------------------------------------------------------


static Void Indent(ostream &s, Int n)
{
    Int i;
    
    for (i = 0; i < n; i++)
        s << ' ';
}

Void RT_Grid::HierPrint(ostream &s, Int indent)
{
    RT_Grid*    subGrid;
    
    Indent(s, indent);
    s << totalPrims << " prims, bbox: " << min << ' ' << max << endl;
    
    for(subGrid = subGrids; subGrid; subGrid = subGrid->next)
        subGrid->HierPrint(s, indent + 2);
}


// --- Memory methods ---------------------------------------------------------


Void RT_Grid::SetMemLimit(Int maxKbs)
{
    sMemLimit = 1024 * maxKbs;
}

GCLReal RT_Grid::GridMemoryUsage()
{
    return(sMem / 1024.0);
}

GCLReal RT_Grid::MemoryUsage()
{
    return((sMem + sTriMem + sGenMem + sPointMem) / 1024.0);
}

Void RT_Grid::DumpMemoryUsage()
{
    RT_Grid     *sgrid;

    Int numGrids = 0;
    for (sgrid = subGrids; sgrid; sgrid = sgrid->next)
        numGrids++;
    
    cout << "Sub grids:           " << numGrids << endl;
    cout << "Grid pages:          " << grid.NumPages() << endl;
    cout << "Grid memory usage:   " << GridMemoryUsage() << "k" << endl;
    cout << "Triangle memory:     " << sTriMem / 1024.0 << "k" << endl;
    cout << "Generic memory:      " << sGenMem / 1024.0 << "k" << endl;
    cout << "Points memory:       " << sPointMem / 1024.0 << "k" << endl;
}

#ifdef RT_VIS

// --- Visualisation ----------------------------------------------------------

#include "gcl/Draw.h"

Void RT_Grid::Draw(Renderer &r, Int level)
{
    r.C(Colour4(HSVCol(hsvYellow + 60 * level, 1, 1), 0.2));
    FrameBox(r, min, max);
    r.C(Colour4(HSVCol(hsvOrange + 60 * level, 1, 1), 0.2));
    
    Vector  d1, d2;
    Point   p;
    Int     x, y, z;
    
    d1 = max - min;
    d2 = d1;
    d2[0] /= numCells[0];
    d2[1] /= numCells[1];
    d2[2] /= numCells[2];

    for (x = 0; x < numCells[0]; x++)
        for (y = 0; y < numCells[1]; y++)
            for (z = 0; z < numCells[2]; z++)
            {
                if (grid[x * xSpan + y * ySpan + z])
                {
                    p = min + Vector(x, y, z) * d2;
                    FrameBox(r, p, p + d2);
                }
            }

    // Draw subgrids.
    RT_Grid *sgrid;
    level++;

    for (sgrid = subGrids; sgrid; sgrid = sgrid->next)
        sgrid->Draw(r, level);
}

#endif

---