dlgfittingprocess.cpp

Hi guys, I have a B-spline curve algorithm by using Genetic algorithm. Any interested?

// dlgFittingProcess.cpp : implementation file
#include "stdafx.h"
#include "GA.h"
#include "dlgFittingProcess.h"
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "matrix.h"
#include <time.h>
#include "cubicBSplSurf.h"
#include "dlgPlotting.h"

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

#define data 400
#define Draw 1000
#define numtemp 10
// global variable
	int countKnot[102],countKnotb[102],
	int countTemp;
	float knotgen[2001];
	int dataspace = 8;// for data showing
	int meshnumber;
	float tknotdata[3000];
	
	extern float C1[4000],CT[4000],MultiMatrix[4000],Dx[data],Dy[data],Dz[data],B1[4000];
	extern float Px[data],Py[data],Pz[data],B2x[data],B2y[data],B2z[data];
	
	extern float pocx0[data],pocy0[data],pocx1[data],pocy1[data];
	extern int countPoc;

	extern int countBest,countBestb;
	extern float nbasis2[100];
	extern float Nmatrix[4000],Nmatrixb[2000],Nmatrix2[4000];
	extern float C[data][data];
	extern float RMS[101],pm[101],tknot[100],tknotb[100];
	extern int countY,countX;
	extern int countYd;
	extern int begincountY,begincountX;
	extern int space;
	extern int scale;

	extern const int number;
	extern float x[data],y[data],z[data];
	extern	int stepTCount;
	extern	float random[101][100],randomb[101][100],IndiCrossBest[2][100],IndiCrossBestb[2][100],RMSMinCurve[50],IndiCrossBestCurve[50][100];
	extern	int area,count;
	extern	float valueT;
	extern	float RMSMin,Pxfix[data],Pyfix[data],Pzfix[data];
	extern int IndividualLoop,generation,times;
	extern float DrawPxfix[Draw],DrawPyfix[Draw],DrawPzfix[Draw],DrawB2x[Draw],DrawB2y[Draw],DrawB2z[Draw];
/////////////////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////////////////
// CdlgFittingProcess dialog


CdlgFittingProcess::CdlgFittingProcess(CWnd* pParent /*=NULL*/)
	: CDialog(CdlgFittingProcess::IDD, pParent)
{
	//{{AFX_DATA_INIT(CdlgFittingProcess)
	m_population = IndividualLoop ;
	m_length = 5;
	//}}AFX_DATA_INIT
}


void CdlgFittingProcess::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CdlgFittingProcess)
	DDX_Text(pDX, IDC_EDIT1, m_population);
	DDX_Text(pDX, IDC_EDIT2, m_length);
	DDX_Control(pDX, IDC_MSHFLEXGRID2, m_grid2);
	//}}AFX_DATA_MAP
}


BEGIN_MESSAGE_MAP(CdlgFittingProcess, CDialog)
	//{{AFX_MSG_MAP(CdlgFittingProcess)
		// NOTE: the ClassWizard will add message map macros here
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CdlgFittingProcess message handlers
//----------------------------------------------------------------------*
//Set individuals table
//----------------------------------------------------------------------*
void CdlgFittingProcess::Refresh()
{
	char strNum[11];
	m_grid2.Clear();
	m_grid2.SetTextMatrix(0,0,"Individuals");
	for (int i = 1; i < 8; i++) {
		sprintf(strNum,"%d",i);		
		m_grid2.SetTextMatrix(0,i,strNum);
	}

	for (i = 1; i < 11; i++) {
		sprintf(strNum,"%d",i);		
		m_grid2.SetTextMatrix(i,0,strNum);
	}

	for (i = 1; i < 11; i++) {
		sprintf(strNum,"%d",0);		
		m_grid2.SetTextMatrix(i,1,strNum);

	}

	for (i = 1; i < 11; i++) {
		sprintf(strNum,"%d",1);		
		m_grid2.SetTextMatrix(i,7,strNum);

	}
}
/////////////////////////////////////////////////////////////////////////////


BOOL CdlgFittingProcess::OnInitDialog()
{
	CDialog::OnInitDialog();
	GAFitness();
	UpdateData(FALSE);
	return TRUE;
}

//------------------------***-------------------------------------------*
// Function
//------------------------***-------------------------------------------*
//----------------------------------------------------------------------*
// basis function 
//----------------------------------------------------------------------*
void CdlgFittingProcess::bbasis2(int c, float t, int npts, float *xv, float *nbasis)
{
	float temp0[200];
   	float b1,b2;
	int nplusc;
	int i,k;

	nplusc=npts+c;
// -------------------------------------------------------------------*
//  zero the temporary arrays
// -------------------------------------------------------------------*
/*	for(i=1;i<=nplusc;i++){
	 temp0[i]=0.;
	}
*/
// -------------------------------------------------------------------*
// calculate the first order basis function nbasis(i,1)
// -------------------------------------------------------------------*
	for(i=1;i<=nplusc-1;i++){
	 if(t>=xv[i] && t<xv[i+1] )  temp0[i]=1.0;
	  else temp0[i]=0.0;
	} 

// -------------------------------------------------------------------*
// calculate the higher order basis functions
//--------------------------------------------------------------------*
	    for(k=2;k<=c;k++){
	      for(i=1;i<=nplusc-k;i++){

			if( temp0[i]!=0.0)/* if the lower order basis function is zero skip the calculation */
				b1=(t-xv[i])*temp0[i] / (xv[i+k-1]-xv[i]);
			else
				b1=0.0;

		    if( temp0[i+1]!=0.0)/* if the lower order basis function is zero skip the calculation */
				b2=( xv[i+k]-t )*temp0[i+1] / (xv[i+k]-xv[i+1]);
			else
			    b2=0.0;

			temp0[i]=b1+b2;
	      } //i
	    } //k
	
		if( t==xv[nplusc] )  temp0[npts]=1.0;/*    pick up last point	*/
	    for(i=1;i<=npts;i++) nbasis[i]=temp0[i];
//		if( t==xv[nplusc] )  nbasis[npts]=1.0;
	return ;
}

//----------------------------------------------------------------------*
// Matrix multiplication
//----------------------------------------------------------------------*
int CdlgFittingProcess::getArea()
{
	return area;
}

//----------------------------------------------------------------------*
//function for cubic Uniform Bspline for the curve after fitness
//----------------------------------------------------------------------*
void CdlgFittingProcess::drawFitness(CDC *pDC)
{
	float dt,b0,b1,b2,b3,xp,yp,zp,xd[202],yd[202],zd[202];
	int nisegs,ntsegs,t,i,k,j,l,m,n;
	int scalex = scale;int scaley = scale;int scalez = scale*1.2;

	//----------------------------------------------------------------------*
	//Fitting points
	//----------------------------------------------------------------------*
/**/for ( i = 1; i <= countX; i++) {
		xp = (x[i])*scalex; zp = z[i]*scalez;
		pDC->MoveTo((int) xp-40,(int) zp+40);
		pDC->LineTo((int) xp+40,(int) zp-40);
		pDC->MoveTo((int) xp+40,(int) zp+40);
		pDC->LineTo((int) xp-40,(int) zp-40);
	}
	//----------------------------------------------------------------------*
	//Sections drawing P[] mesh
	//----------------------------------------------------------------------*
	for (k = 1; k<=countPoc; k++){	
		xp = (pocx0[k])*scalex;zp = pocy0[k]*scalez;			
		if (k==1){
			xp = (Dx[1])*scalex;zp = Dz[1]*scalez;
			pDC->MoveTo((int) xp,(int) zp);
		}
		else{
			if (k==countPoc){
				xp = (Dx[countX])*scalex;zp = Dz[countX]*scalez; 
			}
			pDC->LineTo((int) xp,(int) zp);
		}
	}

/**/for (k = 1; k<=countPoc; k++){			//curvature drawing
		xp = (pocx0[k])*scalex;zp = pocy0[k]*scalez;
		pDC->MoveTo((int) xp,(int) zp);
		xp = (pocx1[k])*scalex;zp = pocy1[k]*scalez;
		pDC->LineTo((int) xp,(int) zp);
	}

	for (k = 1; k<=countBest-4; k++){			//draw control points		
		xp = (DrawB2x[k])*scalex;zp = DrawB2z[k]*scalez;			
		pDC->Ellipse(xp-30,zp-40,xp+30,zp+40);
	}
}

//----------------------------------------------------------------------*
//function for cubic Uniform Bspline for original curve
//----------------------------------------------------------------------*
void CdlgFittingProcess::FitnessData()
{
FILE *f1, *f2;
	f1 = fopen("test.db","wb");
	f2 = fopen("FitnessData.txt","w+t");
	
	char title[18];int i,j;
	strcpy(title,"Fitness Generation");

	fwrite(&title,1,18,f1);
	fprintf(f2,"%s\n",title);

//Knots vector in generation
	fprintf(f2,"Knots vector in final generation for curve: ");
	fprintf(f2,"\n");
	for ( i=1;i<=IndividualLoop;i++){
		for ( j=1;j<=space;j++){
			fwrite(&random[i][j],4,1,f1);
			fprintf(f2,"%10.4f",random[i][j]);
		}
		fprintf(f2,"\n");
	}

//RMS min	
	fprintf(f2,"RMS min in Curve: ");fprintf(f2,"\n");
	fwrite(&RMSMin,4,1,f1);
	fprintf(f2,"%10.4f",RMSMin);
	fprintf(f2,"\n");

// Knot vector in Curve above:
	fprintf(f2,"Knot vector in Curve above: ");fprintf(f2,"\n");
	for (j=1;j<=space;j++){
		fwrite(&IndiCrossBest[1][j],4,1,f1);
		fprintf(f2,"%10.4f",IndiCrossBest[1][j]);
	}
	fprintf(f2,"\n");

//Data fitting
/*	fprintf(f2,"Data fitting in each Curve: ");fprintf(f2,"\n");
	for (i=1;i<=space;i++){
		fwrite(&Pxfix[i],4,1,f1);
		fprintf(f2,"%10.4f",Pxfix[i]);
		fwrite(&Pzfix[i],4,1,f1);
		fprintf(f2,"%10.4f",Pzfix[i]);
		if (fmod(i,space) == 0)	fprintf(f2,"\n");
	}
	fprintf(f2,"\n");
*/
// Vertex point
	fprintf(f2,"Vertex point in Curve: ");fprintf(f2,"\n");
	for (i=1;i<=space;i++){
		fwrite(&DrawB2x[i],4,1,f1);
		fprintf(f2,"%10.4f",DrawB2x[i]);
		fwrite(&DrawB2z[i],4,1,f1);
		fprintf(f2,"%10.4f",DrawB2z[i]);
		if (fmod(i,space) == 0)	fprintf(f2,"\n");
	}
	fprintf(f2,"\n");


	fclose(f1);
	fclose(f2);
}
/////////////////////////////////////



void CdlgFittingProcess::MultimatrixF2(float *u, float *v, float *multi, int m,int n)
{
	int r;
	for (int i=1; i<=m; i++){
		multi[i] = 0;
		for (int q=1; q<=n; q++){
			r=q+(i-1)*n;
			multi[i] = multi[i] + u[r] * v[q];
		}
	}
}

void CdlgFittingProcess::GAFitness()
{
	float check,checkb,RMSAverage;
	float aCr,bCr;
	float numS1,numS2,numS3,aC,bC;
	int numS3b,numS4,numS4b,CountMutation,CountMStep3,CountMStep3b;
	int i,j,k,l,q,s,jb,kb;
	float stepT,stepTb;
	float det = 1.0;
	countTemp = 0;

	int cuttingPointBestCount,cuttingPointBestCountb ;
	int beginBestCount,beginBestCountb ;
	int cuttingPointCount,cuttingPointCountb ;
	int beginCount,beginCountb,countS3 ;
	int criterion;
	srand( (unsigned)time( NULL ) );
	area = 3;
	Refresh();
	count= 0;

	float sum = 0;
	int recursive;int n;int m;int temp1;
	int same =0;
	int counttimes = 0;
	begincountX = 0;
	begincountY = 0;
	for (i=0;i<=generation;i++) knotgen[i] = 0;
//----------------------------------------------------------------------*
//Loop for initial population, loop: individuals in initial population
//----------------------------------------------------------------------*
Start:
for (recursive=1;recursive<=IndividualLoop;recursive++){
	//----------------------------------------------------------------------*
	//initial knots for sections & waterlines
	//----------------------------------------------------------------------*
	//----------------------------------------------------------------------*
	//number of data points in each sections 
	//----------------------------------------------------------------------*	
	countKnot[recursive] = 0;	
	for (j=1;j<=99;j++){
		if (random[recursive][j] == 1) break;
		else if (random[recursive][j] < 1 && random[recursive][j] > 0) countKnot[recursive] = countKnot[recursive] + 1;
	}
	for (i = 1;i<=countKnot[recursive]+8;i++){
		tknot[i] = random[recursive][i];
	}
	for (i = countKnot[recursive]+9;i<=99;i++){
		random[recursive][i] = 0;
		tknot[i] = random[recursive][i];
	}
	countKnot[recursive] =countKnot[recursive] + 8 - (4);
	//----------------------------------------------------------------------*
	//nbasis function for sections 
	//----------------------------------------------------------------------*
	check = 0;	
	for (i=begincountX+1; i<=begincountX+countX-1; i++){
		check = check + 
			sqrt((x[i+1] - x[i])*(x[i+1] - x[i]) + (z[i+1] - z[i])*(z[i+1] - z[i]));
	}
	valueT = 0;	n=1;
	for (i=1; i<=countX; i++){
		bbasis2(4,valueT,countKnot[recursive],tknot,nbasis2);
		for (j=1; j<=countKnot[recursive]; j++)	{
			Nmatrix[n] = nbasis2[j];	n=n+1;
		}
		stepT = sqrt((x[i+1] - x[i])*(x[i+1] - x[i]) 
				+    (z[i+1] - z[i])*(z[i+1] - z[i]));
			valueT = valueT + stepT/check;
			if (i==begincountX+countX-1) valueT = 1.0;
	}

	//----------------------------------------------------------------------*
	//Matrix multiplication	
	//----------------------------------------------------------------------*
	for (i=0; i<=(countKnot[recursive])*countX; i++){
		MultiMatrix[i] = 0;
	}

	temp1 = 1;
	for (i=1; i<=(countKnot[recursive])*countX; i++){
		CT[temp1] = Nmatrix[i];
		temp1 = temp1+countX;
		if (fmod(i,countKnot[recursive])==0) temp1=i/countKnot[recursive]+1;
	}	
	
	
	temp1 = 1;k=1;m=1;
	for (i=1; i<=((countKnot[recursive])*(countKnot[recursive])); i++){
		for (j=1; j<=countX; j++){
			q=j+(k-1)*(countX);
			s=j+(m-1)*(countX);
			MultiMatrix[temp1] = MultiMatrix[temp1]+ CT[q] * CT[s];			
		}
		l=(countKnot[recursive]);
		if (fmod(i,l) == 0){
			m=1;
			k=k+1;
		}
		else{ k=k;	m=m+1;}
		temp1 = temp1 + 1;
	}
	//----------------------------------------------------------------------*
	// generate inverse matrix [C-1] 
	//----------------------------------------------------------------------*
	temp1 = temp1 - 1;
	check = 1.0;
	float dett;
	for (i=1; i<=temp1; i++){
		MultiMatrix[i-1] = MultiMatrix[i];
	}
	n=(countKnot[recursive]);
	invmat(MultiMatrix,n,&dett);
	for (i = temp1;i>=0;i--){
		MultiMatrix[i+1] = MultiMatrix[i];
	}
	//----------------------------------------------------------------------*
	//Vertex points for cubic Non-Uniform Bspline
	//----------------------------------------------------------------------*
	for (k=1;k<=(countKnot[recursive])*countX;k++){
		B1[k] = 0;
	}

	temp1 = 1;k=1;m=1;
	for (i=1; i<=countX*(countKnot[recursive]); i++){
		for (j=1; j<=(countKnot[recursive]); j++){	
			q=j+(k-1)*(countKnot[recursive]);
			s=j+(m-1)*(countKnot[recursive]);
			B1[temp1] = B1[temp1] + MultiMatrix[q]*Nmatrix[s];
		}
		if (fmod(i,countX) == 0){
			m=1;
			k=k+1;
		}
		else{ k=k;	m=m+1;}
		temp1 = temp1 + 1;
	}