bloomenthal.cpp

快速fft变换的c实现

// Bloomenthal.cpp: implementation of the Bloomenthal class. 
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
//#include "RBF3D.h"
#include "Bloomenthal.h"

/*
#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
*/

using namespace BLOOMENTHAL;

//static INTLISTS *cubetable[256];
static struct Bloomenthal::intlists *cubetable[256];

/*			edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
static int corner1[12]	   = {LBN,LTN,LBN,LBF,RBN,RTN,RBN,RBF,LBN,LBF,LTN,LTF};
static int corner2[12]	   = {LBF,LTF,LTN,LTF,RBF,RTF,RTN,RTF,RBN,RBF,RTN,RTF};
static int leftface[12]	   = {B,  L,  L,  F,  R,  T,  N,  R,  N,  B,  T,  F};
			     /* face on left when going corner1 to corner2 */
static int rightface[12]   = {L,  T,  N,  L,  B,  R,  R,  F,  B,  F,  N,  T};
			     /* face on right when going corner1 to corner2 */

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

Bloomenthal::Bloomenthal()
{

}

Bloomenthal::~Bloomenthal()
{

}


double Bloomenthal::torus(double x, double y, double z)
{
	double x2 = x*x, y2 = y*y, z2 = z*z;
    double a = x2+y2+z2+(0.5*0.5)-(0.1*0.1);
    return a*a-4.0*(0.5*0.5)*(y2+z2);
}

double Bloomenthal::sphere(double x, double y, double z)
{
	double rsq = x*x+y*y+z*z;
    return 1.0/(rsq < 0.00001? 0.00001 : rsq);
}


double Bloomenthal::blob(double x, double y, double z)
{
	return 4.0-sphere(x+1.0,y,z)-sphere(x,y+1.0,z)-sphere(x,y,z+1.0);
}


int Bloomenthal::triproc(int i1, int i2, int i3, VERTICES vertices)
{
	gvertices = vertices;
	gntris++;

	if(FILEOUT){
		fprintf(tri_file, "%d %d %d\n", i1, i2, i3);
	}
	else{
		TRILIST* n_tris = new TRILIST;
		n_tris->next = tris;
		n_tris->i1 = i1;
		n_tris->i2 = i2;
		n_tris->i3 = i3;
		tris = n_tris;
	}
    return 1;
}	

/**** An Implicit Surface Polygonizer ****/


/* polygonize: polygonize the implicit surface function
 *   arguments are:
 *	 double function (x, y, z)
 *		 double x, y, z (an arbitrary 3D point)
 *	     the implicit surface function
 *	     return negative for inside, positive for outside
 *	 double size
 *	     width of the partitioning cube
 *	 int bounds
 *	     max. range of cubes (+/- on the three axes) from first cube
 *	 double x, y, z
 *	     coordinates of a starting point on or near the surface
 *	     may be defaulted to 0., 0., 0.
 *	 int triproc (i1, i2, i3, vertices)
 *		 int i1, i2, i3 (indices into the vertex array)
 *		 VERTICES vertices (the vertex array, indexed from 0)
 *	     called for each triangle
 *	     the triangle coordinates are (for i = i1, i2, i3):
 *		 vertices.ptr[i].position.x, .y, and .z
 *	     vertices are ccw when viewed from the out (positive) side
 *		 in a left-handed coordinate system
 *	     vertex normals point outwards
 *	     return 1 to continue, 0 to abort
 *	 int mode
 *	     TET: decompose cube and polygonize six tetrahedra
 *	     NOTET: polygonize cube directly
 *   returns error or NULL
 */

int counter = 0;
char* Bloomenthal::polygonize(Mesh* mesh, double size, float bounds[6], double x, double y, double z, int mode)
{
	gntris = 0;
	if(FILEOUT){
		ver_file = fopen(VERFILE, "w");
		tri_file = fopen(TRIFILE, "w");
	}
	else{
		tris = new TRILIST;
		tris->i1 = -1;
		tris->i2 = -1;
		tris->i3 = -1;
		tris->next = NULL;
	}

	PROCESS p;
    int n, noabort;
    //CORNER *setcorner();
    TEST in, out;// find();

    //p.function = this->function;//function;
    //p.triproc = triproc;
    p.size = size;
    //p.bounds = bounds;
	p.boundsX1 = bounds[0];
	p.boundsX2 = bounds[1];
	p.boundsY1 = bounds[2];
	p.boundsY2 = bounds[3];
	p.boundsZ1 = bounds[4];
	p.boundsZ2 = bounds[5];
    p.delta = size/(double)(RES*RES);

    /* allocate hash tables and build cube polygon table: */
    p.centers = (CENTERLIST **) mycalloc(HASHSIZE,sizeof(CENTERLIST *));
    p.corners = (CORNERLIST **) mycalloc(HASHSIZE,sizeof(CORNERLIST *));
    p.edges =	(EDGELIST   **) mycalloc(2*HASHSIZE,sizeof(EDGELIST *));
    makecubetable();

    /* find point on surface, beginning search at (x, y, z): */
    srand(1);
    in = find1(1, &p, x, y, z);
    out = find1(0, &p, x, y, z);
    if (!in.ok || !out.ok) return "can't find starting point";
    converge(&in.p, &out.p, in.value, &p.start);

    /* push initial cube on stack: */
    p.cubes = (CUBES *) mycalloc(1, sizeof(CUBES)); /* list of 1 */
    p.cubes->cube.i = p.cubes->cube.j = p.cubes->cube.k = 0;
    p.cubes->next = NULL;

    /* set corners of initial cube: */
    for (n = 0; n < 8; n++)
	p.cubes->cube.corners[n] = setcorner1(&p, BIT(n,2), BIT(n,1), BIT(n,0));

    p.vertices.count = p.vertices.max = 0; /* no vertices yet */
    p.vertices.ptr = NULL;

    setcenter(p.centers, 0, 0, 0);

    while (p.cubes != NULL) { /* process active cubes till none left */
	CUBE c;
	CUBES *temp = p.cubes;
	c = p.cubes->cube;

	noabort = mode == TET?
	       /* either decompose into tetrahedra and polygonize: */
	       dotet(&c, LBN, LTN, RBN, LBF, &p) &&
	       dotet(&c, RTN, LTN, LBF, RBN, &p) &&
	       dotet(&c, RTN, LTN, LTF, LBF, &p) &&
	       dotet(&c, RTN, RBN, LBF, RBF, &p) &&
	       dotet(&c, RTN, LBF, LTF, RBF, &p) &&
	       dotet(&c, RTN, LTF, RTF, RBF, &p)
	       :
	       /* or polygonize the cube directly: */
	       docube(&c, &p);
	if (! noabort) return "aborted";

	/* pop current cube from stack */
	p.cubes = p.cubes->next;
	free((char *) temp);

	/* test six face directions, maybe add to stack: */
	testface(c.i-1, c.j, c.k, &c, L, LBN, LBF, LTN, LTF, &p);
	testface(c.i+1, c.j, c.k, &c, R, RBN, RBF, RTN, RTF, &p);
	testface(c.i, c.j-1, c.k, &c, B, LBN, LBF, RBN, RBF, &p);
	testface(c.i, c.j+1, c.k, &c, T, LTN, LTF, RTN, RTF, &p);
	testface(c.i, c.j, c.k-1, &c, N, LBN, LTN, RBN, RTN, &p);
	testface(c.i, c.j, c.k+1, &c, F, LBF, LTF, RBF, RTF, &p);

	for(n=0; n<8; n++)
		free((char *)c.corners[n]);
    }

	//free
	for(int i=0; i<HASHSIZE; i++){
		while(p.centers[i] != NULL){
			CENTERLIST* current = p.centers[i];
			p.centers[i] = current->next;
			free((char*)current);
		}
	}
	free((char*)p.centers);

	for(i=0; i<HASHSIZE; i++){
		while(p.corners[i] != NULL){
			CORNERLIST* current = p.corners[i];
			p.corners[i] = current->next;
			free((char*)current);
		}
	}
	free((char*)p.corners);

	for(i=0; i<2*HASHSIZE; i++){
		while(p.edges[i] != NULL){
			EDGELIST* current = p.edges[i];
			p.edges[i] = current->next;
			free((char*)current);
		}
	}
    free((char*)p.edges);

	//out
	int vertex_N = gvertices.count;
	int face_N = gntris;
	if(FILEOUT)
	{
		fclose(ver_file);
		fclose(tri_file);
		ver_file = fopen(VERFILE, "r");
		tri_file = fopen(TRIFILE, "r");
		
		FILE* mesh_file = fopen(MESHFILE, "w");
		fprintf(mesh_file, "%d\n%d\n", vertex_N, face_N);
		for(i=0; i<vertex_N; i++)
		{
			//float *v = mesh->vertex[i];
			//fscanf(ver_file, "%f %f %f", &v[0], &v[1], &v[2]);

			float vx, vy, vz;
			fscanf(ver_file, "%f %f %f", &vx, &vy, &vz);
			fprintf(mesh_file, "%f %f %f\n", vx, vy, vz);
		}
		for(i=0; i<face_N; i++)
		{
			//mesh->setPolygonCount(i, 3);
			//int* f = mesh->face[i];
			//fscanf(tri_file, "%d %d %d", &f[0], &f[1], &f[2]);

			int i0, i1, i2;
			fscanf(tri_file, "%d %d %d", &i0, &i1, &i2);
			fprintf(mesh_file, "3 %d %d %d\n", i0, i1, i2);
		}
		fclose(ver_file);
		fclose(tri_file);

		fclose(mesh_file);
	}
	else
	{
		mesh->SetFaceCount(face_N);
		mesh->SetVertexCount(vertex_N);
		for(i=0; i<vertex_N; i++)
		{
//			VERTEX v;
//			v = gvertices.ptr[i];
			float x = (float)gvertices.ptr[i].position.x;
			float y = (float)gvertices.ptr[i].position.y;
			float z = (float)gvertices.ptr[i].position.z;
			mesh->SetVertex(i, x, y, z);
//			mesh->vertex[i][0] = (float)v.position.x;
//			mesh->vertex[i][1] = (float)v.position.y;
//			mesh->vertex[i][2] = (float)v.position.z;
		}
		if (gvertices.ptr != NULL) free((char *)gvertices.ptr);

		TRILIST* current = tris;
		for(i=0; i<face_N; i++)
		{
			mesh->SetFace(i, current->i1, current->i2, current->i3);
//			mesh->setPolygonCount(i, 3);
//			mesh->face[i][0] = current->i1;
//			mesh->face[i][1] = current->i2;
//			mesh->face[i][2] = current->i3;
			TRILIST* previous = current;
			current = current->next;
			delete previous;
		}
	}

    return NULL;
}

/* testface: given cube at lattice (i, j, k), and four corners of face,
 * if surface crosses face, compute other four corners of adjacent cube
 * and add new cube to cube stack */

void Bloomenthal::testface(int i, int j, int k, CUBE *old, 
						   int face, int c1, int c2, int c3, int c4, PROCESS *p)
{
	CUBE new1;
    CUBES *oldcubes = p->cubes;
    //CORNER *setcorner();
    static int facebit[6] = {2, 2, 1, 1, 0, 0};
    int n, pos = old->corners[c1]->value > 0.0 ? 1 : 0, bit = facebit[face];

    /* test if no surface crossing, cube out of bounds, or already visited: */
    if ((old->corners[c2]->value > 0) == pos &&
	(old->corners[c3]->value > 0) == pos &&
	(old->corners[c4]->value > 0) == pos) return;
    //if (abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) return;
	float size = (float)p->size;
	if(i*size + p->start.x > p->boundsX1 
		|| i*size + p->start.x < p->boundsX2 
		|| j*size + p->start.y > p->boundsY1 
		|| j*size + p->start.y < p->boundsY2 
		|| k*size + p->start.z > p->boundsZ1 
		|| k*size + p->start.z < p->boundsZ2)
		return;
    if (setcenter(p->centers, i, j, k)) return;

    /* create new cube: */
    new1.i = i;
    new1.j = j;
    new1.k = k;
	/*
    for (n = 0; n < 8; n++) new1.corners[n] = NULL;
    new1.corners[FLIP(c1, bit)] = old->corners[c1];
    new1.corners[FLIP(c2, bit)] = old->corners[c2];
    new1.corners[FLIP(c3, bit)] = old->corners[c3];
    new1.corners[FLIP(c4, bit)] = old->corners[c4];*/
    for (n = 0; n < 8; n++)
	//if (new1.corners[n] == NULL)
	    new1.corners[n] = setcorner1(p, i+BIT(n,2), j+BIT(n,1), k+BIT(n,0));

    /*add cube to top of stack: */
    p->cubes = (CUBES *) mycalloc(1, sizeof(CUBES));
    p->cubes->cube = new1;
    p->cubes->next = oldcubes;
}

/* setcorner: return corner with the given lattice location
   set (and cache) its function value */
Bloomenthal::CORNER * Bloomenthal::setcorner1(PROCESS *p, int i, int j, int k)
{
	/* for speed, do corner value caching here */
    CORNER *c = (CORNER *) mycalloc(1, sizeof(CORNER));
    int index = HASH(i, j, k);
    CORNERLIST *l = p->corners[index];
    c->i = i; c->x = p->start.x+((double)i-.5)*p->size;
    c->j = j; c->y = p->start.y+((double)j-.5)*p->size;
    c->k = k; c->z = p->start.z+((double)k-.5)*p->size;
    for (; l != NULL; l = l->next)