---

Cluster.cc

Go to the documentation of this file.

/*
    File:       Cluster.cc
    
    Purpose:    Implements volume clustering for wavelet radiosity
    
    Author:     Andrew Willmott
    
    Notes:      

*/

#include "Cluster.h"
#include "RadMethod.h"
#include "Samples.h"

#include "gcl/VecUtil.h"
#include "gcl/Draw.h"

const Int kMinClusterElems = 32;
const Int kMaxClusterLevels = 16;

// Cluster #defines
#define RAD_NON_ISO_SOURCE
#define RAD_NON_ISO_DEST
// the above should always be defined, except for experiments showing how
// bad the isotropic assumption is =)
#define RAD_SAMPLE_ERROR
// find error by sampling
#define RAD_PTB
// calculate intra-cluster visibility; don't include inter-cluster objects.

#define DBG_COUT if (0) cerr

Cluster::Cluster() : HRElem(), elems()
{
    Int i;
    
    centre = vl_0;
    
    for (i = 0; i < 8; i++)
        child[i] = 0;       

    flags.Set(hrCluster);
}

Cluster::~Cluster()
{
    Int i;

    DBG_COUT << "deleting cluster" << endl;

    for (i = 0; i < 8; i++)
        delete child[i];
}

Void Cluster::Reset()
{
    Int     i;
    
    HRElem::Reset();
    for (i = 0; i < 8; i++)
        if (child[i])
            child[i]->Reset();
}

Void Cluster::CreateBounds()
{
    Int     i;
    
    min = max = elems[0]->EltCentre();
    for (i = 0; i < elems.NumItems(); i++)
        elems[i]->EltUpdateBounds(min, max);

    centre = (min + max) / 2.0;
}

Void Cluster::SetElems(HRElemList *elemsIn)
{
    elems = *elemsIn;
    CreateBounds();

    clusterVol = BoxVol(min, max);
    level = 0;
#ifdef RAD_VIS
    eltParent = 0;
#endif
}

Void Cluster::PushElems()
{
    Vector      d1, d2;
    Int         i, j;
    Point       emin, emax, eCentroid;
    HRElemList  newElems;

    for (j = 0; j < elems.NumItems(); j++)
    {
        emin = emax = elems[j]->EltCentre();
        elems[j]->EltUpdateBounds(emin, emax);
        d1 = emax - emin;
        d2 = (max - min) / 2.0;
        
        DBG_COUT << "size check: " << d1 << " < " << d2 << endl;

        // if the element is bigger in size than an octant, add it to this
        // level of the cluster hierarchy.
        if (d1[0] > d2[0] || d1[1] > d2[1] || d1[2] > d2[2])
        {
            newElems.Append(elems[j]);
            continue;
        }
            
        // find the octant the element's centroid belongs to...
        eCentroid = (emin + emax) / 2.0;
        d1 = eCentroid - centre;
        if (sqrlen(d1) < 1e-8)
        {
            // if it's close to the centre, add it to this level after all.
            newElems.Append(elems[j]);
            continue;
        }
        
        i = 0;
        if (d1[0] >= 0.0) i += 1;
        if (d1[1] >= 0.0) i += 2;
        if (d1[2] >= 0.0) i += 4;
        
        // add it to the given octant
        if (child[i] == 0)
            child[i] = new Cluster;
        child[i]->elems.Append(elems[j]);
    }

    for (i = 0; i < 8; i++)
        if (child[i] && child[i]->elems.NumItems() <= 1)
        {
            DBG_COUT << "REDACTING cluster with only 1 item" << endl;
            newElems.Append(child[i]->elems[0]);
            delete child[i];
            child[i] = 0;
        }
    
    elems.Replace(newElems);
}
 
Cluster *Cluster::CreateHierarchy()
{
    Int         i, n = 0, lightElems;
    Cluster     *lastChild = 0, *curChild;  

    E = vl_0;
    clusterArea = 0.0;
    maxRho = cBlack;
    clusterVol = BoxVol(min, max);
    flags.UnSet(hrMixed);

    if (elems.NumItems() >= kMinClusterElems && level < kMaxClusterLevels)
    {
        PushElems();

        for (i = 0; i < 8; i++) 
            if (child[i]) 
            {   
                
                curChild = child[i];

                child[i]->CreateBounds();
                if (sqrlen(child[i]->max - child[i]->min) < 1e-15)
                    cerr << "WARNING: zero-volume cluster" << endl;
                child[i]->level = level + 1;
#ifdef RAD_VIS
                child[i]->eltParent = this;
#endif
                child[i] = child[i]->CreateHierarchy();
                if (child[i] != curChild)
                    delete curChild;

                if (child[i])
                {
                    E += child[i]->E * child[i]->clusterArea; // sum power
                    clusterArea += child[i]->clusterArea;
                    FindMaxCmpts(child[i]->EltRho(), maxRho, maxRho);
                    flags.Set(child[i]->flags & hrMixed);
                    lastChild = child[i];
                    n++;
                }
            }
    }
    
    // if there are no elements at this level, and we only have one child,
    // we can replace this node with its child.
    if (elems.NumItems() == 0 && n <= 1)
    {
        DBG_COUT << "REDACTING: cluster node with one child cluster" << endl;
        return(lastChild);
    }

    // Add in emitted power of child elements at this level
    
    lightElems = 0;
    for (i = 0; i < elems.NumItems(); i++)
    {
        if (len(elems[i]->EltE()) > 0.0)
        {
            E += elems[i]->EltE() * elems[i]->EltArea();
            // NOTE -- E currently unused. Should only be needed at
            // leaves.
            lightElems++;
        }
        clusterArea += elems[i]->EltArea();
        FindMaxCmpts(elems[i]->EltRho(), maxRho, maxRho);
#ifdef RAD_VIS
        elems[i]->eltParent = this; 
#endif
    }

    if (clusterArea < 1e-15)
    {
        cerr << "REDACTING: cluster with zero surface area" << endl;
        return(0);
    }
    
    if (lightElems > 0)
        flags.Set(hrHasLights);
    if (lightElems < elems.NumItems())
        flags.Set(hrHasNonLights);
    
    // divide by area to give equivalent point src...
    if (clusterArea > 0.0)
        E /= clusterArea;
    else
        E = vl_0;

    DBG_COUT << "--- Cluster at level " << level << " has E = " << 
        E << ", " << elems.NumItems() << " elems, " <<
        n << " sub-clusters." << endl;

    return(this);
}

Void Cluster::Print(ostream &s)
{
    HRLinkIter i;
    
    s << "level " << level << " cluster: " << endl;
    s << elems.NumItems() << " patch children" << endl;
    s << "E = " << E << " B = " << B << endl;
    s << "links: ";
    for (i.Begin(links); !i.AtEnd(); i.Inc())
        s << i.Data();
        
    s << endl;
}

Void Cluster::DrawNodeElem(Renderer &r)
{
    Int     i;

    DrawLeafElem(r);
    for (i = 0; i < 8; i++)
        if (child[i])
            child[i]->DrawElem(r);
}

Void Cluster::DrawLeafElem(Renderer &r)
{
#ifdef RAD_VIS
    if (eltHighlight != 0)
    {
        if (eltHighlight == 1)
            r.C(Colour4(cRed, 0.3));
        else if (eltHighlight == 2)
            r.C(Colour4(cYellow, 0.3));
        else if (eltHighlight == 3)
            r.C(Colour4(cGreen, 0.3));
        PaintBox(r, min, max);

        r.C(HSVCol(hsvYellow + level * 60, 1, 1));
        FrameBox(r, min, max);
    }
    else if (gRadControl->outlineClusters)
    {
        r.C(HSVCol(hsvYellow + level * 60, 1, 1));
        FrameBox(r, min, max);
        
        DrawPatches(r);
    }
#endif
}

Void Cluster::DrawPatches(Renderer &r)
{
    Int i;

    for (i = 0; i < elems.NumItems(); i++)
        elems[i]->DrawLeafElem(r);

    for (i = 0; i < 8; i++)
        if (child[i])
            child[i]->DrawPatches(r);
}

// --- Radiosity transfer -----------------------------------------------------

Colour Cluster::lastB;

Void Cluster::Add()
{
    DBG_COUT << "Cluster add: " << B;
    B += R;
    DBG_COUT << " -> " << B << endl;
}

Void Cluster::Push()
{
    Int i;
    
    for (i = 0; i < 8; i++)
        if (child[i])
            child[i]->B = B;

#ifdef RAD_NON_ISO_DEST
    for (i = 0; i < elems.NumItems(); i++)
        elems[i]->ClearB();
#else
    for (i = 0; i < elems.NumItems(); i++)
        elems[i]->B_Coeffs()[0] = B;
#endif
}

Void Cluster::Pull()
{
    Int     i;
    
    B = vl_0;
    for (i = 0; i < 8; i++)
        if (child[i])
            B += child[i]->EltBA();

    Assert(vl_is_finite(sqrlen(B)), "bad cluster power");

    for (i = 0; i < elems.NumItems(); i++)
    {
        HRElem *fuckYouGDB = elems[i];
        B += elems[i]->EltBA();
        Assert(vl_is_finite(sqrlen(B)), "bad cluster power");
    }
    
    Assert(vl_is_finite(sqrlen(B)), "bad cluster power");

    Expect(vl_is_finite(clusterArea) && (clusterArea > 0), "trivial cluster");
    if (clusterArea > 0.0)
        B /= clusterArea;

    Assert(vl_is_finite(sqrlen(B)), "bad cluster B");
}

Void Cluster::ClearB()
{
    lastB = B;
    B = vl_0;
}

Void Cluster::CalcLeafRadiosity()
{
    _Error("Cluster can't be a leaf");
    // this might change when/if we handle volume-based radiosity
}

Void Cluster::InitRad()
{
    Int     i;
    
    for (i = 0; i < elems.NumItems(); i++)
        elems[i]->InitRad();
        
    for (i = 0; i < 8; i++)
        if (child[i])
            child[i]->InitRad();

    Pull(); 
}

GCLReal Cluster::Error()
{
    return(dot(kRadRGBToLum, B - lastB));
}

Void Cluster::ClearR()
{
    R = vl_0;
}

// --- Link handling ----------------------------------------------------------


Void Cluster::MakeChildLinks(HRElem *other, HRLink *link,
                             Int which, Int levels)
// Splits a link into sublinks that point to the children of
// either 'from' or 'to'.
{
    Int i;

    if (which == kSubTo)
    {
        // add to cluster children
        for (i = 0; i < 8; i++)
            if (child[i])
                child[i]->RefineLink(
                    link->CreateSubLink(other, child[i]),
                levels);
        // add to patch children
        for (i = 0; i < elems.NumItems(); i++)
            elems[i]->RefineLink(
                link->CreateSubLink(other, elems[i]),
            levels);                
    }
    else if (which == kSubFrom)
    {
        // add links to cluster children
        for (i = 0; i < 8; i++)
            if (child[i])
                other->RefineLink(
                    link->CreateSubLink(child[i], other),
                    levels);

        // add to patch children
        for (i = 0; i < elems.NumItems(); i++)
            other->RefineLink(
                link->CreateSubLink(elems[i], other),
                levels);
    }
    else if (which == kSubBoth)
    {
        for (i = 0; i < 8; i++)
            if (child[i])
                other->MakeChildLinks(child[i], link, kSubFrom, levels);
        for (i = 0; i < elems.NumItems(); i++)
                other->MakeChildLinks(elems[i], link, kSubFrom, levels);
    }
    else
        _Error("split type not implemented");
}


// --- children abstraction ---------------------------------------------------


Bool Cluster::IsLeaf()
{
    return(false);
}

Void Cluster::ApplyToChildren(Void (HRElem::*method)(Void*), Void *ref)
{
    Int i;

    for (i = 0; i < 8; i++)
        if (child[i])
            (child[i]->*method)(ref);

    for (i = 0; i < elems.NumItems(); i++)
        (elems[i]->*method)(ref);
}


// --- transport calculations -------------------------------------------------


Void Cluster::EltSampleTransport(
                    Int         numSamples, 
                    Point       p[], 
                    Vector      n[],
                    Matd        &coeffs
                )
{
    Int     i;
    GCLReal result, rLen;
    Vector  r;
#ifdef RAD_NON_ISO_SOURCE
    Vector  dir = vl_0;
#endif

    for (i = 0; i < numSamples; i++)
    {
        r = centre;
        r -= p[i];

        rLen = len(r);
        if (rLen > 0.0)
            r /= rLen;
        result = dot(n[i], r);

        if (result <= 0.0)
        {
            coeffs[i][0] = 0.0;
            continue;
        }

        result /= (sqr(rLen) * vl_pi);
#ifdef RAD_NON_ISO_SOURCE
        // use the sample to weight our average direction normal...
        dir += r;
#endif

        coeffs[i][0] = result;
        Assert(vl_is_finite(coeffs[i][0]), "TS: bad transport coeff!!!");
    }

#ifdef RAD_NON_ISO_SOURCE
    // now that we have an average direction normal, multiply
    // in projected area.
    rLen = len(dir);
    if (rLen > 0.0)
        result = EltProjArea(-dir) / rLen;
    else
        result = 0.0;

    coeffs *= result;
#else
    // if we're doing isotropic, the equivalent projected area
    // is the total projected area over 4. (We're assuming a sphere,
    // basically.)
    coeffs *= clusterArea / 4.0;
#endif

    DBG_COUT << "ClusEST: sampled coeffs = " << coeffs << endl;
}

const Int kNumCubeSamples = 16;

static Int gRotPick = 0;

Void Cluster::EltGetSamples(Int numSamples, Point p[])
{
    Vector      d, *offsets = (Vector*) tCubeSamples;
    Int         i;
    
    d = max - min;

    for (i = 0; i < numSamples; i++)
        p[i] = min + offsets[i] * d;
}

GCLReal Cluster::EltCalcTransport(HRElem *from, Matd &coeffs)
{
    Vector      r;
    GCLReal     rLen, error;

    if (from == this)
    {
        coeffs[0][0] = 1.0;
        error = 1.0;
    }
    else
    {
        r = from->EltCentre();
        r -= centre;
        rLen = len(r);
        if (rLen > 0.0)
            r /= rLen;
            
#ifndef RAD_NON_ISO_DEST
        // Assume isotropic, and thus one quarter of total surface
        // area facing the source.
        r *= 0.25;
#endif

#ifndef RAD_SAMPLE_ERROR
        from->EltSampleTransport(1, &centre, &r, coeffs);
        error = coeffs[0][0];
#else
        Point       p[kNumCubeSamples];
        Vector      n[kNumCubeSamples];
        static Matd trans;
        Int         i, samples;

        EltGetSamples(kNumCubeSamples, p);
        for (i = 0; i < kNumCubeSamples; i++)
            n[i] = r;

        samples = kNumCubeSamples;
        trans.SetSize(samples, from->NumCoeffs());

        // sample transport factors to our set of sample points
        from->EltSampleTransport(samples, p, n, trans);

        error = VecError(col(trans, 0));
        DBG_COUT << "cluster sample error: " << error << endl;

        coeffs[0] = trans[0];
        for (i = 1; i < trans.Rows(); i++)
            coeffs[0] += trans[i];
        coeffs[0] /= trans.Rows();
        Assert(vl_is_finite(coeffs[0][0]), "TC: bad transport coeff!!!");
#endif
    }

    DBG_COUT << "ClusECT: coeffs = " << coeffs << endl;

    return(error);
}

// --- visibility estimation --------------------------------------------------


Void Cluster::EltSetVisPoints(HRElem *to, Point *pts)
{
    Int         i;
    Point       dstPoint[16];
    Vector      d, r;
    Vector      *samples;

    EltGetSamples(16, pts);
    
#ifdef RAD_PTB
    // proj samples to sides of volume facing the source/dest
    if (to && to != this)
    {
        r = to->EltCentre() - centre;
        for (i = 0; i < 16; i++)
            pts[i] += ProjectToBox(min, max, pts[i], r) * r;
    }       
#endif

    if (gRadControl->jitterRot)
        if (++gRotPick == kNumSampArrays)
            gRotPick = 0;
}

Void Cluster::AddIrradiance(const Colour &srcR, const Vector &r)
{
    Int     i;
    
    for (i = 0; i < elems.NumItems(); i++)
        elems[i]->AddIrradiance(srcR, r);

    for (i = 0; i < 8; i++)
        if (child[i])
            child[i]->AddIrradiance(srcR, r);
}

Colour Cluster::GetPower(const Vector &m)
{
    return(B * EltProjArea(m));
}

GCLReal Cluster::EltProjArea(const Vector &n)
{
#ifdef RAD_NON_ISO_DEST
    Int         i;
    GCLReal     eltArea = 0.0;

    for (i = 0; i < elems.NumItems(); i++)
    {
        eltArea += elems[i]->EltProjArea(n);
        Assert(vl_is_finite(eltArea), "FPA: proj area is NAN");
    }

    for (i = 0; i < 8; i++)
    {
        if (child[i])
            eltArea += child[i]->EltProjArea(n);

        Assert(vl_is_finite(eltArea), "FPA: proj area is NAN");
    }

#ifdef DEBUG
    if (elems.NumItems() > 0 && elems[0]->IsFaceClus())
    {
        DBG_COUT << "min eltArea " << n << " = " << eltArea << endl;
    }
#endif

    return(eltArea);
#else
    return(0.25 * clusterArea)
#endif
}

GCLReal Cluster::EltMaxProjArea(const Vector &n)
{
    Int         i;
    GCLReal     eltArea = 0.0;

    for (i = 0; i < elems.NumItems(); i++)
    {
        eltArea += elems[i]->EltMaxProjArea(n);
        Assert(vl_is_finite(eltArea), "FPA: proj area is NAN");
    }

    for (i = 0; i < 8; i++)
    {
        if (child[i])
            eltArea += child[i]->EltMaxProjArea(n);

        Assert(vl_is_finite(eltArea), "FPA: proj area is NAN");
    }

#ifdef DEBUG
    if (elems.NumItems() > 0 && elems[0]->IsFaceClus())
    {
        DBG_COUT << "max eltArea " << n << " = " << eltArea << endl;
        DBG_COUT << "rho = " << EltRho() << endl;
    }
#endif

    return(eltArea);
}

Void Cluster::DistributeColours()
{
    Int     i;

    for (i = 0; i < elems.NumItems(); i++)
        elems[i]->DistributeColours();

    for (i = 0; i < 8; i++)
        if (child[i])
            child[i]->DistributeColours();
}

Void Cluster::DistributeColoursBest(ShadeInfo &shadeInfo)
{
    Int     i;

    if (!links.IsEmpty())
        shadeInfo.linkStack.Push(&links);

    for (i = 0; i < elems.NumItems(); i++)
        elems[i]->DistributeColoursBest(shadeInfo);

    for (i = 0; i < 8; i++)
        if (child[i])
            child[i]->DistributeColoursBest(shadeInfo);

    if (!links.IsEmpty())
        shadeInfo.linkStack.Pop();
}

---