radgen.cpp

使用stl技术,(还没看,是听说的)

// ---------------------------------------------------------------------------------------------------------------------------------
//  _____            _  _____                                 
// |  __ \          | |/ ____|                                
// | |__) | __ _  __| | |  __  ___ _ __       ___ _ __  _ __  
// |  _  / / _` |/ _` | | |_ |/ _ \ '_ \     / __| '_ \| '_ \ 
// | | \ \| (_| | (_| | |__| |  __/ | | | _ | (__| |_) | |_) |
// |_|  \_\\__,_|\__,_|\_____|\___|_| |_|(_) \___| .__/| .__/ 
//                                               | |   | |    
//                                               |_|   |_|    
//
// Description:
//
//   Radiosity generation main code
//
// Notes:
//
//   Best viewed with 8-character tabs and (at least) 132 columns
//
// History:
//
//   08/17/2001 by Paul Nettle: Original creation
//
// Restrictions & freedoms pertaining to usage and redistribution of this software:
//
//   This software is 100% free. If you use this software (in part or in whole) you must credit the author. This software may not be
//   re-distributed (in part or in whole) in a modified form without clear documentation on how to obtain a copy of the original
//   work. You may not use this software to directly or indirectly cause harm to others. This software is provided as-is and without
//   warrantee -- Use at your own risk. For more information, visit HTTP://www.FluidStudios.com/
//
// Copyright 2002, Fluid Studios, Inc., all rights reserved.
// ---------------------------------------------------------------------------------------------------------------------------------

#include "stdafx.h"
#include "FSRad.h"
#include "ProgressDlg.h"
#include "GeomDB.h"
#include "BSP.h"
#include "SOctree.h"
#include "RadGen.h"
#include "BeamTree.h"
#include "PolarQuadtree.h"
#include "LMapGen.h"

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

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

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

static	bool	_processEnergy(SOctree & node, void * userData)
{
	RadGen *	rg = static_cast<RadGen *>(userData);
	if (!rg) return false;

	return rg->processEnergy(node);
}

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

void	RadGen::calcVisiblePieces(const RadPatch & emitter, const RadPrimPointerListGrainy & potentialReceivers, RadPrimListGrainy & visiblePieces)
{
#if 0
	// Insert the polygons into the polar QT for visibility

	geom::Plane3		axis = emitter.plane();
	axis.origin() = emitter.origin();
	PolarQuadtree		pqt(axis);

	// The potential receivers is an ordered list from back-to-front, so go through it backwards...

	for (RadPrimPointerListGrainy::node * i = potentialReceivers.tail(); i; i = i->prev())
	{
		pqt.insert(*i->data());
	}

	visiblePieces += pqt.visiblePieces();
#else
	// Insert the polygons into the beamtree for visibility

	BeamTree		bt;
	geom::Plane3Array	pa;
	pa += emitter.plane();
	bt.init(pa);
	geom::Plane3Array	planes;

	// The potential receivers is an ordered list from back-to-front, so go through it backwards...

	for (RadPrimPointerListGrainy::node * i = potentialReceivers.tail(); i; i = i->prev())
	{
		bt.insert(emitter.origin(), *i->data(), visiblePieces);
	}
#endif
}

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

void	RadGen::emit(RadPatch & emitter, const RadPrimPointerListGrainy & potentialReceivers, RadLMapArray & storedEnergy)
{
	// What's visible from this emitter?

	RadPrimListGrainy	visiblePieces;
	calcVisiblePieces(emitter, potentialReceivers, visiblePieces);

	// Clear out the emitter's energy before we shoot. Theoretically, we never should need to do this, because an emitter
	// should never be able to emit energy onto itself, but it won't hurt to do this for completeness.

	geom::Color3	shootingEnergy = emitter.energy();
	geom::Color3	reflectedEnergy(0,0,0);
	geom::Color3	receivedEnergy(0,0,0);
	emitter.energy() = geom::Color3(0,0,0);

	// Light up the visible pieces
	{
		for (RadPrimListGrainy::node * i = visiblePieces.head(); i; i = i->next())
		{
			RadPrim &	prim = i->data();
			RadPrim &	originalPrim = *i->data().originalPrimitive();
			RadLMap &	se = storedEnergy[originalPrim.textureID()];

			prim.prepareNoPatches();
			unsigned int	fIndex = 0;

			float	totalArea = (prim.uXFormVector() % prim.vXFormVector()).length();

			geom::Point3	elementCenter = prim.minXYZ() + prim.uXFormVector() * 0.5 + prim.vXFormVector() * 0.5;

			unsigned int	minU = static_cast<unsigned int>(prim.minUV().u());
			unsigned int	minV = static_cast<unsigned int>(prim.minUV().v());
			unsigned int	maxU = static_cast<unsigned int>(prim.maxUV().u());
			unsigned int	maxV = static_cast<unsigned int>(prim.maxUV().v());

			unsigned int	ominU = static_cast<unsigned int>(originalPrim.minUV().u());
			unsigned int	omaxU = static_cast<unsigned int>(originalPrim.maxUV().u());
			unsigned int	ominV = static_cast<unsigned int>(originalPrim.minUV().v());

			for (unsigned int v = minV; v <= maxV; ++v)
			{
				geom::Point3	ec = elementCenter;
				for (unsigned int u = minU; u <= maxU; ++u, ++fIndex, ec += prim.uXFormVector())
				{
					// Receiver must have area

					float	area_j = prim.elementAreas()[fIndex];
					if (area_j <= 0) continue;

					// Ignore same planes

					if ((prim.plane().normal() ^ emitter.plane().normal()) > 0.009f) continue;

					// Get the whole area of this patch from the original polygon

					unsigned int	pu = u - ominU;
					unsigned int	pv = v - ominV;
					unsigned int	elementIndex = pv * (omaxU-ominU+1) + pu;

					// Here we calculate the radiosity equation for a differential area form factor.
					// It looks like this:
					//
					// Fij = (cos(Theta_i) * cos(Theta_j)) / (PI * distanceSquared) * Hij * dAj
					//
					// Where:
					//	Fij is the form factor from the transmission patch (i) to the receiving patch (j)
					//	theta_i is the angle between the normal of (i) and the direction from (i) to (j)
					//	theta_j is the angle between the normal of (j) and the direction from (j) to (i)
					//	Hij is the parametric scalar representing the amount of (j) that is visible to (i)
					//	dAj is the differential area (i.e. the area of (i) divided by area of (j))

					// Vector from i to j

					geom::Vector3	ijVect = ec - emitter.origin();
					float		distanceSquared = ijVect.lengthSquared();
					ijVect.normalize();

					float	cosTheta_i = emitter.plane().normal() ^ (ijVect);
					float	cosTheta_j = prim.plane().normal() ^ (-ijVect);
					float	Hij = area_j / totalArea;
					float	dAj = totalArea;

					// Our differential form factor

					float	Fij = (cosTheta_i * cosTheta_j) / (fstl::pi<float>() * distanceSquared) * Hij * dAj;

					// The energy received

					geom::Color3	energy = shootingEnergy * Fij;
					geom::Color3	reflected = energy * originalPrim.reflectanceColor();

					// Keep track of the energy going back into the scene

					receivedEnergy += energy;
					reflectedEnergy += reflected;

					// Store the energy in the element

					se.data()[v*se.width()+u] += reflected * static_cast<float>(areaLightMultiplier());

					// Transfer the energy up to the patch

					if (!directLightOnly())
					{
						unsigned int	pu = u - ominU;
						unsigned int	pv = v - ominV;
						unsigned int	patchIndex = pv/subdivisionV() * originalPrim.uPatches() + pu/subdivisionU();
						RadPatch & patch = originalPrim.patches()[patchIndex];
						patch.energy() += reflected;
					}
				}

				elementCenter += prim.vXFormVector();
			}
		}
	}

	if (directLightOnly()) reflectedEnergy = geom::Color3(0,0,0);

	energyEscapedThisIteration() = shootingEnergy - receivedEnergy;
	energyAbsorbedThisIteration() = shootingEnergy - energyEscapedThisIteration() - reflectedEnergy;
}

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

void	RadGen::emitNusselt(RadPatch & emitter, const RadPrimPointerListGrainy & potentialReceivers, RadLMapArray & storedEnergy)
{
	// What's visible from this emitter?

	RadPrimListGrainy	visiblePieces;
	calcVisiblePieces(emitter, potentialReceivers, visiblePieces);

	// Clear out the emitter's energy before we shoot. Theoretically, we never should need to do this, because an emitter
	// should never be able to emit energy onto itself, but it won't hurt to do this for completeness.

	geom::Color3	shootingEnergy = emitter.energy();
	geom::Color3	reflectedEnergy(0,0,0);
	geom::Color3	receivedEnergy(0,0,0);
	emitter.energy() = geom::Color3(0,0,0);

	// Light up the visible pieces
	{
		for (RadPrimListGrainy::node * i = visiblePieces.head(); i; i = i->next())
		{
			RadPrim &	prim = i->data();
			RadPrim &	originalPrim = *i->data().originalPrimitive();
			RadLMap &	se = storedEnergy[originalPrim.textureID()];

			prim.prepareNoPatches();
			unsigned int	fIndex = 0;

			float	totalArea = (prim.uXFormVector() % prim.vXFormVector()).length();

			// We are doing our form factors a bit differently, which allows us to precalculate this :)

			geom::Point3	elementCenter = prim.minXYZ() + prim.uXFormVector() * 0.5 + prim.vXFormVector() * 0.5;

			unsigned int	minU = static_cast<unsigned int>(prim.minUV().u());
			unsigned int	minV = static_cast<unsigned int>(prim.minUV().v());
			unsigned int	maxU = static_cast<unsigned int>(prim.maxUV().u());
			unsigned int	maxV = static_cast<unsigned int>(prim.maxUV().v());

			unsigned int	ominU = static_cast<unsigned int>(originalPrim.minUV().u());
			unsigned int	omaxU = static_cast<unsigned int>(originalPrim.maxUV().u());
			unsigned int	ominV = static_cast<unsigned int>(originalPrim.minUV().v());

			for (unsigned int v = static_cast<unsigned int>(prim.minUV().v()); v <= static_cast<unsigned int>(prim.maxUV().v()); ++v)
			{
				geom::Point3	ec = elementCenter;
				for (unsigned int u = static_cast<unsigned int>(prim.minUV().u()); u <= static_cast<unsigned int>(prim.maxUV().u()); ++u, ++fIndex, ec += prim.uXFormVector())
				{
					// Receiver must have area

					float	area_j = prim.elementAreas()[fIndex];
					if (area_j <= 0) continue;

					// Ignore elements behind the emitter

					if (((ec - emitter.origin()) ^ emitter.plane().normal()) < 0) continue;

					// Ignore same planes

					if ((prim.plane().normal() ^ emitter.plane().normal()) > 0.009f) continue;

					// Get the whole area of this patch from the original polygon

					unsigned int	pu = u - ominU;
					unsigned int	pv = v - ominV;
					unsigned int	elementIndex = pv * (omaxU-ominU+1) + pu;

					// Here we calculate Nusselt's analogy. It is defined as:
					//
					// The form factor of an element (i.e. amount of light the element receives from this patch)
					// is equivalent to the fraction of the unit circle that is formed by taking the projection of
					// the element onto the hemisphere surface and projecting it down onto the circle.
					//
					// Without a diagram, this is hard to understand, so here it is, a bit clearer:
					//
					// We want to know how much energy (from this patch) should be attributed to an element.
					// So you take a unit sphere and center it on the patch. Note that the patch is a planar
					// entity, so we're only concerned with the hemisphere that is on the front-side of the
					// patch. Note also, that since this is a unit hemisphere, the "base" of the hemisphere
					// (where it touches the patch) forms a unit circle.
					//
					// So, Nussalt's analogy is to project a patch onto the surface of the hemisphere. However,
					// doing this is only part of the Equation. This would give elements at sharp angles an
					// equal amount of light as those near the top of the hemisphere (i.e. directly in front of
					// the patch.) This is incorrect; what we need is a non-linear ratio of how much energy each
					// patch gets. So we perform another step -- these elements are then projected orthogonally
					// downward, from the surface of the hemisphere onto its base.
					//
					// Because of this orthogonal projection, we end up with a non-linear scaling of the size
					// of the elements, where those elements near the top of the unit hemisphere will have a larger
					// orthogonal projection than those that were along the sides of the unit hemisphere.
					//
					// We now calculate the area of those final projected elements. If you add up all the areas,
					// you will find that they total the area of a single unit circle, because they were all
					// projected onto a unit circle. So we divide these areas by PI (i.e. the area of a unit circle)
					// and we have the ratio for our form factor.
					//
					// If this is confusing, just peek at a diagram of Nusselt's Analogy, and you'll see it's actually
					// quite simple.

					float	nusseltEnergy = 0;
					{
						geom::Point3	newec = ec - prim.uXFormVector() * 0.5 - prim.vXFormVector() * 0.5;

						geom::Vector3	v0 = newec - emitter.origin();
						v0.normalize();
						geom::Vector3	v1 = newec + prim.uXFormVector() - emitter.origin();
						v1.normalize();
						geom::Vector3	v2 = newec + prim.uXFormVector() + prim.vXFormVector() - emitter.origin();
						v2.normalize();
						geom::Vector3	v3 = newec + prim.vXFormVector() - emitter.origin();
						v3.normalize();

						v0 = emitter.plane().closest(v0 + emitter.origin());
						v1 = emitter.plane().closest(v1 + emitter.origin());
						v2 = emitter.plane().closest(v2 + emitter.origin());
						v3 = emitter.plane().closest(v3 + emitter.origin());

						RadPrim orthoProjection;
						orthoProjection.xyz() += v0;
						orthoProjection.xyz() += v1;
						orthoProjection.xyz() += v2;
						orthoProjection.xyz() += v3;
						nusseltEnergy = orthoProjection.calcArea() / fstl::pi<float>();

						// We now account for the fraction of this patch as it clips to the primitive

						nusseltEnergy *= area_j / totalArea;
					}

					// The energy received

					geom::Color3	energy = shootingEnergy * nusseltEnergy;
					geom::Color3	reflected = energy * originalPrim.reflectanceColor();

					// Keep track of the energy going back into the scene

					receivedEnergy += energy;
					reflectedEnergy += reflected;

					// Store the energy in the element

					se.data()[v*se.width()+u] += reflected * static_cast<float>(areaLightMultiplier());

					// Transfer the energy up to the patch

					if (!directLightOnly())
					{
						unsigned int	pu = u - ominU;
						unsigned int	pv = v - ominV;
						unsigned int	patchIndex = pv/subdivisionV() * originalPrim.uPatches() + pu/subdivisionU();
						RadPatch & patch = originalPrim.patches()[patchIndex];
						patch.energy() += reflected;
					}
				}

				elementCenter += prim.vXFormVector();
			}