sph_surf.c

游戏编程精粹６第中关于粒子的实时流体仿真系统,对入门的游戏开发者很有帮助．

{
    v->x = x;
    v->y = y;
    v->z = z;
}


GLfloat ABS(GLfloat x)
{
    return (x > 0 ? x : -x);
}



_inline void draw_tri(GLfloat* n0, GLfloat* p0,
					  GLfloat* n1, GLfloat* p1,
					  GLfloat* n2, GLfloat* p2)
{
	glNormal3fv(n0);
	glVertex3fv(p0);
	glNormal3fv(n1);
	glVertex3fv(p1);
	glNormal3fv(n2);
	glVertex3fv(p2);
}



#ifndef FORCE_FLAT
void lerp_vertex(vector3f* p, vector3f* n,
                 cellvertex* v0, cellvertex* v1, GLfloat isolevel)
{
    GLfloat mu;

    if (ABS(isolevel - v0->val) < EPSILON)
    {
        *n = v0->n;
        *p = v0->p;
        return;
    }

    if (ABS(isolevel - v1->val) < EPSILON)
    {
        *n = v1->n;
        *p = v1->p;
        return;
    }
    
    if (ABS(v0->val - v1->val) < EPSILON)
    {
        *n = v0->n;
        *p = v0->p;
        return;
    }

    mu = (isolevel - v0->val) / (v1->val - v0->val);

    p->x = v0->p.x + mu * (v1->p.x - v0->p.x);
    p->y = v0->p.y + mu * (v1->p.y - v0->p.y);
    p->z = v0->p.z + mu * (v1->p.z - v0->p.z);

    n->x = v0->n.x + mu * (v1->n.x - v0->n.x);
    n->y = v0->n.y + mu * (v1->n.y - v0->n.y);
    n->z = v0->n.z + mu * (v1->n.z - v0->n.z);
}

#else

void lerp_vertex(vector3f* p, vector3f* n,
                 cellvertex* v0, cellvertex* v1, GLfloat isolevel)
{
    GLfloat mu;
	GLfloat mag;
	GLfloat inv_mag;

    if (ABS(isolevel - v0->val) < EPSILON)
    {
        *n = v0->n;
        *p = v0->p;
        return;
    }

    
    if (ABS(isolevel - v1->val) < EPSILON)
    {
        *n = v1->n;
        *p = v1->p;
        return;
    }
    
    if (ABS(v0->val - v1->val) < EPSILON)
    {
        *n = v0->n;
        *p = v0->p;
        return;
    }

    mu = (isolevel - v0->val) / (v1->val - v0->val);

    p->x = v0->p.x + mu * (v1->p.x - v0->p.x);
    p->y = v0->p.y + mu * (v1->p.y - v0->p.y);
    p->z = v0->p.z + mu * (v1->p.z - v0->p.z);


    n->x = v0->n.x + mu * (v1->n.x - v0->n.x);
    n->y = v0->n.y + mu * (v1->n.y - v0->n.y);
    n->z = v0->n.z + mu * (v1->n.z - v0->n.z);

	mag = glg_sqrt(n->x*n->x + n->y*n->y + n->z*n->z);
	inv_mag = 1.0f/mag;

	n->x *= inv_mag;
	n->y *= inv_mag;
	n->z *= inv_mag;

	if (n->z > (0.9f))
	{
		n->x = 0.0f;
		n->y = 0.0f;
		n->z = 1.0f;
	}
}
#endif


void glLerpVertex(cellvertex* v0, cellvertex* v1, GLfloat isolevel)
{
	vector3f p;
	vector3f n;

	lerp_vertex(&p, &n, v0, v1, isolevel);

	glNormal3fv((GLfloat*)&n);
	glVertex3fv((GLfloat*)&p);
}



void polygonize(gridcell* grid, GLfloat isolevel)
{
    vector3f vertlist[12];
    vector3f normlist[12];
    int i;
    int cubeindex = 0;
    int edgeFlags = 0;

    if (grid->v[0]->val < isolevel) cubeindex |= 1;
    if (grid->v[1]->val < isolevel) cubeindex |= 2;
    if (grid->v[2]->val < isolevel) cubeindex |= 4;
    if (grid->v[3]->val < isolevel) cubeindex |= 8;
    if (grid->v[4]->val < isolevel) cubeindex |= 16;
    if (grid->v[5]->val < isolevel) cubeindex |= 32;
    if (grid->v[6]->val < isolevel) cubeindex |= 64;
    if (grid->v[7]->val < isolevel) cubeindex |= 128;

    edgeFlags = edgeTable[cubeindex];

    if (edgeFlags == 0)
        return;

   if (edgeFlags & 1)
       lerp_vertex(&vertlist[0], &normlist[0],
                   grid->v[0], grid->v[1], isolevel);
   if (edgeFlags & 2)
       lerp_vertex(&vertlist[1], &normlist[1],
                   grid->v[1], grid->v[2], isolevel);       
   if (edgeFlags & 4)
       lerp_vertex(&vertlist[2], &normlist[2],
                   grid->v[2], grid->v[3], isolevel);       
   if (edgeFlags & 8)
       lerp_vertex(&vertlist[3], &normlist[3],
                   grid->v[3], grid->v[0], isolevel);       
   if (edgeFlags & 16)
       lerp_vertex(&vertlist[4], &normlist[4],
                   grid->v[4], grid->v[5], isolevel);       

   if (edgeFlags & 32)
       lerp_vertex(&vertlist[5], &normlist[5],
                   grid->v[5], grid->v[6], isolevel);       

   if (edgeFlags & 64)
       lerp_vertex(&vertlist[6], &normlist[6],
                   grid->v[6], grid->v[7], isolevel);       

   if (edgeFlags & 128)
       lerp_vertex(&vertlist[7], &normlist[7],
                   grid->v[7], grid->v[4], isolevel);       

   if (edgeFlags & 256)
       lerp_vertex(&vertlist[8], &normlist[8],
                   grid->v[0], grid->v[4], isolevel);       

   if (edgeFlags & 512)
       lerp_vertex(&vertlist[9], &normlist[9],
                   grid->v[1], grid->v[5], isolevel);       

   if (edgeFlags & 1024)
       lerp_vertex(&vertlist[10], &normlist[10],
                   grid->v[2], grid->v[6], isolevel);       

   if (edgeFlags & 2048)
       lerp_vertex(&vertlist[11], &normlist[11],
                   grid->v[3], grid->v[7], isolevel);       


    for (i = 0; triTable[cubeindex][i] != -1; i += 3)
    {
        glNormal3fv((GLfloat*)&normlist[triTable[cubeindex][i]]);
        glVertex3fv((GLfloat*)&vertlist[triTable[cubeindex][i]]);

        glNormal3fv((GLfloat*)&normlist[triTable[cubeindex][i + 2]]);
        glVertex3fv((GLfloat*)&vertlist[triTable[cubeindex][i + 2]]);

        glNormal3fv((GLfloat*)&normlist[triTable[cubeindex][i + 1]]);
        glVertex3fv((GLfloat*)&vertlist[triTable[cubeindex][i + 1]]);

    }
}


void set_cellvertex(cellvertex* v, 
                    GLfloat* density, GLfloat stride,
                    GLuint width, GLuint height,
                    GLuint xn, GLuint yn, GLuint zn)
{
    GLuint wh = width*height;
    GLuint base = xn + yn*width + zn*width*height;
    v->n.x = -(density[base + 1] - density[base - 1]);
    v->n.y = -(density[base + width] - density[base - width]);
    v->n.z = -(density[base + wh] - density[base - wh]);

    v->val = density[base];
    v->p.x = xn*stride;
    v->p.y = yn*stride;
    v->p.z = zn*stride;
}

int clamp(int x, int mi, int ma)
{
	if (x <= mi)
		return mi;
	if (x >= ma)
		return ma;
	return x;
}

void render_marching_cubes(const implicit_surface* v,
						   int x0, int x1,
						   int y0, int y1,
						   int z0, int z1,
						   float threshold)
{
	int i;
	int xn;
	int yn;
	int zn;

	int xn1;
	int yn1;
	int zn1;
	int width;
	int height;
	int depth;
	GLfloat stride;
	GLfloat* density;
	gridcell cell;

	density = v->density;
	width = v->width;
	height = v->height;
	depth = v->depth;
	stride = v->voxelsize;

	x0 = clamp(x0, 0, width - 1);
	x1 = clamp(x1, 0, width - 1);
	y0 = clamp(y0, 0, height - 1);
	y1 = clamp(y1, 0, height - 1);
	z0 = clamp(z0, 0, depth - 1);
	z1 = clamp(z1, 0, depth - 1);


    glBegin(GL_TRIANGLES);
	for (i = 0; i < 8; i++)
		cell.v[i] = &cell.vmem[i];

	for (zn = z0 + 1; zn <= z1 - 2; zn++)
		for (yn = y0 + 1; yn <= y1 - 2; yn++)
		{
			xn = x0 + 1;
			yn1 = yn + 1;
			zn1 = zn + 1;
			set_cellvertex(cell.v[0], density, stride, width, height,
				xn, yn, zn);

			set_cellvertex(cell.v[1], density,	stride, width, height,
				xn, yn1, zn);

			set_cellvertex(cell.v[4], density,	stride, width, height,
				xn, yn, zn1);

			set_cellvertex(cell.v[5], density,	stride, width, height,
				xn, yn1, zn1);

			for (xn = x0 + 1; xn <= x1 - 2; xn++)
			{
				xn1 = xn + 1;

                set_cellvertex(cell.v[2], density,
                               stride, width, height,
                               xn1, yn1, zn);

                set_cellvertex(cell.v[3], density,
                               stride, width, height,
                               xn1, yn, zn);


                set_cellvertex(cell.v[6], density,
                               stride, width, height,
                               xn1, yn1, zn1);

                set_cellvertex(cell.v[7], density,
                               stride, width, height,
                               xn1, yn, zn1);

                polygonize(&cell, threshold);
				SWAP_V(cell.v[0], cell.v[3]);
				SWAP_V(cell.v[1], cell.v[2]);
				SWAP_V(cell.v[4], cell.v[7]);
				SWAP_V(cell.v[5], cell.v[6]);
            }
		}

		glEnd();
}