Radiosity.cc

Go to the documentation of this file.

#include <gra/reflection/Radiosity.h>
#include <esg/explorer/ShadowExplorer.h>
#include <esg/explorer/ObjsExplorer.h>
#include <assert.h>

using namespace gra;


Radiosity::Patch::Patch(PolygonalEnergy * e,
                        Mesh::Plane     * p,
                        const Vector3   & re,
                        const Vector3   & em,
                        const Vector3   & c,
                        float             a)
    : pEnergy(e),
      pMeshPlane(p),
      reflectance(re),
      emittance(em),
      normal(p->a, p->b, p->c),
      centre(c),
      area(a)
{
}



float Radiosity::_edge_area(const Vector3& p,const Vector3& q,const Vector3& n)
{ 
    /*
     * Numerically stable edge-area calculator. (Including degenerate cases.)
     * Andrew J. Willmott, 1992
     *
     * 4 dots, 1 cross, 1 atan per edge.
     */

    const float epsilon = 1e-6;
    float sinFactor, qdotp, projArea;
    Vector3 aux;

    qdotp =  p.dot(q);
    sinFactor = p.lengthSquared() * q.lengthSquared() - Math::sqr(qdotp);

    aux.cross(p,q);
    aux.negate();
    projArea = aux.dot(n);

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

float Radiosity::_est_patch_factor(Patch&         patch,
                                   const Vector3& point,
                                   const Vector3& normal)
{
    /*
     * Elemental factor from the patch "patch" to point
     * essentially cos theta * cos phi * Ap / (pi * r^2)
     */

    float   result, rad;
    Vector3 temp(point);

    temp.sub(patch.centre);

    result = (patch.normal.x * temp.x +
              patch.normal.y * temp.y +
              patch.normal.z * temp.z);
        
    if (result <= 0) result = 0;
    else {
        result *= -normal.dot(temp);
        if (result <= 0) result = 0;
        else {
            rad = Math::sqr(temp.lengthSquared());
            result *= patch.area;
            result /= (rad * PI); 
        }
    }
        
    return result;
}

float Radiosity::_patch_factor(Patch&         patch,
                               const Vector3& point,
                               const Vector3& normal)
{
    /*
     * 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.
     */

    assert(patch.pMeshPlane);

    Vector3 x1, x2, x3;
    float   sum     = 0.0;
    bool    isFirst = true;
    
    Mesh::Edge * pE = patch.pMeshPlane->any_edge; assert(pE);
    do {
        // x2 = first vertex of the actual edge
        if (pE->LeftLoop == patch.pMeshPlane)
            x2.set(pE->V1->x, pE->V1->y, pE->V1->z);
        else
            x2.set(pE->V2->x, pE->V2->y, pE->V2->z);
        x2.sub(point);

        if (isFirst) { // remember the first vertex
            x3.set(x2); 
            isFirst = false;
        } else 
            sum += _edge_area(x1, x2, normal);
            
        x1.set(x2);  // x1 becames the first vertex of next edge
        
        // go to next the edge in loop
        if (pE->LeftLoop == patch.pMeshPlane) pE = pE->LeftOut;
        else                                  pE = pE->RightOut;
    } while (pE != patch.pMeshPlane->any_edge);

    sum += _edge_area(x1, x3, normal);  // end loop

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

#define _MAX(a,b) (((a)<(b))?(b):(a))
float Radiosity::_est_form_factor(Patch& from, Patch& to)
{
    float result, rad, npt, res2;
    Vector3 temp(to.centre);

    temp.sub(from.centre);
    
    result = -(to.normal.x*temp.x + to.normal.y*temp.y + to.normal.z*temp.z);
    if (result <= 0.0) return 0.0;
        
    npt = (from.normal.x*temp.x + from.normal.y*temp.y + from.normal.z*temp.z);
    if (npt <= 0.0) return 0.0;
    
    rad = Math::sqr(temp.lengthSquared());
    result *= from.area / PI;
    res2 = result * _MAX(npt, 0.25 * sqrt(from.area));

    if (_quadLevel == 0 || rad * _dFError >= res2) return result*(npt/rad);
    
    // result is blowing up!
    // let's bail to the patch factor formula

    Vector3 v;
    float   sum  = 0.0;
    int     numV = 0;
    
    if (_quadLevel > 1) {
        // use more-accurate PF estimate:
        // sample corners
        Mesh::Edge * pE = to.pMeshPlane->any_edge; assert(pE);
        do {
            if (pE->LeftLoop == to.pMeshPlane)
                v.set(pE->V1->x, pE->V1->y, pE->V1->z);
            else
                v.set(pE->V2->x, pE->V2->y, pE->V2->z);
            sum += _patch_factor(from, v, to.normal);
            numV++;
            // go to next the edge in loop
            if (pE->LeftLoop == to.pMeshPlane) pE = pE->LeftOut;
            else                               pE = pE->RightOut;
        } while (pE != to.pMeshPlane->any_edge);
        
        if (_quadLevel > 2) {
            // estimate ff from tetrahedron constructed from
            // corner samples and centre. this will
            // underestimate the true ff slightly.
            sum += 2.0 * _patch_factor(from, to.centre, to.normal);
            sum /= numV + 2;
        } else // underestimate curve
            sum /= numV;
    } else
        sum = _patch_factor(from, to.centre, to.normal);
        
    return sum;
}
#undef _MAX

void Radiosity::_make_form_factors(Patch& fromPatch, Vector3 result[])
{
    float rho, factor, visibility;
    int   i = 0;

    for (Patch * curPatch = _scenePatches.firstItem();
         curPatch;
         curPatch = _scenePatches.nextItem())
    {
        rho = curPatch->reflectance.length();
        if (rho < 0.0001) { result[i++].set(0,0,0); continue; }

        factor = _est_form_factor(fromPatch, *curPatch);
        
#ifdef RAD_SPEEDY
        if (rho * factor < 1e-5) { result[i++].set(0,0,0); continue; }
#endif

        visibility = _visibility(fromPatch, *curPatch);
                
        if (visibility > 0.0001) {
            result[i].set(curPatch->reflectance);
            result[i++].scale(factor * visibility);
        } else 
            result[i++].set(0,0,0);
    }
}

float Radiosity::_visibility(Patch& from, Patch& to)
// Returns visibility of the centre of this patch from p.
{
    float   d;

    // Is 'to.centre' completely behind this patch?
    d = from.normal.dot(from.centre);
    if (from.normal.dot(to.centre) <= d) return 0.0;

    // Is this patch completely behind 'to' plane?
    d = to.normal.dot(to.centre);
    if (from.centre.dot(to.normal) <= d) return 0.0;

    // If neither is the case, we have at least partial
    // visibility, so ray cast to estimate it
    assert(_pIntersector);
    
    Vector3 dir(to.centre);
    dir.sub(from.centre);
    float  maxDist = dir.length();
    dir.normalize();
    
    _pIntersector->init(ENV_WANT_INTERFERENCE, maxDist);
#warning "FIX ME: Light OID"
    ShadowExplorer explorer(from.centre, dir, maxDist); 
    if (explorer.result()) return 0.0;
    else                   return 1.0;
}


// ------------- public --------------

void Radiosity::setScene(Scene * pScene)
{
    ReflectionModel::setScene(pScene);
    if (!_pScene || !_pScene->root()) return;

    ObjsExplorer          explorer;
    Matrix4               trMat;
    bool                  transformed;
    Mesh*                 pMesh;
    //EnergyCoat*           pE;
    PolygonalEnergy*      pE;
    AutoPtr<EnergyCoat> * pAE;
    AutoPtr<Mesh>       * pAMesh;
    Vector3               emission;
    Vector3               reflectance;
    MatVisitor            visitor;

    _scenePatches.deleteAll();

    explorer.explore(*(_pScene->root()));
    
    for (SceneGraphObject * pObj = explorer.result(trMat, transformed);
         pObj;
         pObj = explorer.result(trMat, transformed))
    {
        // inspect polygonal energy
        if (!pObj->supportsEnergy()) continue;
        
        pAE = pObj->getEnergyState();
        if (!pAE) continue;

        if (!pAE->origObject() ||
            pAE->origObject()->type() != EnergyCoat::POLYGONAL) continue;
        
        pE = (PolygonalEnergy*) pAE->origObject();
        if (pE->getRegionsDescr() != PolygonalEnergy::FACET_BASED) continue;
        
        pAMesh = pE->getMesh();
        if (!pAMesh) continue;
        
        pMesh = pAMesh->origObject();
        if (!pMesh) continue;

        // inspect diffuse coeficient replacing BRDF
        visitor.init();
        pObj->inspectMaterials(visitor);
        reflectance.set(visitor.diffuse());
        
        // initialize patches for current energy coat
        pMesh->resetActSolid();
        do {
            pMesh->resetActPlane();
            do {
                if (!(pE->getRegionEnergy(pMesh->getActPlaneID(), emission)))
                    //emission.set(0., 0., 0.);
                    fprintf(stderr,"Radiosity::setScene(): Uninitilized facet energy or facet ID mismatch\n");

                //fprintf(stderr,"%i: %f %f %f\n",pMesh->getActPlaneID(),emission.x,emission.y,emission.z);
                _scenePatches.append(new Patch((PolygonalEnergy*) pE,
                                               pMesh->getActPlane(),
                                               reflectance,
                                               emission,
                                               pMesh->getActPlaneCentroid(),
                                               pMesh->getActPlaneArea()));
            } while (pMesh->goToNextPlane());
        } while (pMesh->goToNextSolid());
    }
}

Color3f* Radiosity::illuminatePoint(PointEnv&)
{
    float    power, maxPower;
    Patch *  pMaxPowerPatch = _scenePatches.firstItem();
    Vector3  rgbToLum(1.0/3.0, 1.0/3.0, 1.0/3.0);
    Vector3  radToShoot;
    Vector3  envPower;
    Vector3  deltaRad;
    float    origShoot = 0;
    int      iterations = 0;
    unsigned i;
    Vector3  envRefl(0,0,0);
    float    envArea = 0.0;
    unsigned maxPowerIndex;
    float    error = 1.0;
    Vector3  envRad(1,1,1);
    bool     finished = false;

    if (!_pIntersector) return NULL;

    Vector3* ffRow = new Vector3 [_scenePatches.length()];
    Vector3* S     = new Vector3 [_scenePatches.length()];

    // setup the ambient term + emittance
    i = 0;
    for (Patch * curPatch = _scenePatches.firstItem();
         curPatch;
         curPatch = _scenePatches.nextItem())
    {
        envRefl.add(curPatch->reflectance * curPatch->area);
        envArea += curPatch->area;
        S[i].set(curPatch->emittance);
        i++;
    }
    // R = 1 / (1 - avgReflectance):
    envRefl.x /= 1 / (1 - envRefl.x / envArea);
    envRefl.y /= 1 / (1 - envRefl.y / envArea);
    envRefl.z /= 1 / (1 - envRefl.z / envArea);

    // distribute energy
    while (!finished) {
        maxPower = 0;
        envPower.set(0,0,0);
        maxPowerIndex = 0;
        pMaxPowerPatch = NULL;

        iterations++;

        // find patch with maximal energy
        i = 0;
        for (Patch * curPatch = _scenePatches.firstItem();
             curPatch;
             curPatch = _scenePatches.nextItem())
        {
            power = rgbToLum.dot(S[i] * curPatch->area);
            if (maxPower < power) {
                maxPower = power;
                pMaxPowerPatch = curPatch;
                maxPowerIndex = i;
            }
            i++;
        }

        assert(pMaxPowerPatch);

        // if (idle()) return false;

        radToShoot.set(S[maxPowerIndex]);
        S[maxPowerIndex].set(0,0,0);

        _make_form_factors(*pMaxPowerPatch, ffRow);

        envPower.set(0,0,0);
        error = 0.0;

        // energy shooting
        i = 0;
        for (Patch * curPatch = _scenePatches.firstItem();
             curPatch;
             curPatch = _scenePatches.nextItem())
        {
            deltaRad.x = ffRow[i].x * radToShoot.x;
            deltaRad.y = ffRow[i].y * radToShoot.y;
            deltaRad.z = ffRow[i].z * radToShoot.z;
            curPatch->emittance.add(deltaRad);
            S[i].add(deltaRad);
            error += deltaRad.lengthSquared();
            envPower.add(S[i] * curPatch->area);
            i++;
        }

        if (_useAmbient) {
            envRad.x = envRefl.x * envPower.x / envArea;
            envRad.y = envRefl.y * envPower.y / envArea;
            envRad.z = envRefl.z * envPower.z / envArea;
        } else
            envRad.set(0,0,0);

        error = sqrt(error);
        
        if (iterations == 1) origShoot = maxPower;
        else
            if (maxPower <= 0.01 * origShoot) finished = true;
    } // while


    if (_useAmbient) {
        for (Patch * curPatch = _scenePatches.firstItem();
             curPatch;
             curPatch = _scenePatches.nextItem())
        {
            curPatch->emittance.x += curPatch->reflectance.x*envRad.x;
            curPatch->emittance.y += curPatch->reflectance.x*envRad.y;
            curPatch->emittance.z += curPatch->reflectance.x*envRad.z;
        }
    }
    
    // set energy of shapes
    for (Patch * curPatch = _scenePatches.firstItem();
         curPatch;
         curPatch = _scenePatches.nextItem())
    {
        //fprintf(stderr,"%i, %f %f %f\n",curPatch->pMeshPlane->jmeno_roviny,curPatch->emittance.x,curPatch->emittance.y,curPatch->emittance.z);
        curPatch->pEnergy->setRegionEnergy(curPatch->pMeshPlane->jmeno_roviny,
                                           curPatch->emittance);
    }

    delete [] ffRow;
    delete [] S;

    return NULL;
}

---