raytracer.cpp

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

		if (shade > 0)
		{
			// calculate diffuse shading
			if (prim->GetMaterial()->GetDiffuse() > 0)
			{
				real dot = DOT( L, N );
				if (dot > 0)
				{
					real diff = dot * prim->GetMaterial()->GetDiffuse() * shade;
					// add diffuse component to ray color
					a_Acc += diff * color * light->GetColor();
				}
			}
			// determine specular component using Schlick's BRDF approximation
			if (prim->GetMaterial()->GetSpecular() > 0)
			{
				// point light source: sample once for specular highlight
				vector3 R = L - 2.0f * DOT( L, N ) * N;
				real dot = DOT( a_Ray.GetDirection(), R );
				if (dot > 0)
				{
					real spec = dot * prim->GetMaterial()->GetSpecular() * shade / (50 - 50 * dot + dot);
					// add specular component to ray color
					a_Acc += spec * light->GetColor();
				}
			}
		}
	}
	// calculate reflection
	real refl = prim->GetMaterial()->GetReflection();
	if ((refl > 0.0f) && (a_Depth < TRACEDEPTH))
	{
		real drefl = prim->GetMaterial()->GetDiffuseRefl();
		if ((drefl > 0) && (a_Depth < 3))
		{
			// calculate diffuse reflection
			vector3 RP = a_Ray.GetDirection() - 2.0f * DOT( a_Ray.GetDirection(), N ) * N;
			vector3 RN1 = vector3( RP.z, RP.y, -RP.x );
			vector3 RN2 = RP.Cross( RN1 );
			refl *= a_SScale;
			for ( int i = 0; i < SAMPLES; i++ )
			{
				real xoffs, yoffs;
				do
				{
					xoffs = (m_Twister.Rand() - 0.5f) * drefl;
					yoffs = (m_Twister.Rand() - 0.5f) * drefl;
				}
				while ((xoffs * xoffs + yoffs * yoffs) > (drefl * drefl));
				vector3 R = RP + RN1 * xoffs + RN2 * yoffs * drefl;
				NORMALIZE( R );
				real dist;
				Color rcol( 0, 0, 0 );
				Raytrace( Ray( pi + R * EPSILON, R ), rcol, a_Depth + 1, a_RIndex, dist, a_Samples * 0.25f, a_SScale * 4 );
				m_RaysCast++;
				a_Acc += refl * rcol * color;
			}
		}
		else
		{
			// calculate perfect reflection
			vector3 N = prim->GetNormal( pi );
			vector3 R = a_Ray.GetDirection() - 2.0f * DOT( a_Ray.GetDirection(), N ) * N;
			Color rcol( 0, 0, 0 );
			real dist;
			Raytrace( Ray( pi + R * EPSILON, R ), rcol, a_Depth + 1, a_RIndex, dist, a_Samples * 0.5f, a_SScale * 2 );
			m_RaysCast++;
			a_Acc += refl * rcol * color;
		}
	}
	// calculate refraction
	real refr = prim->GetMaterial()->GetRefraction();
	if ((refr > 0) && (a_Depth < TRACEDEPTH))
	{
		real rindex = prim->GetMaterial()->GetRefrIndex();
		real n = a_RIndex / rindex;
		vector3 N = prim->GetNormal( pi ) * (real)result;
		real cosI = -DOT( N, a_Ray.GetDirection() );
		real cosT2 = 1.0f - n * n * (1.0f - cosI * cosI);
		if (cosT2 > 0.0f)
		{
			vector3 T = (n * a_Ray.GetDirection()) + (n * cosI - _sqrt( cosT2 )) * N;
			Color rcol( 0, 0, 0 );
			real dist;
			Raytrace( Ray( pi + T * EPSILON, T ), rcol, a_Depth + 1, rindex, dist, a_Samples * 0.5f, a_SScale * 2 );
			m_RaysCast++;
			// apply Beer's law
			Color absorbance = prim->GetMaterial()->GetColor() * 0.15f * -dist;
			Color transparency = Color( _exp( absorbance.r ), _exp( absorbance.g ), _exp( absorbance.b ) );
			a_Acc += rcol * transparency;
		}
	}
	// return pointer to primitive hit by primary ray
	return prim;
}

// -----------------------------------------------------------
// Engine::CalcShade
// Determines the light intensity received from a point light
// (in case of a SPHERE primitive) or an area light (in case
// of an AABB primitive)
// -----------------------------------------------------------
real Engine::CalcShade( Light* a_Light, vector3 a_IP, vector3& a_Dir, real a_Samples, real a_SScale )
{
	real retval;
	Primitive* prim = 0;
	if (a_Light->GetType() == Light::POINT)
	{
		// handle point light source
		retval = 0;
		a_Dir = a_Light->GetPos() - a_IP;
		real tdist = LENGTH( a_Dir );
		a_Dir *= (1.0f / tdist);
		tdist *= 1 - 4 * EPSILON;
		m_RaysCast++;
		if (!FindOccluder( Ray( a_IP + a_Dir * EPSILON, a_Dir ), tdist )) return 1;
	}
	else if (a_Light->GetType() == Light::AREA)
	{
		// Monte Carlo rendering
		retval = 0;
		a_Dir = a_Light->GetPos() - a_IP;
		NORMALIZE( a_Dir );
		vector3 deltax = a_Light->GetCellX(), deltay = a_Light->GetCellY();
		for ( int i = 0; i < a_Samples; i++ )
		{
			vector3 lp = a_Light->GetGrid( i & 15 ) + m_Twister.Rand() * deltax + m_Twister.Rand() * deltay;
			vector3 dir = lp - a_IP;
			real ldist = LENGTH( dir );
			dir *= 1.0f / ldist;
			ldist *= 1 - 4 * EPSILON;
			m_RaysCast++;
			if (!FindOccluder( Ray( a_IP + dir * EPSILON, dir ), ldist )) retval += a_SScale;
		}
	}
	return retval;
}

// -----------------------------------------------------------
// Engine::InitRender
// Initializes the renderer, by resetting the line / tile
// counters and precalculating some values
// -----------------------------------------------------------
void Engine::InitRender( vector3& a_Pos, vector3& a_Target )
{
	// set firts line to draw to
	m_CurrLine = 20;
	// set pixel buffer address of first pixel
	m_PPos = m_CurrLine * m_Width;
	// set eye and screen plane position
	m_Origin = vector3( 0, 0, -5 );
	m_P1 = vector3( -4,  3, 0 );
	m_P2 = vector3(  4,  3, 0 );
	m_P3 = vector3(  4, -3, 0 );
	m_P4 = vector3( -4, -3, 0 );
	// calculate camera matrix
	vector3 zaxis = a_Target - a_Pos;
	zaxis.Normalize();
	vector3 up( 0, 1, 0 );
	vector3 xaxis = up.Cross( zaxis );
	vector3 yaxis = xaxis.Cross( -zaxis );
	matrix m;
	m.cell[0] = xaxis.x, m.cell[1] = xaxis.y, m.cell[2] = xaxis.z;
	m.cell[4] = yaxis.x, m.cell[5] = yaxis.y, m.cell[6] = yaxis.z;
	m.cell[8] = zaxis.x, m.cell[9] = zaxis.y, m.cell[10] = zaxis.z;
	m.Invert();
	m.cell[3] = a_Pos.x, m.cell[7] = a_Pos.y, m.cell[11] = a_Pos.z;
	// move camera
	m_Origin = m.Transform( m_Origin );
	m_P1 = m.Transform( m_P1 );
	m_P2 = m.Transform( m_P2 );
	m_P3 = m.Transform( m_P3 );
	m_P4 = m.Transform( m_P4 );
	// calculate screen plane interpolation vectors
	m_DX = (m_P2 - m_P1) * (1.0f / m_Width);
	m_DY = (m_P4 - m_P1) * (1.0f / m_Height);
	// setup the tile renderer
	m_CurrCol = 0;
	m_CurrRow = 20 / TILESIZE;
	m_XTiles = m_Width / TILESIZE;
	m_YTiles = (m_Height - 40) / TILESIZE;
	// reset counters
	m_Intersections = 0;
	m_RaysCast = 0;
}

// -----------------------------------------------------------
// Engine::RenderRay
// Helper function, fires one ray in the regular grid
// -----------------------------------------------------------
Primitive* Engine::RenderRay( vector3 a_ScreenPos, Color& a_Acc )
{
	aabb e = m_Scene->GetExtends();
	vector3 dir = a_ScreenPos - m_Origin;
	NORMALIZE( dir );
	Color acc( 0, 0, 0 );
	Ray r( m_Origin, dir );
	m_RaysCast++;
	real dist;
	// trace ray
	return Raytrace( r, a_Acc, 1, 1.0f, dist, SAMPLES, 1.0f / SAMPLES );
}

// -----------------------------------------------------------
// Engine::Render
// Fires rays in the scene in a tile based fashion
// -----------------------------------------------------------
bool Engine::RenderTiles()
{
	// render scene in a tile based fashion
	aabb e = m_Scene->GetExtends();
	// initialize timer
	int msecs = GetTickCount();
	// render remaining tiles
	int tx = m_CurrCol, ty = m_CurrRow;
	int tdest = tx * TILESIZE + (ty * TILESIZE) * m_Width;
	vector3 tdir = m_P1 + (real)(tx * TILESIZE) * m_DX + (real)(ty * TILESIZE) * m_DY;
	while (1)
	{
		int dest = tdest;
		vector3 ldir = tdir;
		for ( int y = 0; y < TILESIZE; y++ )
		{
			vector3 pdir = ldir;
			for ( int x = 0; x < TILESIZE; x++ )
			{
				Color acc( 0, 0, 0 );
				Primitive* prim = RenderRay( pdir, acc );
				int red, green, blue;
				red = (int)(acc.r * 256);
				green = (int)(acc.g * 256);
				blue = (int)(acc.b * 256);
				if (red > 255) red = 255;
				if (green > 255) green = 255;
				if (blue > 255) blue = 255;
				m_Dest[dest++] = (red << 16) + (green << 8) + blue;
				pdir += m_DX;
			}
			ldir += m_DY;
			dest += (m_Width - TILESIZE);
		}
		tdest += TILESIZE;
		tdir += m_DX * TILESIZE;
		if (++tx == m_XTiles)
		{
			tx = 0;
			ty++;
			tdest = ty * TILESIZE * m_Width;
			tdir = m_P1 + (real)(ty * TILESIZE) * m_DY;
		}
		if (ty < m_YTiles)
		{
			if ((GetTickCount() - msecs) > 200) 
			{
				m_CurrCol = tx;
				m_CurrRow = ty;
				return false;
			}
		}
		else break;
	}
	return true;
}

}; // namespace Raytracer