metaballs.cpp

To review what a metaball is, and how to render them and to describe a new optimization I ve made

#include "metaballs.h"
#include "marchingcubes.h"
#include <memory.h>
#include "graphics.h"
#include <d3dx8.h>

//=============================================================================
CMetaballs::CMetaballs()
{
	m_fLevel    = 100.0f;
	m_nNumBalls = 20;

	m_nGridSize = 0;

	m_nMaxOpenVoxels = 32;
	m_pOpenVoxels    = new int[m_nMaxOpenVoxels*3];
   
	m_nNumOpenVoxels    = 0;
	m_pfGridEnergy      = 0;
	m_pnGridPointStatus = 0;
	m_pnGridVoxelStatus = 0;

	m_nNumVertices = 0;
	m_nNumIndices  = 0;
	m_pVertices    = 0;
	m_pIndices     = 0;

	srand(timeGetTime());

	for( int i = 0; i < m_nNumBalls; i++ )
	{
		m_Balls[i].p[0] = 0;
		m_Balls[i].p[1] = 0;
		m_Balls[i].p[2] = 0;
		m_Balls[i].v[0] = 0;
		m_Balls[i].v[1] = 0;
		m_Balls[i].v[2] = 0;
		m_Balls[i].a[0] = (float(rand())/RAND_MAX*2-1)/2;
		m_Balls[i].a[1] = (float(rand())/RAND_MAX*2-1)/2;
		m_Balls[i].a[2] = (float(rand())/RAND_MAX*2-1)/2;
		m_Balls[i].t = float(rand())/RAND_MAX;
		m_Balls[i].m = 1;
	}
}

//=============================================================================
void CMetaballs::Update(float dt)
{
	for( int i = 0; i < m_nNumBalls; i++ )
	{
		m_Balls[i].p[0] += dt*m_Balls[i].v[0];
		m_Balls[i].p[1] += dt*m_Balls[i].v[1];
		m_Balls[i].p[2] += dt*m_Balls[i].v[2];

		m_Balls[i].t -= dt;
		if( m_Balls[i].t < 0 )
		{
			// When is the next time to act?
			m_Balls[i].t = float(rand())/RAND_MAX;

			// Use a new attraction point
			m_Balls[i].a[0] = (float(rand())/RAND_MAX*2-1)/2;
			m_Balls[i].a[1] = (float(rand())/RAND_MAX*2-1)/2;
			m_Balls[i].a[2] = (float(rand())/RAND_MAX*2-1)/2;
		}

		// Accelerate towards the attraction point
		float x = m_Balls[i].a[0] - m_Balls[i].p[0];
		float y = m_Balls[i].a[1] - m_Balls[i].p[1];
		float z = m_Balls[i].a[2] - m_Balls[i].p[2];
		float fDist = 1/sqrtf(x*x + y*y + z*z);

		x *= fDist;
		y *= fDist;
		z *= fDist;
		
		m_Balls[i].v[0] += 0.1f*x*dt;
		m_Balls[i].v[1] += 0.1f*y*dt;
		m_Balls[i].v[2] += 0.1f*z*dt;

		fDist = m_Balls[i].v[0]*m_Balls[i].v[0] + 
		        m_Balls[i].v[1]*m_Balls[i].v[1] + 
		        m_Balls[i].v[2]*m_Balls[i].v[2];
		if( fDist > 0.040f )
		{
			fDist = 1/sqrtf(fDist);
			m_Balls[i].v[0] = 0.20f*m_Balls[i].v[0]*fDist;
			m_Balls[i].v[1] = 0.20f*m_Balls[i].v[1]*fDist;
			m_Balls[i].v[2] = 0.20f*m_Balls[i].v[2]*fDist;
		}

		if( m_Balls[i].p[0] < -1+m_fVoxelSize ) 
		{
			m_Balls[i].p[0] = -1+m_fVoxelSize;
			m_Balls[i].v[0] = 0;
		}
		if( m_Balls[i].p[0] >  1-m_fVoxelSize ) 
		{
			m_Balls[i].p[0] =  1-m_fVoxelSize;
			m_Balls[i].v[0] = 0;
		}
		if( m_Balls[i].p[1] < -1+m_fVoxelSize ) 
		{
			m_Balls[i].p[1] = -1+m_fVoxelSize;
			m_Balls[i].v[1] = 0;
		}
		if( m_Balls[i].p[1] >  1-m_fVoxelSize ) 
		{
			m_Balls[i].p[1] =  1-m_fVoxelSize;
			m_Balls[i].v[1] = 0;
		}
		if( m_Balls[i].p[2] < -1+m_fVoxelSize ) 
		{
			m_Balls[i].p[2] = -1+m_fVoxelSize;
			m_Balls[i].v[2] = 0;
		}
		if( m_Balls[i].p[2] >  1-m_fVoxelSize ) 
		{
			m_Balls[i].p[2] =  1-m_fVoxelSize;
			m_Balls[i].v[2] = 0;
		}
	}
}

//=============================================================================
void CMetaballs::Render()
{
	int nCase,x,y,z;
	bool bComputed;

	// Lock the vertex and index buffer
	CGraphics *pGfx = CGraphics::GetInstance();

	pGfx->LockVertexBuffer(MAX_VERTICES, (BYTE**)&m_pVertices);
	pGfx->LockIndexBuffer(MAX_INDICES, (BYTE**)&m_pIndices);
	m_nNumIndices = 0;
	m_nNumVertices = 0;

	// Clear status grids
	memset(m_pnGridPointStatus, 0, (m_nGridSize+1)*(m_nGridSize+1)*(m_nGridSize+1));
	memset(m_pnGridVoxelStatus, 0, m_nGridSize*m_nGridSize*m_nGridSize);

	for( int i = 0; i < m_nNumBalls; i++ )
	{
		x = ConvertWorldCoordinateToGridPoint(m_Balls[i].p[0]);
		y = ConvertWorldCoordinateToGridPoint(m_Balls[i].p[1]);
		z = ConvertWorldCoordinateToGridPoint(m_Balls[i].p[2]);

		// Work our way out from the center of the ball until the surface is
		// reached. If the voxel at the surface is already computed then this
		// ball share surface with a previous ball.
		bComputed = false;
		while(1)
		{
			if( IsGridVoxelComputed(x,y,z) )
			{
				bComputed = true;
				break;
			}

			nCase = ComputeGridVoxel(x,y,z);
			if( nCase < 255 )
				break;

			z--;
		}

		if( bComputed )
			continue;

		// Compute all voxels on the surface by computing neighbouring voxels
		// if the surface goes into them.
		AddNeighborsToList(nCase,x,y,z);

		while( m_nNumOpenVoxels )
		{
			m_nNumOpenVoxels--;
			x = m_pOpenVoxels[m_nNumOpenVoxels*3];
			y = m_pOpenVoxels[m_nNumOpenVoxels*3 + 1];
			z = m_pOpenVoxels[m_nNumOpenVoxels*3 + 2];

			nCase = ComputeGridVoxel(x,y,z);

			AddNeighborsToList(nCase,x,y,z);
		}
	}

	// Unlock buffers
	pGfx->UnlockVertexBuffer();
	pGfx->UnlockIndexBuffer();

	// Render the last triangles
	IDirect3DDevice8 *pDev = pGfx->GetD3DDevice();
	pDev->DrawIndexedPrimitive(D3DPT_TRIANGLELIST, 0, m_nNumVertices, 
		                                           0, m_nNumIndices/3);
}

//=============================================================================
void CMetaballs::AddNeighborsToList(int nCase, int x, int y, int z)
{
	if( CMarchingCubes::m_CubeNeighbors[nCase] & (1<<0) )
		AddNeighbor(x+1, y, z);

	if( CMarchingCubes::m_CubeNeighbors[nCase] & (1<<1) )
		AddNeighbor(x-1, y, z);

	if( CMarchingCubes::m_CubeNeighbors[nCase] & (1<<2) )
		AddNeighbor(x, y+1, z);

	if( CMarchingCubes::m_CubeNeighbors[nCase] & (1<<3) )
		AddNeighbor(x, y-1, z);

	if( CMarchingCubes::m_CubeNeighbors[nCase] & (1<<4) )
		AddNeighbor(x, y, z+1);

	if( CMarchingCubes::m_CubeNeighbors[nCase] & (1<<5) )
		AddNeighbor(x, y, z-1);
}

//=============================================================================
void CMetaballs::AddNeighbor(int x, int y, int z)
{
	if( IsGridVoxelComputed(x,y,z) || IsGridVoxelInList(x,y,z) )
		return;

	// Make sure the array is large enough
	if( m_nMaxOpenVoxels == m_nNumOpenVoxels )
	{
		m_nMaxOpenVoxels *= 2;
		int *pTmp = new int[m_nMaxOpenVoxels*3];
		memcpy(pTmp, m_pOpenVoxels, m_nNumOpenVoxels*3*sizeof(int));
		delete[] m_pOpenVoxels;
		m_pOpenVoxels = pTmp;
	}

	m_pOpenVoxels[m_nNumOpenVoxels*3] = x;
	m_pOpenVoxels[m_nNumOpenVoxels*3+1] = y;
	m_pOpenVoxels[m_nNumOpenVoxels*3+2] = z;

	SetGridVoxelInList(x,y,z);

	m_nNumOpenVoxels++;
}

//=============================================================================
float CMetaballs::ComputeEnergy(float x, float y, float z)
{
	float fEnergy = 0;
	float fSqDist;

	for( int i = 0; i < m_nNumBalls; i++ )
	{
		// The formula for the energy is 
		// 
		//   e += mass/distance^2 

		fSqDist = (m_Balls[i].p[0] - x)*(m_Balls[i].p[0] - x) +
		          (m_Balls[i].p[1] - y)*(m_Balls[i].p[1] - y) +
		          (m_Balls[i].p[2] - z)*(m_Balls[i].p[2] - z);

		if( fSqDist < 0.0001f ) fSqDist = 0.0001f;

		fEnergy += m_Balls[i].m / fSqDist;