scene.cpp

蒙特卡罗方法可以有效地解决复杂的工程问题

// -----------------------------------------------------------
// scene.cpp
// 2004 - Jacco Bikker - [e-mail] - www.bik5.com -   <><
// -----------------------------------------------------------

#include "common.h"
#include "string.h"
#include "scene.h"
#include "raytracer.h"
#include "stdio.h"
#include "memory.h"
#include "surface.h"

namespace Raytracer {

MManager* KdTree::s_MManager = 0;

// -----------------------------------------------------------
// Texture class implementation
// -----------------------------------------------------------

Texture::Texture( Color* a_Bitmap, int a_Width, int a_Height ) :
	m_Bitmap( a_Bitmap ),
	m_Width( a_Width ), m_Height( a_Height )
{
}

Texture::Texture( char* a_File )
{
	FILE* f = fopen( a_File, "rb" );
	if (f)
	{
		// extract width and height from file header
		unsigned char buffer[20];
		fread( buffer, 1, 20, f );
		m_Width = *(buffer + 12) + 256 * *(buffer + 13);
		m_Height = *(buffer + 14) + 256 * *(buffer + 15);
		fclose( f );
		// read pixel data
		f = fopen( a_File, "rb" );
		unsigned char* t = new unsigned char[m_Width * m_Height * 3 + 1024];
		fread( t, 1, m_Width * m_Height * 3 + 1024, f );
		fclose( f );
		// convert RGB 8:8:8 pixel data to realing point RGB
		m_Bitmap = new Color[m_Width * m_Height];
		real rec = 1.0f / 256;
		for ( int size = m_Width * m_Height, i = 0; i < size; i++ )
			m_Bitmap[i] = Color( t[i * 3 + 20] * rec, t[i * 3 + 19] * rec, t[i * 3 + 18] * rec );
		delete t;
	}
}

Color Texture::GetTexel( real a_U, real a_V )
{
	// fetch a bilinearly filtered texel
	real fu = (a_U + 1000.5f) * m_Width;
	real fv = (a_V + 1000.0f) * m_Width;
	int u1 = ((int)fu) % m_Width;
	int v1 = ((int)fv) % m_Height;
	int u2 = (u1 + 1) % m_Width;
	int v2 = (v1 + 1) % m_Height;
	// calculate fractional parts of u and v
	real fracu = fu - _floor( fu );
	real fracv = fv - _floor( fv );
	// calculate weight factors
	real w1 = (1 - fracu) * (1 - fracv);
	real w2 = fracu * (1 - fracv);
	real w3 = (1 - fracu) * fracv;
	real w4 = fracu *  fracv;
	// fetch four texels
	Color c1 = m_Bitmap[u1 + v1 * m_Width];
	Color c2 = m_Bitmap[u2 + v1 * m_Width];
	Color c3 = m_Bitmap[u1 + v2 * m_Width];
	Color c4 = m_Bitmap[u2 + v2 * m_Width];
	// scale and sum the four colors
	return c1 * w1 + c2 * w2 + c3 * w3 + c4 * w4;
}

// -----------------------------------------------------------
// Material class implementation
// -----------------------------------------------------------

Material::Material() :
	m_Color( Color( 0.2f, 0.2f, 0.2f ) ),
	m_Refl( 0 ), m_Diff( 0.2f ), m_Spec( 0.8f ), 
	m_RIndex( 1.5f ), m_DRefl( 0 ), m_Texture( 0 ),
	m_UScale( 1.0f ), m_VScale( 1.0f )
{
}

void Material::SetUVScale( real a_UScale, real a_VScale )
{ 
	m_UScale = a_UScale; 
	m_VScale = a_VScale; 
	m_RUScale = 1.0f / a_UScale;
	m_RVScale = 1.0f / a_VScale;
}

void Material::SetParameters( real a_Refl, real a_Refr, Color& a_Col, real a_Diff, real a_Spec )
{
	m_Refl = a_Refl;
	m_Refr = a_Refr;
	m_Color = a_Col;
	m_Diff = a_Diff;
	m_Spec = a_Spec;
}

// -----------------------------------------------------------
// Primitive methods
// -----------------------------------------------------------

Primitive::Primitive( int a_Type, vector3& a_Centre, real a_Radius )
{
	m_Centre = a_Centre;
	m_SqRadius = a_Radius * a_Radius;
	m_Radius = a_Radius;
	m_RRadius = 1.0f / a_Radius;
	m_Type = a_Type;
	m_Material = new Material();
	// set vectors for texture mapping
	m_Vn = vector3( 0, 1, 0 );
	m_Ve = vector3( 1, 0, 0 );
	m_Vc = m_Vn.Cross( m_Ve );
}

Primitive::Primitive( int a_Type, Vertex* a_V1, Vertex* a_V2, Vertex* a_V3 )
{
	m_Type = a_Type;
	m_Material = 0;
	m_Vertex[0] = a_V1;
	m_Vertex[1] = a_V2;
	m_Vertex[2] = a_V3;
	// init precomp
	vector3 A = m_Vertex[0]->GetPos();
	vector3 B = m_Vertex[1]->GetPos();
	vector3 C = m_Vertex[2]->GetPos();
	vector3 c = B - A;
	vector3 b = C - A;
	m_N = b.Cross( c );
	int u, v;
	if (_fabs( m_N.x ) > _fabs( m_N.y))
	{
		if (_fabs( m_N.x ) > _fabs( m_N.z )) k = 0; else k = 2;
	}
	else
	{
		if (_fabs( m_N.y ) > _fabs( m_N.z )) k = 1; else k = 2;
	}
	u = (k + 1) % 3;
	v = (k + 2) % 3;
	// precomp
	real krec = 1.0f / m_N.cell[k];
	nu = m_N.cell[u] * krec;
	nv = m_N.cell[v] * krec;
	nd = m_N.Dot( A ) * krec;
	// first line equation
	real reci = 1.0f / (b.cell[u] * c.cell[v] - b.cell[v] * c.cell[u]);
	bnu = b.cell[u] * reci;
	bnv = -b.cell[v] * reci;
	// second line equation
	cnu = c.cell[v] * reci;
	cnv = -c.cell[u] * reci;
	// finalize normal
	m_N.Normalize();
	m_Vertex[0]->SetNormal( m_N );
	m_Vertex[1]->SetNormal( m_N );
	m_Vertex[2]->SetNormal( m_N );
}

Primitive::~Primitive()
{
	if (m_Type == SPHERE) delete m_Material;
}

unsigned int modulo[] = { 0, 1, 2, 0, 1 };
int Primitive::Intersect( Ray& a_Ray, real& a_Dist )
{
	if (m_Type == SPHERE)
	{
		vector3 v = a_Ray.GetOrigin() - m_Centre;
		real b = -DOT( v, a_Ray.GetDirection() );
		real det = (b * b) - DOT( v, v ) + m_SqRadius;
		int retval = MISS;
		if (det > 0)
		{
			det = _sqrt( det );
			real i1 = b - det;
			real i2 = b + det;
			if (i2 > 0)
			{
				if (i1 < 0) 
				{
					if (i2 < a_Dist) 
					{
						a_Dist = i2;
						retval = INPRIM;
					}
				}
				else
				{
					if (i1 < a_Dist)
					{
						a_Dist = i1;
						retval = HIT;
					}
				}
			}
		}
		return retval;
	}
	else
	{
		#define ku modulo[k + 1]
		#define kv modulo[k + 2]
		vector3 O = a_Ray.GetOrigin(), D = a_Ray.GetDirection(), A = m_Vertex[0]->GetPos();
		const real lnd = 1.0f / (D.cell[k] + nu * D.cell[ku] + nv * D.cell[kv]);
		const real t = (nd - O.cell[k] - nu * O.cell[ku] - nv * O.cell[kv]) * lnd;
		if (!(a_Dist > t && t > 0)) return MISS;
		real hu = O.cell[ku] + t * D.cell[ku] - A.cell[ku];
		real hv = O.cell[kv] + t * D.cell[kv] - A.cell[kv];
		real beta = m_U = hv * bnu + hu * bnv;
		if (beta < 0) return MISS;
		real gamma = m_V = hu * cnu + hv * cnv;
		if (gamma < 0) return MISS;
		if ((m_U + m_V) > 1) return MISS;
		a_Dist = t;
		return (DOT( D, m_N ) > 0)?INPRIM:HIT;
	}
}

vector3 Primitive::GetNormal( vector3& a_Pos ) 
{ 
	if (m_Type == SPHERE) 
	{
		return (a_Pos - m_Centre) * m_RRadius; 
	}
	else 
	{
		vector3 N1 = m_Vertex[0]->GetNormal();
		vector3 N2 = m_Vertex[1]->GetNormal();
		vector3 N3 = m_Vertex[2]->GetNormal();
		vector3 N = N1 + m_U * (N2 - N1) + m_V * (N3 - N1);
		NORMALIZE( N );
		return N;
	}
}
	
Color Primitive::GetColor( vector3& a_Pos )
{
	Color retval;
	if (!m_Material->GetTexture()) retval = m_Material->GetColor(); else
	{
		if (m_Type == SPHERE)
		{
			vector3 vp = (a_Pos - m_Centre) * m_RRadius;
			real phi = _acos( -DOT( vp, m_Vn ) );
			real u, v = phi * m_Material->GetVScaleReci() * (1.0f / PI);
			real theta = (_acos( DOT( m_Ve, vp ) / _sin( phi ))) * (2.0f / PI);
			if (DOT( m_Vc, vp ) >= 0) u = (1.0f - theta) * m_Material->GetUScaleReci();
								 else u = theta * m_Material->GetUScaleReci();
			retval = m_Material->GetTexture()->GetTexel( u, v ) * m_Material->GetColor();
		}
		else
		{
			real U1 = m_Vertex[0]->GetU(), V1 = m_Vertex[0]->GetV();		
			real U2 = m_Vertex[1]->GetU(), V2 = m_Vertex[1]->GetV();		
			real U3 = m_Vertex[2]->GetU(), V3 = m_Vertex[2]->GetV();		
			real u = U1 + m_U * (U2 - U1) + m_V * (U3 - U1);
			real v = V1 + m_U * (V2 - V1) + m_V * (V3 - V1);
			retval = m_Material->GetTexture()->GetTexel( u, v ) * m_Material->GetColor();
		}
	}
	return retval;
}

#define FINDMINMAX( x0, x1, x2, min, max ) \
  min = max = x0; if(x1<min) min=x1; if(x1>max) max=x1; if(x2<min) min=x2; if(x2>max) max=x2;
// X-tests
#define AXISTEST_X01( a, b, fa, fb )											\
	p0 = a * v0.cell[1] - b * v0.cell[2], p2 = a * v2.cell[1] - b * v2.cell[2]; \
    if (p0 < p2) { min = p0; max = p2;} else { min = p2; max = p0; }			\
	rad = fa * a_BoxHalfsize.cell[1] + fb * a_BoxHalfsize.cell[2];				\
	if (min > rad || max < -rad) return 0;
#define AXISTEST_X2( a, b, fa, fb )												\
	p0 = a * v0.cell[1] - b * v0.cell[2], p1 = a * v1.cell[1] - b * v1.cell[2];	\
    if (p0 < p1) { min = p0; max = p1; } else { min = p1; max = p0;}			\
	rad = fa * a_BoxHalfsize.cell[1] + fb * a_BoxHalfsize.cell[2];				\
	if(min>rad || max<-rad) return 0;
// Y-tests
#define AXISTEST_Y02( a, b, fa, fb )											\
	p0 = -a * v0.cell[0] + b * v0.cell[2], p2 = -a * v2.cell[0] + b * v2.cell[2]; \
    if(p0 < p2) { min = p0; max = p2; } else { min = p2; max = p0; }			\
	rad = fa * a_BoxHalfsize.cell[0] + fb * a_BoxHalfsize.cell[2];				\
	if (min > rad || max < -rad) return 0;
#define AXISTEST_Y1( a, b, fa, fb )												\
	p0 = -a * v0.cell[0] + b * v0.cell[2], p1 = -a * v1.cell[0] + b * v1.cell[2]; \
    if (p0 < p1) { min = p0; max = p1; } else { min = p1; max = p0; }			\
	rad = fa * a_BoxHalfsize.cell[0] + fb * a_BoxHalfsize.cell[2];				\
	if (min > rad || max < -rad) return 0;
// Z-tests
#define AXISTEST_Z12( a, b, fa, fb )											\
	p1 = a * v1.cell[0] - b * v1.cell[1], p2 = a * v2.cell[0] - b * v2.cell[1]; \
    if(p2 < p1) { min = p2; max = p1; } else { min = p1; max = p2; }			\
	rad = fa * a_BoxHalfsize.cell[0] + fb * a_BoxHalfsize.cell[1];				\
	if (min > rad || max < -rad) return 0;
#define AXISTEST_Z0( a, b, fa, fb )												\
	p0 = a * v0.cell[0] - b * v0.cell[1], p1 = a * v1.cell[0] - b * v1.cell[1];	\
    if(p0 < p1) { min = p0; max = p1; } else { min = p1; max = p0; }			\
	rad = fa * a_BoxHalfsize.cell[0] + fb * a_BoxHalfsize.cell[1];				\
	if (min > rad || max < -rad) return 0;

bool Primitive::PlaneBoxOverlap( vector3& a_Normal, vector3& a_Vert, vector3& a_MaxBox )
{
	vector3 vmin, vmax;
	for( int q = 0; q < 3; q++ )
	{
		float v = a_Vert.cell[q];
		if (a_Normal.cell[q] > 0.0f)
		{
			vmin.cell[q] = -a_MaxBox.cell[q] - v;