raytracer.cpp

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

// -----------------------------------------------------------
// raytracer.cpp
// 2004 - Jacco Bikker - [e-mail] - www.bik5.com -   <><
// -----------------------------------------------------------
// LATEST PERFORMANCE: 512x384, shark, 256 balls, 1.232 seconds
// (Intel compiler, full power)
// -----------------------------------------------------------

#include "raytracer.h"
#include "scene.h"
#include "common.h"
#include "windows.h"
#include "winbase.h"
#include "memory.h"

namespace Raytracer {

Ray::Ray( vector3& a_Origin, vector3& a_Dir ) : 
	m_Origin( a_Origin ), 
	m_Direction( a_Dir )
{
}

Engine::Engine()
{
	m_Scene = new Scene();
	KdTree::SetMemoryManager( new MManager() );
	m_Mod = new int[64];
	m_Mod = (int*)((((unsigned long)m_Mod) + 32) & (0xffffffff - 31));
	m_Mod[0] = 0, m_Mod[1] = 1, m_Mod[2] = 2, m_Mod[3] = 0, m_Mod[4] = 1;
	m_Stack = new kdstack[64];
	m_Stack = (kdstack*)((((unsigned long)m_Stack) + 32) & (0xffffffff - 31));
}

Engine::~Engine()
{
	delete m_Scene;
}

// -----------------------------------------------------------
// Engine::SetTarget
// Sets the render target canvas
// -----------------------------------------------------------
void Engine::SetTarget( Pixel* a_Dest, int a_Width, int a_Height )
{
	// set pixel buffer address & size
	m_Dest = a_Dest;
	m_Width = a_Width;
	m_Height = a_Height;
}

// -----------------------------------------------------------
// Engine::FindNearest
// Finds the nearest intersection in a KdTree for a ray
// -----------------------------------------------------------
int Engine::FindNearest( Ray& a_Ray, real& a_Dist, Primitive*& a_Prim )
{
	real tnear = 0, tfar = a_Dist, t;
	int retval = 0;
	vector3 p1 = m_Scene->GetExtends().GetPos();
	vector3 p2 = p1 + m_Scene->GetExtends().GetSize();
	vector3 D = a_Ray.GetDirection(), O = a_Ray.GetOrigin();
	for ( int i = 0; i < 3; i++ ) if (D.cell[i] < 0) 
	{
		if (O.cell[i] < p1.cell[i]) return 0;
	}
	else if (O.cell[i] > p2.cell[i]) return 0;
	// clip ray segment to box
	for ( i = 0; i < 3; i++ )
	{
		real pos = O.cell[i] + tfar * D.cell[i];
		if (D.cell[i] < 0)
		{
			// clip end point
			if (pos < p1.cell[i]) tfar = tnear + (tfar - tnear) * ((O.cell[i] - p1.cell[i]) / (O.cell[i] - pos));
			// clip start point
			if (O.cell[i] > p2.cell[i]) tnear += (tfar - tnear) * ((O.cell[i] - p2.cell[i]) / (tfar * D.cell[i]));
		}
		else
		{
			// clip end point
			if (pos > p2.cell[i]) tfar = tnear + (tfar - tnear) * ((p2.cell[i] - O.cell[i]) / (pos - O.cell[i]));
			// clip start point
			if (O.cell[i] < p1.cell[i]) tnear += (tfar - tnear) * ((p1.cell[i] - O.cell[i]) / (tfar * D.cell[i]));
		}
		if (tnear > tfar) return 0;
	}
	// init stack
	int entrypoint = 0, exitpoint = 1;
	// init traversal
	KdTreeNode* farchild, *currnode;
	currnode = m_Scene->GetKdTree()->GetRoot();
	m_Stack[entrypoint].t = tnear;
	if (tnear > 0.0f) m_Stack[entrypoint].pb = O + D * tnear;
				 else m_Stack[entrypoint].pb = O;
	m_Stack[exitpoint].t = tfar;
	m_Stack[exitpoint].pb = O + D * tfar;
	m_Stack[exitpoint].node = 0;
	// traverse kd-tree
	while (currnode)
	{
		while (!currnode->IsLeaf())
		{
			real splitpos = currnode->GetSplitPos();
			int axis = currnode->GetAxis();
			if (m_Stack[entrypoint].pb.cell[axis] <= splitpos)
			{
				if (m_Stack[exitpoint].pb.cell[axis] <= splitpos)
				{
					currnode = currnode->GetLeft();
					continue;
				}
				if (m_Stack[exitpoint].pb.cell[axis] == splitpos)
				{
					currnode = currnode->GetRight();
					continue;
				}
				currnode = currnode->GetLeft();
				farchild = currnode + 1; // GetRight();
			}
			else
			{
				if (m_Stack[exitpoint].pb.cell[axis] > splitpos)
				{
					currnode = currnode->GetRight();
					continue;
				}
				farchild = currnode->GetLeft();
				currnode = farchild + 1; // GetRight();
			}
			t = (splitpos - O.cell[axis]) / D.cell[axis];
			int tmp = exitpoint++;
			if (exitpoint == entrypoint) exitpoint++;
			m_Stack[exitpoint].prev = tmp;
			m_Stack[exitpoint].t = t;
			m_Stack[exitpoint].node = farchild;
			m_Stack[exitpoint].pb.cell[axis] = splitpos;
			int nextaxis = m_Mod[axis + 1];
			int prevaxis = m_Mod[axis + 2];
			m_Stack[exitpoint].pb.cell[nextaxis] = O.cell[nextaxis] + t * D.cell[nextaxis];
			m_Stack[exitpoint].pb.cell[prevaxis] = O.cell[prevaxis] + t * D.cell[prevaxis];
		}
		ObjectList* list = currnode->GetList();
		real dist = m_Stack[exitpoint].t;
		while (list)
		{
			Primitive* pr = list->GetPrimitive();
			int result;
			m_Intersections++;
			if (result = pr->Intersect( a_Ray, dist ))
			{
				retval = result;
				a_Dist = dist;
				a_Prim = pr;
			}
			list = list->GetNext();
		}
		if (retval) return retval;
		entrypoint = exitpoint;
		currnode = m_Stack[exitpoint].node;
		exitpoint = m_Stack[entrypoint].prev;
	}
	return 0;
}

// -----------------------------------------------------------
// Engine::FindOccluder
// Finds any occluder between the origin and a light source
// -----------------------------------------------------------
int Engine::FindOccluder( Ray& a_Ray, real& a_Dist )
{
	real tnear = EPSILON, t;
	vector3 O, D = a_Ray.GetDirection();
	// init stack
	int entrypoint = 0, exitpoint = 1;
	// init traversal
	KdTreeNode* farchild, *currnode = m_Scene->GetKdTree()->GetRoot();
	m_Stack[entrypoint].t = tnear;
	m_Stack[entrypoint].pb = O = a_Ray.GetOrigin();
	m_Stack[exitpoint].t = a_Dist;
	m_Stack[exitpoint].pb = O + D * a_Dist;
	m_Stack[exitpoint].node = 0;
	// traverse kd-tree
	while (currnode)
	{
		while (!currnode->IsLeaf())
		{
			real splitpos = currnode->GetSplitPos();
			int axis = currnode->GetAxis();
			if (m_Stack[entrypoint].pb.cell[axis] <= splitpos)
			{
				if (m_Stack[exitpoint].pb.cell[axis] <= splitpos)
				{
					currnode = currnode->GetLeft();
					continue;
				}
				if (m_Stack[exitpoint].pb.cell[axis] == splitpos)
				{
					currnode = currnode->GetRight();
					continue;
				}
				currnode = currnode->GetLeft();
				farchild = currnode + 1; // GetRight();
			}
			else
			{
				if (m_Stack[exitpoint].pb.cell[axis] > splitpos)
				{
					currnode = currnode->GetRight();
					continue;
				}
				farchild = currnode->GetLeft();
				currnode = farchild + 1; // GetRight();
			}
			t = (splitpos - O.cell[axis]) / D.cell[axis];
			int tmp = exitpoint;
			if (++exitpoint == entrypoint) exitpoint++;
			m_Stack[exitpoint].prev = tmp;
			m_Stack[exitpoint].t = t;
			m_Stack[exitpoint].node = farchild;
			m_Stack[exitpoint].pb.cell[axis] = splitpos;
			int nextaxis = m_Mod[axis + 1];
			int prevaxis = m_Mod[axis + 2];
			m_Stack[exitpoint].pb.cell[nextaxis] = O.cell[nextaxis] + t * D.cell[nextaxis];
			m_Stack[exitpoint].pb.cell[prevaxis] = O.cell[prevaxis] + t * D.cell[prevaxis];
		}
		ObjectList* list = currnode->GetList();
		real dist = a_Dist; // m_Stack[exitpoint].t;
		while (list)
		{
			m_Intersections++;
			if (list->GetPrimitive()->Intersect( a_Ray, dist )) return 1;
			list = list->GetNext();
		}
		entrypoint = exitpoint;
		currnode = m_Stack[exitpoint].node;
		exitpoint = m_Stack[entrypoint].prev;
	}
	return 0;
}

// -----------------------------------------------------------
// Engine::Raytrace
// Naive ray tracing: Intersects the ray with every primitive
// in the scene to determine the closest intersection
// -----------------------------------------------------------
Primitive* Engine::Raytrace( Ray& a_Ray, Color& a_Acc, int a_Depth, real a_RIndex, real& a_Dist, real a_Samples, real a_SScale )
{
	// trace primary ray
	a_Dist = 10000.0f;
	Primitive* prim = 0;
	int result;
	// find the nearest intersection
	if (!(result = FindNearest( a_Ray, a_Dist, prim ))) return 0;
	// determine color at point of intersection
	vector3 pi = a_Ray.GetOrigin() + a_Ray.GetDirection() * a_Dist;
	Color color = prim->GetColor( pi );
	vector3 N = prim->GetNormal( pi );
	// trace lights
	for ( int l = 0; l < m_Scene->GetNrLights(); l++ )
	{
		Light* light = m_Scene->GetLight( l );
		// handle point light source
		vector3 L;
		real shade = CalcShade( light, pi, L, a_Samples, a_SScale );