main.c

完整的3D 模型检索程序

#include "glut.h"
#include <gl/gl.h>
#include <gl/glu.h>

#include <stdio.h>
#include <malloc.h>
#include <memory.h>
#include <time.h>
#include <float.h>
#include <string.h>
#include <limits.h>

#include "ds.h"
#include "RWObj.h"
#include "Bitmap.h"
#include "TranslateScale.h"
#include "Rotate.h"
#include "RegionShape.h"
#include "RecovAffine.h"
#include "Refine.h"
#include "edge.h"
#include "convert.h"
#include "ColorDescriptor.h"
#include "Circularity.h"
#include "FourierDescriptor.h"
#include "Eccentricity.h"

#define abs(a) (a>0)?(a):-(a)

#define	QUANT8				256		// 2^8
#define FD_SCALE			2		// *2 first, and then quantization
#define CIR_SCALE			2.318181818		// the range of circularity is [0~110], so *2.318 to be [0~255]
#define ECC_SCALE			25.5			// the range of circularity is [0~10], so *25.5 to be [0~255]

unsigned char	CamMap[CAMNUM_2]={0,1,2,3,4,5,6,7,8,9,5,6,7,8,9,2,3,4,0,1};

char srcfn[100];
char destfn[100];

int			winw = WIDTH, winh = HEIGHT;

pVer		vertex=NULL;
pTri		triangle=NULL;
int			NumVer=0, NumTri=0;		// total number of vertex and triangle.

pVer		vertex1, vertex2;
pTri		triangle1, triangle2;
int			NumVer1, NumTri1, NumVer2, NumTri2;		// total number of vertex and triangle.

// translate and scale of model 1 and 2
Ver				Translate1, Translate2;
double			Scale1, Scale2;

void FindCenter(unsigned char *srcBuff, int width, int height, double *CenX, double *CenY)
{
	int					x, y, count;
	unsigned char		*pImage;
	int					maxX, minX, maxY, minY;
	int					MeanX, MeanY; 

	count = 0;
	pImage = srcBuff;

	// ***********************************************************************************************
	// if use "mean" to each 2D shape independnetly, the origin will be moved a lot in 3D
	// if ues "center" to each 2D shape independnetly, the origin will be moved only a little in 3D
	// if center can be defined in 3D, the origin will not be moved any more.
	// But this will not very robust in 3D similarity transformation
	// In addition, to make center of each 2D shape more close to user drawn 2D shapes,
	// it's better to define center for each 2D shape independently

	// uee center of max and min to be center
	maxX = maxY = -1;
	minX = minY = INT_MAX;
	for (y=0 ; y<height ; y++)
	for (x=0 ; x<width; x++)
	{
		if( *pImage < 255 )
		{
			if( x > maxX ) maxX = x;
			if( x < minX ) minX = x;
			if( y > maxY ) maxY = y;
			if( y < minY ) minY = y;
		}
		pImage++;
	}

	if( maxX > 0 )
	{
		*CenX = (maxX+minX) / 2.0;
		*CenY = (maxY+minY) / 2.0;
	}
	else
		*CenX = *CenY = -1;		// nothing to be rendered

	// use mean to be center
/*	count = 0;
	MeanX = MeanY = 0;
	for (y=0 ; y<height ; y++)
	for (x=0 ; x<width; x++)
	{
		if( *pImage < 255 )
		{
			MeanX += x;
			MeanY += y;
			count ++;
		}
		pImage++;
	}

	if( count > 0 )
	{
		*CenX = (double)MeanX / (double)count;
		*CenY = (double)MeanY / (double)count;
	}
	else
		*CenX = *CenY = -1;		// nothing to be rendered
*/
}

void display(void)
{
	int				i, j;
	glClear (GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

	glPushMatrix();				
//		glColor3f((GLfloat)0.0, (GLfloat)0.0, (GLfloat)0.0);
		for(i = 0; i<NumTri; i++)
		{
			glColor3f((GLfloat)triangle[i].r, (GLfloat)triangle[i].g, (GLfloat)triangle[i].b);
			glBegin(GL_POLYGON);
				for(j=0; j<triangle[i].NodeName; j++)
					glVertex3d(vertex[triangle[i].v[j]].coor[0], vertex[triangle[i].v[j]].coor[1], vertex[triangle[i].v[j]].coor[2]);
			glEnd();
		}
	glPopMatrix();				

	glutSwapBuffers();
}

void RenderToMem(unsigned char *bmBits, unsigned char *bmColor, pVer CamVertex, pVer v, pTri t, int nv, int nt)
{
/*	
	// this is the same with reshape(), so this didn't do again
	glMatrixMode (GL_PROJECTION);
	glLoadIdentity ();
	gluPerspective(100.0, (GLfloat) winw/(GLfloat) winh, 0.1, 10.0);
//	glOrtho(-1, 1, -1, 1, 0.1, 100.0);
	glViewport (0, 0, (GLsizei) winw, (GLsizei) winh);
*/
	// test LookAt Up
/*	vector	v1, v2, v3;
	v1.x = CamVertex->coor[0];
	v1.y = CamVertex->coor[1];
	v1.z = CamVertex->coor[2];
	v2.x = 0;
	v2.y = 1;
	v2.z = 0;
	v3 = cross(v1, v2);*/

	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
//	gluLookAt(CAMSCALE*CamVertex->coor[0], CAMSCALE*CamVertex->coor[1], CAMSCALE*CamVertex->coor[2],
	gluLookAt(CamVertex->coor[0], CamVertex->coor[1], CamVertex->coor[2],
				0, 0, 0,
				0, 1, 0);
//				v3.x, v3.y, v3.z);

	vertex = v;
	triangle = t;
	NumVer = nv;
	NumTri = nt;
	display();

	glReadBuffer(GL_BACK);
//	glReadPixels(0, 0, winw, winh, GL_RED, GL_UNSIGNED_BYTE, bmBits);
	glReadPixels(0, 0, winw, winh, GL_DEPTH_COMPONENT, GL_UNSIGNED_BYTE, bmBits);
	if( bmColor )
		glReadPixels(0, 0, winw, winh, GL_RGB, GL_UNSIGNED_BYTE, bmColor);
}

void keyboard (unsigned char key, int x, int y)
{
	unsigned char	*srcBuff[CAMNUM], *destBuff[CAMNUM], *EdgeBuff, *ColorBuff[CAMNUM], *YuvBuff;
	char			filename[400];
	pVer			CamVertex[ANGLE];
	pTri			CamTriangle[ANGLE];
	int				CamNumVer[ANGLE], CamNumTri[ANGLE];		// total number of vertex and triangle.
	FILE			*fpt, *fpt1, *fpt2, *fpt3, *fpt4, *fpt_art_q8, *fpt_art_q4, *fpt_fd_q8, *fpt_fd, *fpt_cir_q8, *fpt_ecc_q8;//, *fpt_ccd;
	int				i, j, k, srcCam, destCam, p, r, a, itmp;
	double			cost[ANGLE][ANGLE][CAMNUM_2][CAMNUM_2];
	double			**matrix;
	static int		UseCam = 2;
	clock_t			start, finish;
	double			Err;
	// for region shape descriptor
	double			src_ArtCoeff[ANGLE][CAMNUM][ART_ANGULAR][ART_RADIAL];
	double			dest_ArtCoeff[ANGLE][CAMNUM][ART_ANGULAR][ART_RADIAL];
	double			MinErr, err;
	int				align[60][CAMNUM_2];
	unsigned char	q8_ArtCoeff[ANGLE][CAMNUM][ART_COEF];
	unsigned char	q4_ArtCoeff[ANGLE][CAMNUM][ART_COEF_2];
	// for color decsriptor
	unsigned __int64 CompactColor[ANGLE][CAMNUM];	// 63 bits for each image
	unsigned __int64 dest_CompactColor[ANGLE][CAMNUM];	// 63 bits for each image
	// for circularity
	double			cir_Coeff[ANGLE][CAMNUM];
	unsigned char	q8_cirCoeff[ANGLE][CAMNUM], dest_cirCoeff[ANGLE][CAMNUM];
	// for fourier descriptor
	double			src_FdCoeff[ANGLE][CAMNUM][FD_COEFF_NO], dest_FdCoeff[ANGLE][CAMNUM][FD_COEFF_NO];
	unsigned char	q8_FdCoeff[ANGLE][CAMNUM][FD_COEFF_NO];
	sPOINT			*Contour;
	unsigned char	*ContourMask;
	// for eccentricity
	double			ecc_Coeff[ANGLE][CAMNUM];
	unsigned char	q8_eccCoeff[ANGLE][CAMNUM], dest_eccCoeff[ANGLE][CAMNUM];
	// for compare
	int				Count, TopNum;
	pMatRes			pSearch, pmr, pmrr, pTop;
	// quantization version
	char			fname[400];
//	char			fn[200];
 	int				high, low, middle;
	double			QuantTable[17] = {	0.000000000, 0.003585473, 0.007418411, 0.011535520, 
										0.015982337, 0.020816302, 0.026111312, 0.031964674, 
										0.038508176, 0.045926586, 0.054490513, 0.064619488, 
										0.077016351, 0.092998687, 0.115524524, 0.154032694, 1.000000000};
	double			CenX[CAMNUM], CenY[CAMNUM];
	int				total;

	switch (key) 
	{
	case 27:
		exit(0);
		break;

// *************************************************************************************************
/*	case 'e':
		// initialize ART
		GenerateBasisLUT();
		// initialize: read camera set
		for(destCam=0; destCam<ANGLE; destCam++)
		{
			sprintf(filename, "12_%d", destCam);
			ReadObj(filename, CamVertex+destCam, CamTriangle+destCam, CamNumVer+destCam, CamNumTri+destCam);
		}

		if( (fpt1 = fopen("pair.txt", "r")) == NULL )
		{
			printf("pair.txt does not exist.\n");
			break;
		}

		while( fscanf(fpt1, "%s", srcfn) != EOF )
		{
			if( fscanf(fpt1, "%s", destfn) == EOF )
				break;

			// get the translatation and scale of the two model
			if( ReadObj(srcfn, &vertex1, &triangle1, &NumVer1, &NumTri1) == 0 )
				continue;

			if( ReadObj(destfn, &vertex2, &triangle2, &NumVer2, &NumTri2) == 0 )
				continue;

			// save original two models
			sprintf(filename, "%s_to_%s_org.obj", destfn, srcfn);
			SaveMergeObj(filename, vertex1, triangle1, NumVer1, NumTri1, vertex2, triangle2, NumVer2, NumTri2);

			// ****************************************************************
			// Corase alignment
			// ****************************************************************

			// record execute time --- start
			start = clock();

			// Translate and scale model 1 and 2
			TranslateScale(vertex2, NumVer2, triangle2, NumTri2, destfn, &Translate2, &Scale2);
			TranslateScale(vertex1, NumVer1, triangle1, NumTri1, srcfn, &Translate1, &Scale1);

			// read RED only, so size is winw*winh
			for(i=0; i<CAMNUM; i++)
				srcBuff[i] = (unsigned char *) malloc (winw * winh * sizeof(unsigned char));
			for(srcCam=0; srcCam<ANGLE; srcCam++)
			{
				// capture CAMNUM silhouette of srcfn to memory
				// 0.24 sec for an example (3012 triangles)
				for(i=0; i<CAMNUM; i++)
					RenderToMem(srcBuff[i], NULL, CamVertex[srcCam]+i, vertex1, triangle1, NumVer1, NumTri1);

				// 0.02 sec for an example
				FindRadius(srcBuff);

				// 0.16 sec for an example
				for(i=0; i<CAMNUM; i++)
					ExtractCoefficients(srcBuff[i], src_ArtCoeff[srcCam][i], NULL);
			}
			for(i=0; i<CAMNUM; i++)
				free(srcBuff[i]);

			// ********************************************************************************
			// capture CAMNUM silhouette of destfn to memory,
			// and get the cost between srcfn silhouette and destfn silhouette

			// read REB only, so size is winw*winh
			for(i=0; i<CAMNUM; i++)
				destBuff[i] = (unsigned char *) malloc (winw * winh * sizeof(unsigned char));
			for(destCam=0; destCam<ANGLE; destCam++)
			{
				// capture CAMNUM silhouette of dest to memory
				// 0.47 sec for an example (18110 triangles)
				for(i=0; i<CAMNUM; i++)
					RenderToMem(destBuff[i], NULL, CamVertex[destCam]+i, vertex2, triangle2, NumVer2, NumTri2);

				// 0.02 sec for an example
				FindRadius(destBuff);

				// 0.16 sec for an example
				for(i=0; i<CAMNUM; i++)
				{
					ExtractCoefficients(destBuff[i], dest_ArtCoeff[destCam][i], NULL);

					for(srcCam=0; srcCam<ANGLE; srcCam++)
						for(j=0; j<CAMNUM; j++)
							cost[srcCam][destCam][j][i] = GetDistance(dest_ArtCoeff[destCam][i], src_ArtCoeff[srcCam][j]);
				}
			}
			for(i=0; i<CAMNUM; i++)
				free(destBuff[i]);

			// ********************************************************************************
			// get the affine transformation from model 2 to model 1
			// initialize
			matrix = (double **) malloc (4 * sizeof(double *));
			for(i=0; i<4; i++)
			{
				matrix[i] = (double *) malloc(4 * sizeof(double));
				memset(matrix[i], 0, 4 * sizeof(double));
			}
			matrix[0][0] = matrix[1][1] = matrix[2][2] = matrix[3][3] = 1;

			// get the transformation matrix of model 2 based on model 1
			Err = RecoverAffine(matrix, cost, &UseCam);

			Rotate(vertex2, NumVer2, matrix);

			Scale1 = 1.0 / Scale1;
			Scale(vertex2, NumVer2, Scale1);

			Translate1.coor[0] = -Translate1.coor[0];
			Translate1.coor[1] = -Translate1.coor[1];
			Translate1.coor[2] = -Translate1.coor[2];
			Translate(vertex2, NumVer2, Translate1);

			sprintf(filename, "%s_to_%s.obj", destfn, srcfn);
			SaveObj(filename, vertex2, triangle2, NumVer2, NumTri2);

			// save two model
			// model 1 should translate and scale back to original
			Scale(vertex1, NumVer1, Scale1);
			Translate(vertex1, NumVer1, Translate1);
			sprintf(filename, "%s_to_%s_merge.obj", destfn, srcfn);
			SaveMergeObj(filename, vertex1, triangle1, NumVer1, NumTri1, vertex2, triangle2, NumVer2, NumTri2);

			for(i=0; i<4; i++)
				free(matrix[i]);
			free(matrix);

			// record execute time --- end
			finish = clock();
			fpt = fopen("time.txt", "a");
			fprintf(fpt, "%s to %s ( V: %d ; T: %d )\t: Corase: %.1f sec;", destfn, srcfn, NumVer1, NumTri1, (double)(finish - start) / CLOCKS_PER_SEC );
			fclose(fpt);

			// write Err to file
			fpt2 = fopen("Coarse_Err.txt", "a");
			fprintf(fpt2, "Err of %s to %s : %f\n", destfn, srcfn, Err);
			fclose(fpt2);

			// ****************************************************************
			// Refine alignment
			// ****************************************************************

			// record execute time --- start
			start = clock();