---

HierRad.cc

Go to the documentation of this file.

/*
    File:           HierRad.cc

    Function:       See header file

    Author(s):      Andrew Willmott

    Copyright:      (c) 1997-2000, Andrew Willmott

    Notes:          
*/

#include "HierRad.h"
#include "GenLink.h"
#include "Haar.h"
#include <stdio.h>
#include <unistd.h>



#ifdef RAD_VIS
#include <strstream.h>
#endif

#ifdef DEBUG
#define DBG_COUT if (1) cerr
#define DBG_COUT2 if (0) cerr
#else
#define DBG_COUT if (0) cerr
#define DBG_COUT2 if (0) cerr
#endif

HierRad::HierRad() : RadMethod(), cluster(0), leaves(0)
{
}

HierRad::~HierRad()
{
    delete cluster;
    delete leaves;
}

Void HierRad::SetScene(scScenePtr theScene)
{
    Int     i;

    DBG_COUT << "setting HR scene..." << endl;


    clusterableElems.Clear();
    
    RadMethod::SetScene(theScene);      // this calls ResetOptions too

    // we don't do any remeshing, so collect stats here.
        
    stats.Init();
    for (i = 0; i < baseElems.NumItems(); i++)
        stats.Update(baseElems[i]);
#ifdef DEBUG
    stats.Report(cerr);
#endif
}

Void HierRad::ResetOptions()
{
    Int     i, j;
    
    gRadControl->finalPass = false;
    
    if (clusterableElems.NumItems() > 0 
        )
    {
        DBG_COUT << "resetting existing elements" << endl;

        // delete any previous meshing/clusters/etc.
        for (i = clusterableElems.NumItems() - 1; i >= 0; i--)
            clusterableElems[i]->Reset();   
        if (cluster)
            cluster->Reset();
    }
    else
    {
        DBG_COUT << "creating elements" << endl;

        delete cluster;
        cluster = 0;

        // XXX cluster + higher order not supported yet.
        if (gRadControl->basis > kHaar)
            gRadControl->cluster = false;

        // add patches to the elements list
        DBG_COUT << "base elems = " << baseElems.NumItems() << endl;
        clusterableElems.SetSize(baseElems.NumItems());
        // upcast from RadElem pointer to HRElem pointer; we can assume that
        // all patches are a subclass of HRMeshElem.
        for (i = 0; i < baseElems.NumItems(); i++)
            clusterableElems[i] = (HRElemPtr) (HRMeshElem*) baseElems[i];
        

        DBG_COUT << "total HR elems = " << clusterableElems.NumItems() << endl;
    }

    numTrees = 0;
    leaves = 0;

    ColourMeshInitial();
    ColourVertices();

#ifdef RAD_VIS
    gRadControl->pvData = 0;
    Field(out1) << clusterableElems.NumItems() << " patch elements";
    Field(out1) << show;
#endif
}


Void HierRad::Draw(Renderer &r)
{
    Bool    saveWire = gRadControl->wire;

    if (gRadControl->patchView && gRadControl->pvData)
    {
        // draw contributors to this element
        gRadControl->outlineClusters = true;
        gRadControl->pvData->DrawContributors(r);               
        gRadControl->outlineClusters = false;
        gRadControl->wire = true;
    }
            
    // draw traditional patches
    RadMethod::Draw(r);


    // draw cluster last if it exists, so as to get transparency effects
    if (gRadControl->cluster && cluster)
        cluster->DrawElem(r);

    gRadControl->wire = saveWire;
}

Void HierRad::DrawMatrix(Renderer &r)
{
#ifdef RAD_VIS
    if (numTrees == 0)
        return;

        
    if (!gRadControl->drawMatrix || gRadControl->stop
            || gRadControl->cluster)
        return;

    Int i;
    
    for (i = 0; i < numTrees; i++)
        trees[i]->DrawMatrix(r, numTrees);
#endif
}

HRLink *HierRad::CreateLink()
{
    HRLink *result;
    
    //  Create a top-level link between two patches
    
    if (gRadControl->cluster)
        result = new GenLink;
    else if (gRadControl->basis == kHaar)
        result = new HaarLink;
    else
        _Error("(HierRad::CreateLink) Unsupported basis");
        
    return(result);
}


// --- Rendering methods ------------------------------------------------------

// NOTE: All these routines return false if the render should
// terminate.

Bool HierRad::CreateClusterLink()
// Cluster the scene, and create the single initial top-level
// link.
{
    HRLink      *newLink;
    
    // Create cluster
    if (Stage(2)) return(false);
    if (!cluster)
    {
        cluster = new Cluster;
        cluster->SetElems(&clusterableElems);
        cluster->CreateHierarchy();
    }
            
    clusTree = cluster;
    trees = &clusTree;
    numTrees = 1;

    // Create single top-level link 
    linkPatch = 0;
    if (Stage(3)) return(false);

    newLink = new GenLink;

    newLink->MakeLink(cluster, cluster);
    cluster->links.Prepend(newLink);

    return(true);
}

Bool HierRad::CreateInitialPatchLinks()
// Creates top-level links between initial base polygons
{
    Int         i, j;
    HRLink      *newLink1, *newLink2;

    // We only want to test visibility once
    // between each pair of patches.

    numTrees = clusterableElems.NumItems();
    trees = (HRElemPtr*) clusterableElems.Ref();
    if (Stage(2)) return(false);

    for (i = 0; i < numTrees; i++)
    {
        linkPatch = i;
        if (Stage(3)) return(false);

        for (j = i + 1; j < numTrees; j++)
        {

            newLink1 = 0;
            if (trees[i]->EltRho() != cBlack)
            {
                // Add link from j -> i
                            
                newLink1 = CreateLink();
                if (newLink1->MakeLink(trees[j], trees[i]))
                {
                    DBG_COUT2 << "link " << trees[j]->id << " to " 
                        << trees[i]->id << endl;
                    trees[i]->links.Prepend(newLink1);
                }
                else
                {
                    newLink1->FreeLink(newLink1);
                    DBG_COUT << "top-level link REDACTED" << endl;
                    continue;   // don't bother with the other direction.
                }
            }       
            
            if (trees[j]->EltRho() != cBlack)
            {
                // Add link from i -> j, reusing i->j link visibility
                // if possible.

                newLink2 = CreateLink();
                if (newLink2->MakeLink(trees[i], trees[j], newLink1, true))
                {
                    DBG_COUT2 << "link " << trees[i]->id 
                        << " to " << trees[j]->id << endl;
                    trees[j]->links.Prepend(newLink2);
                }
                else
                {
                    DBG_COUT << "top-level link REDACTED" << endl;
                    newLink2->FreeLink(newLink2);
                }
            }       
        }           
    }
    
    return(true);
}

Bool HierRad::TwoStageSolve()
{
    Int         i;
    GCLReal     origErr;
    GCLReal     resErr;
    
    origErr = 1.0;
    error = 1.0;
    gRadControl->refineLevels = 32;
    if (gRadControl->solver != sv_scheduled)
        iterations = 0;
    if (Stage(5)) return(false);    

    //  Refine the links in the scene
    if (Stage(6)) return(false);
    for (i = 0; i < numTrees; i++)
    {
        linkPatch = i;
        if (Stage(7)) return(false);
        trees[i]->Refine();
    }

#ifdef RAD_RESTART_FROM_EMITTERS
    if (gRadControl->solver == sv_scheduled && iterations > 0)
    {
        DBG_COUT << "resetting radiosities" << endl;

        for (i = 0; i < numTrees; i++)
            trees[i]->InitRad();
    }
#endif
    
    // Solve for radiosities... 
    do
    {
        resErr = 0.0;
        if (Stage(5)) return(false);

        //  Do a gather over the whole scene
        if (Stage(8)) return(false);
        for (i = 0; i < numTrees; i++)
        {
            poly = i;
            if (Stage(9)) return(false);
#ifdef RAD_VIS
//          if (doUpdate && gRadControl->showLinks)
//              display->Clear().Draw(gRadControl->scene);
#endif
            trees[i]->GatherAll();
#ifdef RAD_VIS
//          if (doUpdate && gRadControl->showLinks)
//              display->Show();
#endif
        }

        //  Now reconstruct the basis functions with a push/pull.
        //  Separating the gather and push/pull stages like
        //  this is the equivalent of Jacobi iteration, which is
        //  more memory and time consuming than Gauss Seidel. Our
        //  experiments with Gauss Seidel have shown that it can
        //  lead to some instability problems however, especially
        //  when there are a small number of very bright lights.
            
        if (Stage(10)) return(false);
        for (i = 0; i < numTrees; i++)
        {
            poly = i;
            if (Stage(11)) return(false);
            trees[i]->ClearB();
            trees[i]->PushPull();
            resErr += abs(trees[i]->Error());
        }
        
        if (origErr == 1.0 && resErr > 0.0)
            origErr = resErr;
        error = resErr / origErr;
        if (!gRadControl->solver == sv_scheduled)
            iterations++;
            
        if (Stage(12)) return(false);       
    }
    while (gRadControl->solver != sv_scheduled && error > gRadControl->error);  // UUU

    return(true);
}

Bool HierRad::MultiGridSolve()
{
    Int         i;
    Bool        subdivided;
    GCLReal     origErr;
    GCLReal     resErr;
    
    // Here we mix a single refine stage with a single solve stage, and
    // repeat until all links have been fully refined.
    
    origErr = 1.0;
    error = 1.0;
    subdivided = true;
    iterations = 0;

    // make sure mixed clusters etc. are refined before we start...
    // a bit hacky...
    if (gRadControl->cluster)
    {
        Double save1 = gRadControl->kAError;
        Double save2 = gRadControl->kFError;

        gRadControl->refineLevels = 32;
        gRadControl->kAError = 1e9;
        gRadControl->kFError = 1e9;
        cluster->Refine();
        gRadControl->kAError = save1;
        gRadControl->kFError = save2;
    }

    // refine only one level at a time.
    gRadControl->refineLevels = 1;
    
    do
    {
        resErr = 0.0;
        iterations++;
        if (Stage(5)) return(false);

        //  Do a gather over the whole scene
        if (Stage(8)) return(false);
        for (i = 0; i < numTrees; i++)
        {
            poly = i;
            if (Stage(9)) return(false);
#ifdef RAD_VIS
//          if (doUpdate && gRadControl->showLinks)
//              display->Clear().Draw(gRadControl->scene);
#endif
            trees[i]->GatherAll();
#ifdef RAD_VIS
//          if (doUpdate && gRadControl->showLinks)
//              display->Show();
#endif
        }
            
        // Push/pull
        if (Stage(10)) return(false);
        for (i = 0; i < numTrees; i++)
        {
            poly = i;
            if (Stage(11)) return(false);
            trees[i]->ClearB();
            trees[i]->PushPull();
            resErr += abs(trees[i]->Error());
        }
        
        if (origErr == 1.0)
            origErr = resErr;
        error = resErr / origErr;
        if (Stage(12)) return(false);

        if (!subdivided)    // don't refine if we've subdivided all we can.
            continue;
                        
        //  Refine the links in the scene
        if (Stage(6)) return(false);        
        subdivided = false;
        for (i = 0; i < numTrees; i++)
        {
            linkPatch = i;
            if (Stage(7)) return(false);
            if (trees[i]->Refine())
                subdivided = true;
        }

        if (iterations > 30) break; // XXX
    }
    while ((gRadControl->refAllLinks && subdivided) || error > gRadControl->error);
    
    return(true);
}

Bool HierRad::ScheduledSolve()
{
    GCLReal     origErr;
    
    // Here we repeat a two stage solve, dropping epsilon by a
    // set amount on each occasion.
    
    origErr = gRadControl->kFError;
    iterations = 0;
        
    do
    {
        if (!TwoStageSolve()) return(false);
        iterations++;
        gRadControl->kFError /= gRadControl->alpha;
#ifdef RAD_VIS
        Field(out1) << "epsilon = " << gRadControl->kFError << show;        
#endif      
    } while (iterations < gRadControl->schedIterations && !Pause());

    gRadControl->kFError = origErr;

    return(true);
}

Bool HierRad::Render()
{
    Bool    result;
    Int     i;
    
    //  Initialisation
    
    gRadControl->useBF = !(gRadControl->solver == sv_twoStage);
    linksCreated = 0;
    linksCreatedLast = 0;

    if (gRadControl->basis != kHaar) 
    // turn off anchoring & balancing for non-haar meshes
    {
        gRadControl->anchor = false;
        gRadControl->graded = false;
    }

    RadMethod::Render();

    if (Stage(1)) return(false);

    //  Set up initial emission of the patches
    if (gRadControl->cluster)
        result = CreateClusterLink();
    else
        result = CreateInitialPatchLinks();

    if (!result)
        return(false);

    for (i = 0; i < numTrees; i++)
        trees[i]->InitRad();

    if (Stage(4)) return(false);

    if (gRadControl->solver == sv_scheduled)
        result = ScheduledSolve();
    else if (gRadControl->solver == sv_interleaved)
        result = MultiGridSolve();
    else if (gRadControl->solver == sv_twoStage)
        result = TwoStageSolve();
    else
        _Error("Unknown solver algorithm");
        
    if (Stage(13)) return(false);
        
    return(result);
}


#ifdef RAD_VIS
Void HierRad::DumpStats()
{
}

Int HierRad::Stage(Int stage)
{
    Int         i;
    ostrstream  outMsg;
    Char        tempStr[256];
    Bool        animate = gRadControl->animate || gRadControl->pause;
    
    switch (stage)
    {
    case 1: // pre setup
        Field(out2) << "Starting patches: " << numTrees << show;
        Field(out1) << "Initialising patches." << show;
        ColourVertices();
        StartUpdate();
        break;

    case 2: // create links (post setup)
        if (gRadControl->cluster)
            Field(out1) << "Creating clusters..." << show;
        else
            Field(out1) << "Creating links... " << show;
        ColourVertices();
        display->Redraw();
        break;

    case 3: // create a link 
        doUpdate = animate || Update();
        if (doUpdate)
        {
            RenderMatrix();
            if (animate)
                Field(out1) << "Creating links for element " << linkPatch 
                    << "/" << numTrees << show;
            else
                Field(out1) << "Created " << linkPatch << " links out of "
                            << numTrees << show;
        }
        break;

    case 4: // solve system with refinement
        Field(out1)  << "Solving system with refinement" << show;
        ColourVertices();
        display->Redraw();
        break;
        
    case 5: // solve system iteration
        RenderMatrix();
        Field(out1)  << "Starting iteration " << iterations << show;
        Field(out2)  << "Iteration " << iterations << show;
        doUpdate = animate || Update();
        break;
        
    case 6: // refine links
        Field(out1)  << "Refining links" << show;
        break;
        
    case 7: // refine 1 link
        doUpdate = animate || Update();
        if (doUpdate)
        {
            Field(out1)  << "Refining links for element " << linkPatch << show;
            Field(out3)  << "Total links created so far: " << linksCreated << show;
            RenderMatrix();
        }
        break;

    case 8: // gather
        Field(out1)  << "Gathering radiosity" << show;
        ColourVertices();
        display->Redraw();
        break;
                
    case 9: // gather across 1 link
        doUpdate = animate || Update();
        if (doUpdate)
        {
            Field(out1)  << "Gathering to element " << poly << show;
            ColourVertices();
            display->Redraw();
        }
        break;
                    
    case 10: // push/pull
        Field(out1) << "Push/Pulling hierarchies" << show;
        break;
        
    case 11: // push/pull 1
        doUpdate = animate || Update();
        if (doUpdate)
        {
            Field(out1) << "Push/Pull element " << poly << show;
            ColourVertices();
            display->Redraw();
        }
        break;  
        
    case 12: // finish one system iteration
        {
            HRStats     stats;
            
            ColourVertices();
            display->Redraw();

            stats.Init();
            for (i = 0; i < numTrees; i++)
                 trees[i]->CalcStats(&stats);

            Field(out2) << "Overall error: " << error << show;

            outMsg << "Iterations: " << iterations << 
                ", patches: " << stats.numPatches;
            if (stats.numVolClusters > 0)
                outMsg << ", vc: " << stats.numVolClusters;
            if (stats.numFaceClusters > 0)
                outMsg << ", fc: " << stats.numFaceClusters;
            outMsg << ", links: " << stats.numLinks;
            memcpy(tempStr, outMsg.str(), outMsg.pcount());
            tempStr[outMsg.pcount()] = 0;
            
            Field(out3) << tempStr << show;
        }
        break;
    case 13: // all done
        Field(out1) << "Colouring models...." << show;
        gRadControl->finalPass = true;
        ColourVertices();
        Field(out1) << "Finished." << show;
        display->Redraw();
        break;
        
    }
    
    return(Pause());
}

#else

static Char *tBasis[] = {"none", "haar", "f2", "f3", "m2", "m3"};
static Int tBasisOrder[] = {0, 1, 2, 3, 2, 3};

static GCLReal gAlpha = 0.0;
static Int gLinks0 = 0;

Void HierRad::DumpStats()
{
    Int         i, m = tBasisOrder[gRadControl->basis];
    GCLReal     mem, linkMem, patchMem, volClusMem;
    HRStats     stats;
    
    stats.Init();
    for (i = 0; i < numTrees; i++)
        trees[i]->CalcStats(&stats);
    
    // Calculate projected memory use
#ifdef RAD_THEORY_MEM
    linkMem = sizeof(Void*) * 2 + sizeof(GCLReal) + sizeof(Colour) * 
        sqr(sqr(m));
    patchMem = sizeof(Void*) * 5 + sizeof(Colour) * sqr(m) * 2;
    // 4 children, one link list pointer, m colour coefficients * 2.
    faceClusMem = sizeof(Void*) * 2 + sizeof(Colour) * 6;
    // 2 children, 3 colour coeffs * 2 
    volClusMem = sizeof(Void*) * 8 +  sizeof(Colour) * 2;
    // 8 children, 1 HRElem pointer, 1 colour coeffs * 2
#else
    linkMem = sizeof(GenLink);
    patchMem = sizeof(HRMeshElem);
    volClusMem = sizeof(Cluster);
#endif

    mem = linkMem * stats.numLinks
        + patchMem * stats.numPatches
        + volClusMem * stats.numVolClusters
        ;   


    mem /= 1024.0;
    
    if (gAlpha <= 0.0)
        gLinks0 = stats.numLinks;

    if (gAlpha < 0.0)
        gAlpha = stats.numLinks / sqr(gRadControl->numPolys);
        
    cout << dumpID
         << ' ' << totTime
         << ' ' << gRadControl->stage
         << ' ' << iterations
         << ' ' << stats.numPatches
         << ' ' << stats.numVolClusters
         << ' ' << stats.numLinks
         << ' ' << linksCreated
         << ' ' << gAlpha
         << ' ' << stats.numLinks - gLinks0

         << ' ' << error
         << ' ' << gRadControl->rays
         << ' ' << mem
         << ' ' << grid->MemoryUsage()
         << ' ' << TotalMemoryUse()
         << endl;

    DumpScene();
}

static Char *gSolverNames[] =
{
    "twoStage",
    "interleaved",
    "scheduled",
    0
};

Int HierRad::Stage(Int stage)
{
    if (CheckTime()) return(1);
    
    gRadControl->stage = stage;

    switch (stage)
    {
    case 2:
        cout << "method hier " << tBasis[gRadControl->basis] << endl;
        if (gRadControl->cluster) cout << "clustering" << endl;
        cout << "solver " << gSolverNames[gRadControl->solver] << endl;
        
        cout << "tErr " << gRadControl->error << endl;
        cout << "fErr " << gRadControl->kFError << endl;
        cout << "aErr " << gRadControl->kAError << endl;
        cout << "dErr " << gRadControl->dFError << endl;
        cout << "vErr " << gRadControl->visError << endl;
        cout << "scene " << sceneName << endl;
        cout << "polys " << gRadControl->numPolys << endl;
        cout << "srcPatches " << numTrees << endl;
        // must match with list in DumpStats
        cout << "format "
             << "ID "
             << "time "
             << "stage "
             << "iterations "
             << "patches "
             << "volClusters "
             << "links "
             << "totLinks "
             << "sparsity "
             << "newLinks "
             << "resErr "
             << "rays "
             << "mem "
             << "rtMem "
             << "rss "
             << endl;
        cout << "----------------------------------------"
             << "---------------------------------------" << endl;

        iterations = 0;
        error = 0;
        gRadControl->finalPass = false;
        DumpStats();
        break;
                
    case 4: 
        // At this point the first-pass linking has been done.
        gAlpha = -1.0;  // signal the DumpStats procedure
        DumpStats();
        break;
        
    case 13:
        if (gRadControl->dumpTree)
        {
            ofstream s;
            cerr << "Saving final hierarchy as " << (gRadControl->outFile) + ".hier" << endl;
            s.open(StrConst(gRadControl->outFile) + ".hier");
            HRElem::out = &s;
            s << "# hierarchy dump for " << sceneName << endl;
            s << "# e <ID> <childID1> <childID2>...   -> define element" << endl;
            s << "# l <fromID> <toID>                 -> define link" << endl;
            for (Int i = 0; i < numTrees; i++)
                trees[i]->DumpHierarchy(0);
            s.close();
            HRElem::out = &cout;
        }
        gRadControl->finalPass = true;
        DumpStats();
        timer.ContTimer();
        GetElements();  // final pass gets performed here.
        timer.StopTimer();
        totTime = timer.GetTimer();
        gRadControl->stage++;
        DumpStats();
        
        break;
    }
    
    if (Idle()) return(0);
    timer.ContTimer();
    return(0);
}

#endif

RadElem *HierRad::NewMesh()
{
    RadElem *elem;

    if (gRadControl->basis == kHaar)
        elem = new HaarElem();
    else
        _Error("(HierRad::NewMesh) Unhandled basis");

    return(elem);
}

Void HierRad::ColourVertices()
{
    Int             i, j, oldRays;
    ProgramTimer    timer;
    // colour underlying mesh & models

    oldRays = gRadControl->rays;
    if (gRadControl->finalPass)
    {
        if (gRadControl->shotDisplay)
            gRadControl->showRays = true;
        timer.StartTimer();

    }

    // colour leaves
    if (gRadControl->finalPass && gRadControl->bestPass)
        for (i = 0; i < numTrees; i++)
        {
            ShadeInfo shadeInfo;

            trees[i]->DistributeColoursBest(shadeInfo);
        }   
    else
        for (i = 0; i < numTrees; i++)
            trees[i]->DistributeColours();

    for (i = 0; i < baseElems.NumItems(); i++)
        ((HRMeshElem*) baseElems[i])->HRCorrectLeaves();


    if (gRadControl->finalPass)
    {
        DBG_COUT << "Done doing final distribute : " << timer.GetTimer() << "s" << endl;
        DBG_COUT << "Cast " << gRadControl->rays - oldRays << " more rays." << endl;
    }
    
    // for meshes, as before...

    RadMethod::ColourVertices();
}


PatchList *HierRad::GetElements()
{
    if (!leaves)
    {
        Int             i, j, n, numSubElems, numSplitElems = 0;
        HierElemPtr     subElems[5];
        
        ColourVertices();
        

        leaves = new PatchList;
        for (i = 0; i < baseElems.NumItems(); i++)
            ((HierElemPtr) baseElems[i])->CollectLeaves(*leaves);

        if (gRadControl->anchor)
        {
            n = leaves->NumItems();

            for (i = 0; i < n; i++)
            {
                ((HierElemPtr) leaves->Item(i))->SplitElement(numSubElems, subElems);

                if (numSubElems > 0)
                {
                    numSplitElems++;
                    
                    leaves->Item(i) = subElems[0];
                    for (j = 1; j < numSubElems; j++)
                        leaves->Append(subElems[j]);
                }
            }
        }
    }

    return(leaves);
}

Void HierRad::PrintCmds(ostream &s) const
{
    Int     i;

    RadMethod::PrintCmds(s);

}

Bool HierRad::ParseCmd(StrConst str, istream &s)
{
    // XXX broken at the mo'
    return(RadMethod::ParseCmd(str, s));
}


HRElem *HierRad::ClosestHRIntersection(Point &start, Point &dir, Point &p)
{
#ifdef RAD_VIS
    if (grid)
    {
        if (grid->ClosestIntersection(start, dir))
        {
            p = grid->hitT * dir + start;

            if (grid->hitPrim->object->id == 0)
            {
                if (rtElemPtrs.NumItems() > 0)
                    return(dynamic_cast<HRMeshElem*>(rtElemPtrs[grid->hitPrim->id]));
            }
        }
    }
#endif

    return(0);
}

---