---

RadMesh.cc

Go to the documentation of this file.

/*
    File:           RadMesh.cc

    Function:       See header file

    Author(s):      Andrew Willmott

    Copyright:      (c) 1997-2000, Andrew Willmott

    Notes:          A RadElem can represent either a 4-sided parallelogram,
                    or a 3-sided triangle.

*/

#include "RadMesh.h"
#include "RadMethod.h"
#include "cl/String.h"
#include "gcl/VecUtil.h"
#include "Samples.h"

#ifdef RAD_VIS
#include "gcl/Forms.h"
#include "gcl/ScenePane.h"
#endif

// --- PatchStats class -------------------------------------------------------

Void PatchStats::Init()
{
    maxArea = 0.0;
    minArea = vl_inf;
    maxRefl = vl_0;
    minRefl.MakeBlock(vl_inf);
    maxRefl = vl_0;
    maxPower = vl_0;
    totalPower = vl_0;
}

Void PatchStats::Update(RadElem *patch)
{
    Colour  powerClr;

    FindMinCmpts(minRefl, patch->props->reflectance, minRefl);
    FindMaxCmpts(maxRefl, patch->props->reflectance, maxRefl);
    maxArea = Max(patch->area, maxArea);    
    minArea = Min(patch->area, maxArea);    

    powerClr = patch->props->emittance * patch->area;
    totalPower += powerClr;
    FindMaxCmpts(maxPower, powerClr, maxPower);
}

Void PatchStats::Report(ostream &s)
{
    s << "Patch statistics" << endl;
    s << "----------------" << endl;
    s << "Area range:        [" << minArea << ", " << maxArea << "]" << endl;
    s << "Reflectance range: [" << minRefl << ", " << maxRefl << "]" << endl;
    s << "Max. Power:         " << maxPower << endl;
    s << "Total Power:        " << totalPower << endl;
    s << endl;
}

// --- RadElem class ----------------------------------------------------------


Int RadElem::sGridMem = 0;
Int RadElem::sGridChildMem = 0;

RadElem::RadElem() : props(0),
#ifdef RAD_VIS
    highlight(0),
#endif
    colour(cPurple)
{
}

RadElem::~RadElem()
{
}

NbRadElem::NbRadElem() : RadElem()
{
    Int     i;
    
    for (i = 0; i < 4; i++)
    {
        nbFace[i] = 0;
        nbEdge[i] = 0xFF;       // for debugging only
    }
}

Void RadElem::DrawHighlight(Renderer &r)
{
#ifdef RAD_VIS
    if (highlight == 1)
        r.C(cRed);
    else if (highlight == 2)
        r.C(cYellow);
    else if (highlight == 3)
        r.C(cGreen);
    
    r.Begin(renPoly);
    SendPoints(r);
    r.End();
#endif  
}

Void RadElem::Draw(Renderer &r)
{
#ifdef RAD_VIS
    if (highlight)
    {
        DrawHighlight(r);
        if (!gRadControl->funcView)
            return;
    }   
#endif

    if (gRadControl->wire)
    {
        if (gRadControl->redWire)
            r.Begin(renLineLoop).SetColour(0.5 * (cRed + colour));
        else
            r.Begin(renLineLoop).SetColour(Mix(cWhite, colour,  0.8));
        
        SendPoints(r);
        r.End();

        return;
    }

    if (gRadControl->funcView)
    {
        r.Begin(renLineLoop).C(cWhite);
        SendPoints(r);
        r.End();
    }
    
    r.Begin(renPoly);

    if (!gRadControl->funcView)
    {
        if (props->texCoords && gRadControl->texture)
        {
            r.T(TexCoord(0)).C(VtxClr(0)).P(Vertex(0));
            r.T(TexCoord(1)).C(VtxClr(1)).P(Vertex(1));
            r.T(TexCoord(2)).C(VtxClr(2)).P(Vertex(2));
            
            if (IsQuad())
                r.T(TexCoord(3)).C(VtxClr(3)).P(Vertex(3));
        }
        else if (gRadControl->gouraud)
        {
            r.C(VtxClr(0)).P(Vertex(0)).C(VtxClr(1)).P(Vertex(1))
                .C(VtxClr(2)).P(Vertex(2));
            if (IsQuad())
                r.C(VtxClr(3)).P(Vertex(3));
        }
        else
        {
            r.C(colour);
            SendPoints(r);
        }
    }
    else
    {
        RaiseVertex(r, 0);
        RaiseVertex(r, 1);
        RaiseVertex(r, 2);
        if (IsQuad())
            RaiseVertex(r, 3);
    }

    r.End();
}

RadElem *RadElem::FindContainer(Coord &coord)
{
    return(this);
}

Colour RadElem::SampleLeaf(Coord c)
{
    if (gRadControl->gouraud)
    {   
        if (IsTri())
        {
            Colour      result;
            ClrReal     s = c[0];
            ClrReal     t = c[1];
            ClrReal     u = 1 - s - t;
            
            result = VtxClr(0) * c[1] + VtxClr(1) * u + VtxClr(2) * c[0];
            return(result);
        }
        else
        {
            Colour      c1, c2, result;
            ClrReal     s = c[0];
            ClrReal     t = c[1];
            
            c1 = (1 - s) * VtxClr(1) + s * VtxClr(2); 
            c2 = (1 - s) * VtxClr(0) + s * VtxClr(3); 
            result = (1 - t) * c1 + t * c2;
            return(result);
        }
    }
    else
        return(colour);
}

Colour RadElem::Sample(Coord c)
{
    RadElem *leaf = FindContainer(c);
    
    if (leaf)
        return(leaf->SampleLeaf(c));
    else
        return(cBlack);
}

Void RadElem::Compare(RadElem *to, GCLReal edgeLen, CompareStats &stats)
{
    // Sample density is controlled by edgelen: sample every edgeLen
    // units in either dimension. Same deal as with patchSubdivs.
    
    Colour  err, c1, c2;
    Coord   b,c;
    Colour  pcSum, pcSqrSum, pRefSum;
    Int     nx, ny, i, j;
    GCLReal xlen, ylen, weight;
    
    // calculate # samples per axis

    xlen = len(Vertex(2) - Vertex(1));
    ylen = len(Vertex(0) - Vertex(1));
    
    nx = (Int) (edgeLen * xlen);
    ny = (Int) (edgeLen * ylen);

    if (nx < 1) nx = 1;
    if (ny < 1) ny = 1;

    pcSum = vl_0;
    pRefSum = vl_0;
    pcSqrSum = vl_0;
    
    // Sample away!!!!
    
    for (i = 0; i < nx; i++)
        for (j = 0; j < ny; j++)
        {
            
            c[0] = (0.5 + i) / GCLReal(nx);
            c[1] = (0.5 + j) / GCLReal(ny);
            if (IsTri())
                ; // XXX need to resample, handle aniso for tri
            
            c1 = Sample(c);
            c2 = to->Sample(c);
    
            err = c1 - c2;                  // colour difference
            pcSum += err;                   // avg. c diff
            pcSqrSum += sqr(err);           // c variance
            pRefSum += (c1 - Emittance());  // reference: c1 - E.
        }
    
    weight = area / GCLReal(nx * ny);
    pcSum *= weight;
    pRefSum *= weight;
    pcSqrSum *= weight;

    // update stats, weighted by area-per-sample.

    stats.areaSum += area;
    stats.cSum += pcSum;
    stats.cSqrSum += pcSqrSum;
    stats.refSum += pRefSum;
}

Void RadElem::SetHighlight(Int h)
{
#ifdef RAD_VIS
    highlight = h;
#endif
}

StrConst RadElem::Name()
{
    return("rad");
}

Void RadElem::Print(ostream &s)
{
    s << index[0] << ' ' << index[1] << ' ' << index[2] << ' ' 
      << index[3] << ' ';
    s << clrIdx[0] << ' ' << clrIdx[1] << ' ' << clrIdx[2] << ' '
      << clrIdx[3] << ' ';
    s << colour << ' ';         
}

Void RadElem::Parse(istream &s)
{
    s >> index[0] >> index[1] >> index[2] >> index[3];
    s >> clrIdx[0] >> clrIdx[1] >> clrIdx[2] >> clrIdx[3];
    s >> colour;
}


// --- Form factor routines ---------------------------------------------------


GCLReal RadElem::EstPatchFactor(const Point &p, const Vector &np)
//  Elemental factor from this patch to p. 
//  essentially cos theta * cos phi * Ap / (pi * r^2)
{
    GCLReal result, rad4;
    Vector  temp;

    temp = p - Centre();
    
    result = dot(Normal(), temp);
    
    if (result <= 0)
        result = 0;
    else 
    {
        result *= -dot(np, temp);
        
        if (result <= 0)
            result = 0;
        else
        {
            rad4 = sqr(sqrlen(temp));
            result *= area;
            result /= (rad4 * vl_pi);
        }
    }
    
    return(result);
}


GCLReal RadElem::EdgeArea(const Vector &p, const Vector &q, const Vector &n)
{ 
    // Numerically stable edge-area calculator. (Including degenerate cases.)
    // Andrew J. Willmott, 1992

    // 4 dots, 1 cross, 1 atan per edge.

    const GCLReal epsilon = 1e-6;
    GCLReal sinFactor, qdotp, projArea;

    qdotp =  dot(p, q);
    sinFactor = sqrlen(p) * sqrlen(q) - sqr(qdotp);

    projArea = dot(-cross(p, q), n);

    if (sinFactor > 0)              
        // If SinFactor > 0 we may apply the formula without problem 
    {
        sinFactor = sqrt(sinFactor);
        return(projArea * (vl_halfPi - atan(qdotp / sinFactor)) / 
               (2 * vl_pi * sinFactor));
    }   
        // If not, we have three remaining cases... 
    else if (qdotp > epsilon)
        // CASE 1: p and q point in the same direction }
        return(0);
    else if (qdotp < -epsilon)
        // CASE 2: Viewpoint lies on the edge }
        return(0.5);
    else
        // CASE 3: Viewpoint lies at either end of the edge }
        return(0.125);
}

GCLReal RadElem::PatchFactor(const Point &p, const Vector &np)
// Common terminology for this is the polygon-to-point form factor. 
// The 'form-factor' mentioned in radiosity papers is the average
// patch factor across a patch.
{
    Vector  x1, x2, x3;
    GCLReal sum;
    
    x1 = Vertex(0);
    x1 -= p;
    x3 = x1;
    x2 = Vertex(1);
    x2 -= p;
    
    sum = EdgeArea(x1, x2, np);
    x1 = Vertex(2);
    x1 -= p;
    sum += EdgeArea(x2, x1, np);
    
    if (IsTri())
        sum += EdgeArea(x1, x3, np);
    else
    {
        x2 = Vertex(3);
        x2 -= p;
        sum += EdgeArea(x1, x2, np);
        sum += EdgeArea(x2, x3, np);
    }

    if (sum > 0.0)
        return(sum);
    else
        return(0.0);
}

GCLReal RadElem::ApproxPatchFactor(const Point &p, const Vector &np)
// Uses point-to-point PF as default, and swaps to poly-to-point PF
// when necessary.
{
    GCLReal result, rad4, npt, res2;
    Vector  temp;

    temp = p - Centre();
    
    result = -dot(np, temp);
    
    if (result <= 0)
        result = 0;
    else 
    {
        npt = dot(Normal(), temp);
        
        if (npt <= 0)
            result = 0;
        else
        {
            rad4 = sqr(sqrlen(temp));
            result *= area / vl_pi;
            res2 = result * Max(npt, 0.25 * sqrt(area));

            if (gRadControl->quadLevel > 0 && rad4 * gRadControl->dFError < res2)
            {
                // result is blowing up!
                // let's bail to the patch factor formula
                return(PatchFactor(p, np));
            }
            else
                result *= npt / rad4;
        }
    }
    
    return(result);
}

GCLReal RadElem::EstFormFactor(RadElem *to)
{
    GCLReal result, rad4, npt, res2;
    Vector  temp;

    temp = to->Centre() - Centre();
    
    result = -dot(to->Normal(), temp);
    
    if (result <= 0)
        result = 0;
    else 
    {
        npt = dot(Normal(), temp);
        
        if (npt <= 0)
            result = 0;
        else
        {
            rad4 = sqr(sqrlen(temp));
            result *= area / vl_pi;
            res2 = result * Max(npt, 0.25 * sqrt(area));

            if (gRadControl->quadLevel == 0 || rad4 * gRadControl->dFError >= res2)
                result *= npt / rad4;
            else
            {
                // result is blowing up!
                // let's bail to the patch factor formula

                GCLReal sum;

                if (gRadControl->quadLevel > 1)
                {
                    // use more-accurate PF estimate:
                    // sample corners
                    
                    sum = PatchFactor(to->Vertex(0), to->Normal());
                    sum += PatchFactor(to->Vertex(1), to->Normal());
                    sum += PatchFactor(to->Vertex(2), to->Normal());
                    sum += PatchFactor(to->Vertex(3), to->Normal());
                    
                    if (gRadControl->quadLevel > 2)
                    {
                        // estimate ff from tetrahedron constructed from
                        // corner samples and centre. this will
                        // underestimate the true ff slightly.
                        
                        sum += 2.0 * PatchFactor(to->Centre(), to->Normal());
                        sum /= 6.0;
                    }
                    else
                    // underestimate curve
                        sum /= 4.0;
                }
                else
                    sum = PatchFactor(to->Centre(), to->Normal());
                
                return(sum);
            }
        }
    }
    
    return(result);
}

GCLReal RadElem::SampledFormFactor(Int n, RadElem *to, GCLReal &error)
{
    GCLReal sum, sample, min, max, rn;
    Point   destPt;
    Vector  dx =     to->Vertex(2) - to->Vertex(1);
    Vector  dy =     to->Vertex(0) - to->Vertex(1);
    Int     i, j;
    
    // use most-accurate PF estimate:
    // sample in a grid
    
    min = vl_inf;
    max = -vl_inf;
    sum = 0;
    rn = n;
    
    for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
        {
            destPt = to->Vertex(1);
            destPt += ((j + 0.5) / rn) * dx;
            destPt += ((i + 0.5) / rn) * dy;
        
            sample = ApproxPatchFactor(destPt, to->Normal());
            sum += sample;
            min = Min(min, sample);
            max = Max(max, sample);
        }
    
    error = max - min;

    return(sum / sqr(n));
}

Int RadElem::OrientInfo(RadElem *to)
// Returns basic orientation info: 0 = centres of
// patches are obscured, -1 = patches close enough
// to warrant good PF approx, 1 = can use standard approx.
{
    GCLReal result, rad4, npt, res2;
    Vector  temp;

    temp = to->Centre() - Centre();
    result = -dot(to->Normal(), temp);
    
    if (result <= 0)
        return(0);
    else 
    {
        npt = dot(Normal(), temp);
        
        if (npt <= 0)
            return(0);
        else
        {
            rad4 = sqr(sqrlen(temp));
            result *= area / vl_pi;
            res2 = result * Max(npt, 0.25 * sqrt(area));

            if (gRadControl->quadLevel > 0 &&
                rad4 * gRadControl->dFError < res2)
                return(-1);
            else
                return(1);
        }
    }
}

Void RadElem::SetProps(RadProps *parentProps)
{
    props = parentProps;
    Centre().MakeZero();
    
    centre += Vertex(0);
    centre += Vertex(1);
    centre += Vertex(2);
    if (IsQuad())
    {
        centre += Vertex(3);
        centre /= 4;
    }
    else
        centre /= 3;

    CalcTriAreaNormal(Vertex(0), Vertex(1), Vertex(2), normal);
    area = len(normal);
    if (area > 0)
        normal /= area;
    if (IsTri())
        area *= 0.5;

#ifdef RAD_TEXTURE
    if (gRadControl->textureRefl && IsTextured())
        localRefl = props->texture->FilterBox(TexCoord(0), TexCoord(2));
    else
        localRefl = props->reflectance;
#endif
}

Void RadElem::Reanimate(RadElem *parent)
{
    if (parent)
        SetProps(parent->props);
}

Void RadElem::ColourVertices(Int weights[])
{
    Int j;
    
    for (j = 0; j < Sides(); j++)
    {
        VtxClr(j) += colour;
        weights[clrIdx[j]] += 1;
    }   
}

Void RadElem::CreatePatches(PatchList &patches)
{
    Assert(false, "(RadElem::CreatePatches) need grid for that.");
}

Void RadElem::RaiseVertex(Renderer &r, Int i)
{
    Colour *c;
    
    if (gRadControl->gouraud)
        c = &VtxClr(i);
    else
        c = &colour;

    r.C(*c);

    if (HasNormals())
        r.P(Vertex(i) + dot(kRadRGBToLum, *c) * Normal(i));
    else    
        r.P(Vertex(i) + dot(kRadRGBToLum, *c) * Normal());
}

ostream &operator << (ostream &s, RadElem &rq)
{
    rq.Print(s);
    return(s);
}

istream &operator >> (istream &s, RadElem &rq)
{
    rq.Parse(s);
    return(s);
}


// --- Inter-elem visibility --------------------------------------------------


// utility routines for patches only

GCLReal RadElem::Visibility16(const Point &p, const Vector &n)
//  Returns the fractional visibility of this patch from p.
{
    Vector  temp;
    GCLReal d;

    // Is 'p' completely behind this patch?
    
    d = dot(Normal(), Centre());
    
    if (dot(p, Normal()) < d)
        return(0.0);
    
    // Is this patch completely behind p's plane?
        
    d = dot(n, p);

    if ((dot(Vertex(0), n) <= d)
            && (dot(Vertex(1), n) <= d)
            && (dot(Vertex(2), n) <= d)
            && (IsTri() || dot(Vertex(3), n) <= d))
        return(0.0);
        
    //  If neither is the case, we have at least partial
    //  visibility, so ray cast to estimate it

    return(RadVis16x1(p, n));   
}

GCLReal RadElem::CentreVis(const Point &p, const Vector &n)
// Returns visibility of the centre of this patch from p.
{
    Vector  temp;
    Point   dstPoint, srcPoint;
    GCLReal d;
    Bool    hit;
    
    // Is 'p' completely behind this patch?
    
    d = dot(Normal(), Centre());
    
    if ((dot(p, Normal())) <= d)
        return(0.0);
    
    // Is the centre of this patch behind p's plane?
        
    d = dot(n, p);

    if (dot(Centre(), n) <= d)
        return(0.0);
        
    //  If neither is the case, we have at least partial
    //  visibility, so ray cast to estimate it

    srcPoint = p;
    srcPoint += n * kRaySurfEps;
    dstPoint = Centre();
    dstPoint += Normal() * kRaySurfEps;

    hit = gRadControl->radObject->IntersectsWithRay(srcPoint, dstPoint);
    gRadControl->rays += 1;

#ifdef RAD_VIS
    if (gRadControl->showRays)
    {
        if (hit)
            gRadControl->radObject->display->
                Begin(renLines).C(cRed).P(srcPoint).P(dstPoint).End();
        else
            gRadControl->radObject->display->
                Begin(renLines).C(cGreen).P(srcPoint).P(dstPoint).End();
    }
#endif

    if (hit)
        return(0.0);
    else
        return(1.0);
}

Bool RadElem::PotentiallyVis(RadElem *to)
// Checks to see whether any part of to is potentially visible
// from this patch, and vice versa. 
// returns true if there is potential visibility, false if not.
{
    Vector  temp;
    GCLReal dFrom, dTo;

    // Is 'to' completely behind this patch?
    
    dFrom = dot(Normal(), Vertex(0));
    
    if ((dot(to->Vertex(0), Normal()) <= dFrom)
            && (dot(to->Vertex(1), Normal()) <= dFrom)
            && (dot(to->Vertex(2), Normal()) <= dFrom)
            && (to->IsTri() || dot(to->Vertex(3), Normal()) <= dFrom))
        return(false);
    
    // Is this patch completely behind 'to'?
        
    dTo =   dot(to->Normal(), to->Vertex(0));

    if ((dot(Vertex(0), to->Normal()) <= dTo)
            && (dot(Vertex(1), to->Normal()) <= dTo)
            && (dot(Vertex(2), to->Normal()) <= dTo)
            && (IsTri() || dot(Vertex(3), to->Normal()) <= dTo))
        return(false);
        
    //  If neither is the case, we have at least potential
    //  visibility.

    return(true);
}

Bool RadElem::PotentiallyVisAndTouching(RadElem *to, Bool &touching)
// Checks to see whether any part of to is potentially visible
// from this patch, and vice versa. 
// returns true if there is potential visibility, false if not.
// Also sets touching to true if the two patches intersect or touch.
{
    Vector  temp;
    GCLReal dFrom, dTo, a;
    Bool    pvis, touch;
    Int     n, z, p;
    
    touching = false;
    dFrom = dot(Normal(), Vertex(0));
    n = p = z = 0;

    // we count p=points above plane, z=points on plane, n = points below plane.
    // pvis = at least 1 p1 coeff +ve && at least 1 p2 coeff +ve
    // touch: p1-touch && p2-touch
    // x-touch: x coeffs not all +ve or all -ve.
    // 
    
    // opt: can do early termination of touch as soon as z>0 or p&n
    //      can do early termination of pvis as soon as p>0
    //      can terminate both iff p>0 && (z>0 | n>0)
    //  currently don't do this 'cause extra testing might
    //  kill the optimisations, and I'm lazy...
    //
            
    a = dot(to->Vertex(0), Normal()) - dFrom;
    if (a > 0) p++; else if (a < 0) n++; else z++;
    a = dot(to->Vertex(1), Normal()) - dFrom;
    if (a > 0) p++; else if (a < 0) n++; else z++;
    a = dot(to->Vertex(2), Normal()) - dFrom;
    if (a > 0) p++; else if (a < 0) n++; else z++;
    if (to->IsQuad())
    {
        a = dot(to->Vertex(3), Normal()) - dFrom;
        if (a > 0) p++; else if (a < 0) n++; else z++;
    }
        
    pvis = p;
    touch = z || (p && n);
        
    if (!pvis && !touch)
        return(false);

    dTo =   dot(to->Normal(), to->Vertex(0));
    n = p = z = 0;

    a = dot(Vertex(0), to->Normal()) - dTo;
    if (a > 0) p++; else if (a < 0) n++; else z++;
    a = dot(Vertex(1), to->Normal()) - dTo;
    if (a > 0) p++; else if (a < 0) n++; else z++;
    a = dot(Vertex(2), to->Normal()) - dTo;
    if (a > 0) p++; else if (a < 0) n++; else z++;
    if (IsQuad())
    {
        a = dot(Vertex(3), to->Normal()) - dTo;
        if (a > 0) p++; else if (a < 0) n++; else z++;
    }
        
    if (touch && (z || (p && n)))
        touching = true;
    return(pvis & p);
}

GCLReal RadElem::Visibility44(RadElem *to)
{
    if (PotentiallyVis(to))
        return(RadVis4x4(to));
    else
        return(0.0);
}


// --- Generic element->element visibility method -----------------------------

GCLReal RadElem::Visibility(RadElem *to)
{
    GCLReal result;
    
#ifdef RAD_VIS
    if (gRadControl->showRays)
        RM_DISPLAY_START;
#endif

    switch(gRadControl->visibility)
    {
    case vis_none:
        result = 1.0;
        break;
    case vis_1:
        result = CentreVis(to->Centre(), to->Normal());
        break;
    case vis_16x1:
        result = Visibility16(to->Centre(), to->Normal());
        break;
    case vis_4x4:
        result = Visibility44(to);
        break;
    default:
        Assert(false, "Illegal visibility method specified!");
    }
    
#ifdef RAD_VIS
    if (gRadControl->showRays)
    {
        RM_DISPLAY_END;
        RM_OUT1("Visibility: " << result);
        if (RM_PAUSE) return(0);
    }
#endif

    return(result);
}


// ----------------------------------------------------------------------------

//  Arrays for jittered sampling; 4 x 4 source, 4 x 4 dest.

static Int gRotPick = 0;

/*  
    Routine: RadVis4x4
    Measures visibility between the patch and q.
    Returns visibility factor between 0 and 1. 
    Uses magic-square jittered 4x4 pattern, a la HR paper.
 */

Void RadElem::SetVisPoints(Point *pts)
{
    Int         i;
    Vector      dstHoriz, dstVert;
    Coord       *d1;
    
    dstVert =  Vertex(0) - Vertex(1);
    dstHoriz = Vertex(2) - Vertex(1);

    if (IsQuad())
        d1 = (Coord *) tSquareSamples16[gRotPick]; 
    else
        d1 = (Coord *) tTriangleSamples16[gRotPick];

    for (i = 0; i < 16; i++)
        pts[i] = Vertex(1);

    for (i = 0; i < 16; i++)
        pts[i] += d1[i][0] * dstHoriz +
            d1[i][1] * dstVert + kRaySurfEps * Normal();
    
    if (gRadControl->jitterRot)
        if (++gRotPick == kNumSampArrays)
            gRotPick = 0;
}

GCLReal RadElem::RadVis4x4(RadElem *to)
{
    Int         i, misses;
    Point       srcPoint[16], dstPoint[16];
    Bool        hit;
            
    misses = 0;
    SetVisPoints(srcPoint);
    to->SetVisPoints(dstPoint);
    
    gRadControl->rays += 16;
    
    for (i = 0; i < 16; i++)
    {
        hit = gRadControl->radObject->IntersectsWithRay(srcPoint[i], dstPoint[i]);
        
        if (hit)
            misses++;

#ifdef RAD_VIS
        if (gRadControl->showRays)
        {
            if (hit)
                gRadControl->radObject->display->Begin(renLines).C(cRed).P(srcPoint[i]).
                    P(dstPoint[i]).End();
            else
                gRadControl->radObject->display->Begin(renLines).C(cGreen).P(srcPoint[i]).
                    P(dstPoint[i]).End();
        }
#endif
    }
                
    return((16 - misses) * 1.0 / 16.0);
}

/*  
    Routine: RadVis16x1
    Measures visibility between the patch and point p by
    casting 16 rays between the two.
    Returns visibility factor between 0 and 1. 
    Uses magic-square jittered 4x4 pattern
 */

GCLReal RadElem::RadVis16x1(const Point &p, const Vector &n)
{
    Int         i, misses;
    Point       srcPoint[16], dstPoint;
    Bool        hit;
            
    SetVisPoints(srcPoint);
    dstPoint = p;
    dstPoint += kRaySurfEps * n;
    misses = 0;

    gRadControl->rays += 16;
    
    for (i = 0; i < 16; i++)
    {
        hit = gRadControl->radObject->IntersectsWithRay(srcPoint[i], dstPoint);
        
        if (hit)
            misses++;

#ifdef RAD_VIS
        if (gRadControl->showRays)
        {
            if (hit)
                gRadControl->radObject->display->
                    Begin(renLines).C(cRed).P(srcPoint[i]).P(dstPoint).End();
            else
                gRadControl->radObject->display->
                    Begin(renLines).C(cGreen).P(srcPoint[i]).P(dstPoint).End();
        }
#endif
    }
        
    return(GCLReal(16 - misses) / 16.0);
}

Coord RadElem::FindCoord(Point &p)
{
    Point       origin;
    Vector      r;
    Coord       result;
    VecTrans    s, t;
    
    // construct a coordinate system about vertex 1 of the triangle/quad
    origin = Vertex(1);
    t[2] = Normal();
    t[1] = Vertex(0) - origin;
    t[0] = Vertex(2) - origin;

    // find transformation into that coord system...
    s = inv(t);

    // transform p
    r = p - origin;
    r *= s;
    
    // drop height and return
    return(Coord(r[0], r[1]));
}

Void RadElem::SetIndexes(
                Int         dstIdx,       
                RadElem     *src,       
                Int         srcIdx
            )
{
    index[dstIdx]   = src->index[srcIdx];   
    clrIdx[dstIdx]  = src->clrIdx[srcIdx];  
    normIdx[dstIdx] = src->normIdx[srcIdx]; 
    texIdx[dstIdx]  = src->texIdx[srcIdx];  
}

Void RadElem::SetAllIndexes(
                RadElem     *s0,
                Int         i0,       
                RadElem     *s1,       
                Int         i1,       
                RadElem     *s2,       
                Int         i2,
                RadElem     *s3,
                Int         i3
            )
{
    index[0]   = s0->index[i0]; 
    index[1]   = s1->index[i1]; 
    index[2]   = s2->index[i2]; 
    index[3]   = s3->index[i3]; 

    clrIdx[0]  = s0->clrIdx[i0];    
    clrIdx[1]  = s1->clrIdx[i1];    
    clrIdx[2]  = s2->clrIdx[i2];    
    clrIdx[3]  = s3->clrIdx[i3];    

    normIdx[0] = s0->normIdx[i0];   
    normIdx[1] = s1->normIdx[i1];   
    normIdx[2] = s2->normIdx[i2];   
    normIdx[3] = s3->normIdx[i3];   

    texIdx[0]  = s0->texIdx[i0];    
    texIdx[1]  = s1->texIdx[i1];    
    texIdx[2]  = s2->texIdx[i2];    
    texIdx[3]  = s3->texIdx[i3];    
}

GCLReal RadElem::MemoryUse()
{
    return(sizeof(SELF));
}

GCLReal NbRadElem::MemoryUse()
{
    return(sizeof(SELF));
}


// --- Grid Subdivision -------------------------------------------------------

Void GridBase::Mesh(GCLReal density)
// Meshes an element into a grid of similar elements of roughly 
// the given density.
// XXX still need to re-fix non-lin density stuff
{
    Int         i, j, in, jn;
    RadElem     child, child2, rowEnds;
    GCLReal     w, h, dw, dh;
    Int         im, jm;
    RadElemPtr  *lastRow, upperTri;
    Vector      u, v;

    //  Setup increments, step sizes, etc.

    if (gRadControl->mesh == mesh_random)
    {
        density = (0.1 + 0.9 * drand48()) * density;
    }

    u = Vertex(2) - Vertex(1);
    v = Vertex(0) - Vertex(1);
    
    in = (Int) ceil(len(v) * density);
    jn = (Int) ceil(len(u) * density);
    
    // For triangles we must subdivide equally in both directions.
    
    if (IsTri())            
    {
        if (jn < in)
            jn = in;
        else
            in = jn;
    }
    
    if (gRadControl->mesh != mesh_nonlin)
    {
        im = 0;
        jm = 0;
        w = 1.0 / jn;
        h = 1.0 / in;
        child.area = len(w) * len(h);
        if (IsTri())
            child.area *= 0.5;
        child2.area = child.area;   
        // child2 is always a triangle: it is not used for quad meshing.
    }
    else
    {
        if ((jn % 2) == 1)
            w = 4.0 / sqr(jn + 1);
        else
            w = 4.0 / (sqr(jn + 1) - 1);
            
        if ((in % 2) == 1)
            h = 4.0 / sqr(in + 1);
        else
            h = 4.0 / (sqr(in + 1) - 1);

        jm = jn;
        im = in;
        dw = w;
        dh = h;
    }
        
    // Set up templates for new elements

    lastRow = new RadElemPtr[jn];   
    child.index[3] = -1;
    child2.index[3] = -1;
            
    // Add the new quads/tri's to the grid in turn, from left to right,
    // and top to bottom. Each new quad or pair of tris will require a
    // new bottom-right vertex.         

    for (i = 0; i < in; i++)        
    {                               
        // find outermost bottom vertices of the row, so we can interpolate
        // between them while filling in the row. rowEnds:0 is the left-most
        // bottom vertex, rowEnds:1 is the right-most.
        rowEnds.props = props;
        if (i == in - 1)                
        // Special-case original bottom-left and bottom-right vertices.
        {
            rowEnds.SetIndexes(0, this, 1); 
            rowEnds.SetIndexes(1, this, 2); 
        }
        else
        {
            props->InterpolateIndexes(this, 0, 1, &rowEnds, 0, 
                GCLReal(i + 1.0) / in);
            if (IsQuad())
                props->InterpolateIndexes(this, 3, 2, &rowEnds, 1, 
                    GCLReal(i + 1.0) / in);
            else
                props->InterpolateIndexes(this, 0, 2, &rowEnds, 1, 
                    GCLReal(i + 1.0) / in);
        }

        //  copy left-most into first quad/tri in the row
        child.SetIndexes(1, &rowEnds, 0);
        if (i == 0)
            child.SetIndexes(0, this, 0);
        else
            child.SetIndexes(0, lastRow[0], 1);
        
        if (jm > 0)
            w = dw;
        if (IsTri())
            jn = i + 1;
                
        for (j = 0; j < jn; j++)
        //  Add one row of quads/tri's                      
        {
            // For each quad/tri, we must add the bottom-right corner as a new point.
            if (j == jn - 1)        
            // Special-case right-most quad/tri
                child.SetIndexes(2, &rowEnds, 1);
            else
                props->InterpolateIndexes(&rowEnds, 0, 1, &child, 2, 
                    Double(j + 1.0) / jn);
            
            if (IsTri())
            {               
                //  For triangles: we add them in pairs. The first triangle...
                child.SetProps(props);
                if (jm > 0)
                {
                    child.area = 0.5 * len(w) * len(h);
                    child2.area = child.area;
                }
                upperTri = lastRow[j];
                lastRow[j] = AddChild(child, i, j, in, jn);
                
                if (j == i)         // if we're at the end of a triangle strip
                    break;          // break to the outer loop
                                
                //  then the second triangle.

                // copy upper-right vertex from last row
                child2.SetAllIndexes(&child, 2, upperTri, 2, &child, 0, this, 3);
                child2.SetProps(props);
                AddChild(child2, i, j, in, jn);

                // Move on to next quad/tri pair: shift the indices left.
                child.SetIndexes(0, &child2, 1);
                child.SetIndexes(1, &child, 2);
            }
            else
            {       
                //  For quads: we first find/copy upper-right vertex

                if (i == 0)         
                // first row -- can't borrow from row above
                {
                    if (j == jn - 1)
                    // Special-case original top-right corner
                    {
                        child.SetIndexes(3, this, 3);
                    }
                    else
                    {
                        // interpolate between upper-left & upper right.
                        props->InterpolateIndexes(this, 0, 3, &child, 3, 
                            Double(j + 1.0) / jn);
                    }
                }
                else
                {
                    // set to bottom-right vertex of quad above us
                    child.SetIndexes(3, lastRow[j], 2);
                }
    
                // And add the new quad...
    
                if (jm > 0)
                    child.area = len(w) * len(h);
                child.SetProps(props);
                lastRow[j] = AddChild(child, i, j, in, jn);             

                // Move on to next quad/tri pair: shift the indices left.
                child.SetIndexes(0, &child, 3);
                child.SetIndexes(1, &child, 2);
            }
            
            if (jm > 0)
            {
                if (2 * (j + 1) < jm)
                    w += dw;
                else if (2 * (j + 1) > jm)
                    w -= dw;
            }
        }

        // Move on to next row, resetting necessary variables
        
        if (im > 0)
        {
            if (2 * (i + 1) < im)
                h += dh;
            else if (2 * (i + 1) > im)
                h -= dh;
        }
    }
    
    rows = in;
    cols = jn;
    delete[] lastRow;
}

Int GridBase::FindChildIndex(Coord &coord)
// which of the grid elements does the coord fall inside?
// find its index, and modify coord to be in its local coord. system.
{
    Int     x, y, xy;
    GCLReal rx, ry;

    rx = coord[0];
    ry = coord[1];

    // Do some bounds checking...
    if (IsTri())
    {
        if ((rx - ry) > 1 || rx < 0 || ry < 0)
            return(-1);
    }
    else
        if (rx < 0 || ry < 0 || rx > 1 || ry > 1)
            return(-1);

    // find rectangular grid coords
    rx *= cols;
    ry *= rows;

    x = (Int) rx;
    if (x == cols)
        x--;
    y = (Int) ry;
    if (y == rows)
        y--;    

    // find cell coords
    // our grid is numbered from the top down, hence the (rows - 1 - y)
    // terms below

    if (IsTri())
    { 
        xy = (Int) (rows * (rx + ry));
        if (xy == rows)
            xy--;   

        if (xy > x + y)
        // flipped triangle
        {
            coord[0] = x + 1 - rx;
            coord[1] = y + 1 - ry;
            x = 2 * x + 1;
        }
        else
        // normal triangle
        {
            coord[0] = rx - x;
            coord[1] = ry - y;
            x = 2 * x;
        }

        // number of picked cell
        return(x + sqr(rows - 1 - y));
    }
    else
    {       
        // quad case is much easier!
        coord[0] = rx - x;
        coord[1] = ry - y;
        
        return(x + (rows - 1 - y) * cols);
    }
}

// --- RadGrid methods ---------------------------------------------------------

RadElemPtr RadGrid::AddChild(RadElem &elem, Int i, Int j, Int in, Int jn)
{
    children.Append(elem);

    return(&children.Last());
}

Void RadGrid::Draw(Renderer &r)
{   
    if (children.NumItems() > 0
#ifdef RAD_VIS
         && !highlight
#endif
     )
    {
        Int i;
        
        for (i = 0; i < children.NumItems(); i++)
            children[i].Draw(r);
    }
    else
        RadElem::Draw(r);
}

Void RadGrid::CreatePatches(PatchList &patches)
{
    Int     i;
    
    if (gRadControl->patchSubdivs == 0.0)
    {
        patches.Append(this);
        rows = cols = 1;
        return;
    }

    children.Clear();
    Mesh(gRadControl->patchSubdivs);    

    for (i = 0; i < children.NumItems(); i++)
        patches.Append(&(children[i]));
}

RadElem *RadGrid::FindContainer(Coord &coord)
{
    Int i = FindChildIndex(coord);
    
    return(children[i].FindContainer(coord));
}

Void RadGrid::ColourVertices(Int weights[])
{
    Int i;
    
    for (i = 0; i < children.NumItems(); i++)
        children[i].ColourVertices(weights);
}

StrConst RadGrid::Name()
{
    return("grid");
}

Void RadGrid::Print(ostream &s)
{
    RadElem::Print(s);
    s << rows << ' ' << cols << ' ' << children << ' ';
}

Void RadGrid::Parse(istream &s)
{
    Char c;
    Int i;
    
    RadElem::Parse(s);
    
    s >> rows >> cols;
    
    children.SetSize(rows * cols);

    ChompWhiteSpace(s); 
    s.get(c);
    for (i = 0; i < children.NumItems(); i++)
        children[i].Parse(s);
    ChompWhiteSpace(s);
    s.get(c);
}

Void RadGrid::Reanimate(RadElem *parent)
{
    Int i;

    RadElem::Reanimate(parent);

    for (i = 0; i < children.NumItems(); i++)
        children[i].Reanimate(this);
}

GCLReal RadGrid::MemoryUse()
{
    sGridMem += sizeof(SELF);
    sGridChildMem += sizeof(RadElem) * children.NumItems();
    return(sizeof(SELF) + sizeof(RadElem) * children.NumItems());
}


// --- RadProps methods -------------------------------------------------------


RadProps::RadProps() : points(0), colours(0), normals(0), texCoords(0),
    texture(0)
{
}

Void RadProps::InterpolateIndexes(
                RadElem *srcElt, 
                Int     srcVtx1, 
                Int     srcVtx2, 
                RadElem *dstElt, 
                Int     dstVtx,
                GCLReal m
            )
{
    GCLReal     mc = 1.0 - m;

    // do interpolation
    
    points->Append(srcElt->Vertex(srcVtx1) * mc + srcElt->Vertex(srcVtx2) * m);
    dstElt->index[dstVtx] = points->NumItems() - 1;

    colours->Append(srcElt->VtxClr(srcVtx1) * mc + srcElt->VtxClr(srcVtx2) * m);
    dstElt->clrIdx[dstVtx] = colours->NumItems() - 1;

    if (normals)
    {
        normals->Append(srcElt->Normal(srcVtx1) * mc + srcElt->Normal(srcVtx2) * m);
        dstElt->normIdx[dstVtx] = normals->NumItems() - 1;
    }
    
    if (texCoords)
    {
        texCoords->Append(srcElt->TexCoord(srcVtx1) * mc
            + srcElt->TexCoord(srcVtx2) * m);
        dstElt->texIdx[dstVtx] = texCoords->NumItems() - 1;
    }
}

---