dlgfittingprocess.cpp

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

	//----------------------------------------------------------------------*
	//Calculate a Cartesian product B-spline surface P[]
	//----------------------------------------------------------------------*
	temp1 = 1;	k=1;
	for (i=1; i<=(countKnot[recursive]); i++)
	{	
		B2x[i] = 0;
		B2y[i] = 0;
		B2z[i] = 0;
	}
	m=(countKnot[recursive]);
	n=countX;
	MultimatrixF2(B1,Dx,B2x,m,n);
	MultimatrixF2(B1,Dz,B2z,m,n);
	//----------------------------------------------------------------------*
	//Calculate a Cartesian product B-spline surface P[]
	//----------------------------------------------------------------------*
	for (i=1; i<=countX; i++){
		Px[i] = 0;			Py[i] = 0;			Pz[i] = 0;
	}
	
	temp1 = 1;
	for (i=1; i<=countX; i++){
		for (j=1; j<=(countKnot[recursive]); j++){	
			m=j+(i-1)*(countKnot[recursive]);
			Px[temp1] = Px[temp1] + Nmatrix[m]*B2x[j];
			Pz[temp1] = Pz[temp1] + Nmatrix[m]*B2z[j];
		}
		temp1 = temp1 + 1;
	}
	//----------------------------------------------------------------------*
	//Fitness function BIC: Bayes Information Criterion
	//----------------------------------------------------------------------*
	RMS[recursive] = 0;
	check = 0;
	for (i = 1;i<=countX;i++){
		check=check +((Px[i] - Dx[i]) * (Px[i] - Dx[i]) 
					+ (Pz[i] - Dz[i]) * (Pz[i] - Dz[i]));
	}

	RMS[recursive] = countX*log(check)+log(countX)*(2*(countKnot[recursive])+4);
	
	countKnot[recursive]  = countKnot[recursive]+4-8;		//From now: change to Interior Knot	
	if (RMS[recursive] < 0) RMS[recursive] = RMS[recursive] * (-1);
	if(RMS[recursive] < RMSMin){
		for (i=1; i<=countX; i++)	{	
			Pxfix[i] = 0;
			Pzfix[i] = 0;
		} 
		for (i=1; i<=countX; i++)	{	
			Pxfix[i] = Px[i];
			Pzfix[i] = Pz[i];
		} 
		countBest = 0;

		for (j=1; j<=(countKnot[recursive])+4+4; j++) {
			IndiCrossBest[1][j] = random[recursive][j];
			countBest = countBest + 1;
		}
		for (i=1; i<=countX; i++){
			DrawB2x[i] = B2x[i];
			DrawB2z[i] = B2z[i];
		}
		RMSMin = RMS[recursive];
	}
}
//	MessageBox("End Matrix","Title",MB_OK);
if (count >= generation) goto End;
else{
	count = count+1; 
	goto Crossover;
}

//----------------------------------------------------------------------*
//Crossover (assume that the loop of Crossover is CrossOverLoop times)
//----------------------------------------------------------------------*

Crossover:

for (i=1; i<=IndividualLoop; i++){
	criterion = 0;

Again1:
	cuttingPointBestCount = 0;	
	beginBestCount = 0;			
	cuttingPointCount = 0;		
	beginCount = 0;		


	if (criterion >= 5000){					// stopping criterion
		aCr = 0.01;
		bCr = 0.99;
	}
	else {

		aCr = rand(); aCr = aCr/RAND_MAX;
		bCr = rand(); bCr = bCr/RAND_MAX;
		if (aCr > bCr) {
			check = aCr;
			aCr = bCr;
			bCr = check;
		}
	}
	//----------------------------------------------------------------------*
	//checking for random condition satisfied knots or not
	//----------------------------------------------------------------------*
	j=1;	k=1;	jb=1;	kb=1;

Again2:	
	if (IndiCrossBest[1][j]>aCr && IndiCrossBest[1][j]<bCr) goto Again2b;
	else{
		j=j+1;
		if (j > countBest){
			criterion = criterion + 1;
			goto Again1;
		}
		else	goto Again2;		
	}

Again2b:
	
Step1:
	if (random[i][k]>aCr && random[i][k]<bCr) goto Step1b;
	else{
		k=k+1;
		if (k > countKnot[i]+8){
			criterion = criterion + 1;
			goto Again1;
		}
		else	goto Step1;
	}

Step1b:


Cont:										//Crossover processing
	//----------------------------------------------------------------------*
	//counting number in knots satisfied random condition
	//----------------------------------------------------------------------*
//for A	
	for(j=1; j<=countBest; j++){
		if (IndiCrossBest[1][j]>=aCr && IndiCrossBest[1][j]<=bCr){
			cuttingPointBestCount = cuttingPointBestCount + 1;
		}
		else{
			if (IndiCrossBest[1][j] <aCr )	beginBestCount = beginBestCount + 1;
		}
	}
//for A
	for(j=1; j<=countKnot[i]+8; j++){
		if (random[i][j]>=aCr && random[i][j]<=bCr){
			cuttingPointCount = cuttingPointCount + 1;
		}
		else{
			if (random[i][j] <aCr )	beginCount = beginCount + 1;
		}
	}

	//----------------------------------------------------------------------*
	//cutting and replacing
	//----------------------------------------------------------------------*
//for A
	if (cuttingPointBestCount + countKnot[i] - cuttingPointCount + 8 >99){
		criterion = 0;
		goto Again1;
	}
	else{
//for A
		if (cuttingPointBestCount < cuttingPointCount){
			for (j=beginCount+1,k=beginBestCount+1;k<=beginBestCount+cuttingPointBestCount;j++,k++){
				random[i][j] = IndiCrossBest[1][k];
			}
		
			for (j=beginCount+cuttingPointCount+1;j<=countKnot[i]+8+1;j++){
				random[i][j-(cuttingPointCount - cuttingPointBestCount)] = random[i][j];
			}
		}
		else {
			if (cuttingPointBestCount > cuttingPointCount){
				for (j=countKnot[i]+8;j>=beginCount+cuttingPointCount+1;j--){
					random[i][j+(cuttingPointBestCount - cuttingPointCount)] = random[i][j];
				}
				for (j=beginCount+1,k=beginBestCount+1;k<=beginBestCount+cuttingPointBestCount;j++,k++){
					random[i][j] = IndiCrossBest[1][k];
				}
			}

			else {
				for (j=beginCount+1,k=beginBestCount+1;k<=beginBestCount+cuttingPointBestCount;j++,k++){
					random[i][j] = IndiCrossBest[1][k];
				}
			}
		}
	}
	}

	goto Mutation;

//----------------------------------------------------------------------*
//Mutation 
//----------------------------------------------------------------------*
Mutation:

	RMSAverage = 0;
	for (i=1; i<=IndividualLoop;i++){
		RMSAverage = RMSAverage + RMS[i];
	}
	RMSAverage = RMSAverage / IndividualLoop;

	for (i=1; i<=IndividualLoop;i++){
		CountMutation = 0;	
		if (RMS[i] <= RMSAverage){
			pm[i] = 0.5*(RMS[i] - RMSMin) / (RMSAverage - RMSMin);
		}
		else pm[i] = 0.5;

Random_Step1:
	//----------------------------------------------------------------------*
	//re_counting interior knots
	//----------------------------------------------------------------------*
//for A
		countKnot[i] = 0;
		for (j=1;j<=99;j++){
			if (random[i][j] == 1) break;
			else if (random[i][j] < 1 && random[i][j] > 0) countKnot[i] = countKnot[i] + 1;
		}

		numS1 = rand();numS1 = numS1/RAND_MAX;
		if (numS1<=pm[i]){ 
			numS2 = rand();numS2 = numS2 / RAND_MAX;
			if (numS2 <= 0.5){
Random_Step3:				
	//----------------------------------------------------------------------*
	//add numS3 into knot
	//----------------------------------------------------------------------*
				numS3  = rand();numS3  = numS3  / RAND_MAX;
				if (numS3 == 0) goto Random_Step3;

// stopping criterion				
				else{
//for A
					for (j=1;j<=countKnot[i]+8;j++){
						if (numS3  == random[i][j]) goto Random_Step3;
					}
				}

				CountMStep3 = 0;	CountMStep3b = 0;
//for A
				if (countKnot[i]+8<99){
					for (j = 1; j<=countKnot[i]+8;j++){
						if(random[i][j] < numS3)	CountMStep3 = CountMStep3 + 1;
					}

					for (j = countKnot[i]+8;j>=CountMStep3+1;j--){
						random[i][j+1] = random[i][j];
					}
					random[i][CountMStep3 + 1] = numS3;
				}
			}			
	
			else{ 
Random_Step4:
	//----------------------------------------------------------------------*
	//minus numS4 from knots
	//----------------------------------------------------------------------*
//for A			
				if (countKnot[i] > 2){
					numS4 = 5 + fmod(rand(),(countKnot[i]));
					for (j = numS4;j<=countKnot[i]+8;j++){
						random[i][j] = random[i][j+1];
					}
				}
			}
		}

Random_Step5:

		CountMutation = CountMutation + 1;				// stopping criterion
		if (CountMutation < countKnot[i]){
			goto Random_Step1;
		}
	}

	goto Start;

End:
//	MessageBox("Finish","Title",MB_OK);
	CreateOpPocupine();
	FitnessData();
}


//----------------------------------------------------------------------*
//Plan View
//----------------------------------------------------------------------*
void CdlgFittingProcess::drawPlan(CDC *pDC)
{
	float dt,b0,b1,b2,b3,xp,yp,xd[202],yd[202];
	int nisegs,ntsegs,i,k,n,l,timestest;

	int scalex = 180;int scaley = 200;int scalez = 1000;

	for (k = 1; k<=times; k++){
		countY = 0;
		begincountY = 0;
		n=1;
		while ( n <= space) {
			if (y[n+space*(k-1)] == 0){
				if (countY == 0) begincountY = begincountY+1;
			}
			else {countY = countY + 1;}
			n = n+1;
		}
// Draw vertex points
		for ( i = 1; i <= countY; i++) {
			if (y[i+space*(k-1)] != 0){
				xp = x[k+space*(i-1)]*scalex; yp = y[i+space*(k-1)]*scaley;
				pDC->MoveTo((int) xp-40,(int) yp+40);
				pDC->LineTo((int) xp+40,(int) yp-40);
				pDC->MoveTo((int) xp+40,(int) yp+40);
				pDC->LineTo((int) xp-40,(int) yp-40);
			}
		}
	}

	timestest = countY;
	int countX,begincountX;
// Draw curve points
	for (k = 1; k<=timestest; k++){
		countX = 0;
		begincountX = 0;
		n=1;
		while ( n <= space) {
			if (x[n+space*(k-1)] == 0){
				if (countX == 0) begincountX = begincountX+1;
			}
			else {countX = countX + 1;}
			n = n+1;
		}
		countYd = 3;	
		xd[1]=x[begincountX+1+space*(k-1)]*scalex;	yd[1]=DrawPyfix[k]*scaley;
		xd[2]=x[begincountX+1+space*(k-1)]*scalex;	yd[2]=DrawPyfix[k]*scaley;

		xd[countX+3]=x[space*(k-1)+countX+begincountX]*scalex;	yd[countX+3]=DrawPyfix[space*(countX+begincountX-1)+k]*scaley;
		xd[countX+4]=x[space*(k-1)+countX+begincountX]*scalex;	yd[countX+4]=DrawPyfix[space*(countX+begincountX-1)+k]*scaley;
		xd[countX+5]=x[space*(k-1)+countX+begincountX]*scalex;	yd[countX+5]=DrawPyfix[space*(countX+begincountX-1)+k]*scaley;

		for ( l=begincountX+1;l<=begincountX+countX;l++){
			xd[countYd]=x[l+space*(k-1)]*scalex;
			yd[countYd]=DrawPyfix[k+space*(l-1)]*scaley;
			countYd = countYd+1;
		}

		nisegs = countX-3; ntsegs = 100;
		dt = 1 / (float) ntsegs;

		for (i=3;i<=nisegs+2+5;i++){
			for (int j=0,t=0.;j<=ntsegs;j++,t+=dt){
				b0 = t*t*t/6.;
				b1 = (1.+3.*t+3.*t*t-3.*t*t*t)/6.;
				b2 = (4.-6.*t*t+3.*t*t*t)/6.;
				b3 = (1.-t)*(1.-t)*(1.-t)/6.;
				xp = b3*xd[i-2]+b2*xd[i-1]+b1*xd[i]+b0*xd[i+1];
				yp = b3*yd[i-2]+b2*yd[i-1]+b1*yd[i]+b0*yd[i+1];
				if (i==3 && j==0)
					pDC->MoveTo((int) xp,(int) yp);
				else
					pDC->LineTo((int) xp,(int) yp);
			}
		}


	}

}

//----------------------------------------------------------------------*
//Plan Profile
//----------------------------------------------------------------------*
void CdlgFittingProcess::drawSurf(CDC *pDC)
{
	float dt,b0,b1,b2,b3,xp,zp,xd[202],zd[202];
	int nisegs,ntsegs,i,k,n,l,timestest;

	int scalex = 180;int scaley = 200;int scalez = 100;
	cubicBSplSurf();
	plotSurface(pDC);
/*	for (k = 1; k<=times; k++){
		countY = 0;
		begincountY = 0;
		n=1;
		while ( n <= space) {
			if (y[n+space*(k-1)] == 0){
				if (countY == 0) begincountY = begincountY+1;
			}
			else {countY = countY + 1;}
			n = n+1;
		}
// Draw vertex points
		for ( i = 1; i <= countY; i++) {
			if (y[i+space*(k-1)] != 0){
				xp = x[k+space*(i-1)]*scalex; zp = z[i+space*(k-1)]*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);
			}
		}
	}

	timestest = countY;
	int countX,begincountX;
// Draw curve points
	for (k = 1; k<=timestest; k++){
		countX = 0;
		begincountX = 0;
		n=1;
		while ( n <= space) {
			if (x[n+space*(k-1)] == 0){
				if (countX == 0) begincountX = begincountX+1;
			}
			else {countX = countX + 1;}
			n = n+1;
		}
		countYd = 3;	
		xd[1]=x[begincountX+1+space*(k-1)]*scalex;	zd[1]=DrawPzfix[k]*scalez;
		xd[2]=x[begincountX+1+space*(k-1)]*scalex;	zd[2]=DrawPzfix[k]*scalez;

		xd[countX+3]=x[space*(k-1)+countX+begincountX]*scalex;	zd[countX+3]=DrawPzfix[space*(countX+begincountX-1)+k]*scalez;
		xd[countX+4]=x[space*(k-1)+countX+begincountX]*scalex;	zd[countX+4]=DrawPzfix[space*(countX+begincountX-1)+k]*scalez;
		xd[countX+5]=x[space*(k-1)+countX+begincountX]*scalex;	zd[countX+5]=DrawPzfix[space*(countX+begincountX-1)+k]*scalez;

		for ( l=begincountX+1;l<=begincountX+countX;l++){
			xd[countYd]=x[l+space*(k-1)]*scalex;
			zd[countYd]=DrawPzfix[k+space*(l-1)]*scalez;
			countYd = countYd+1;
		}

		nisegs = countX-3; ntsegs = 100;
		dt = 1 / (float) ntsegs;

		for (i=3;i<=nisegs+2+5;i++){
			for (int j=0,t=0.;j<=ntsegs;j++,t+=dt){
				b0 = t*t*t/6.;
				b1 = (1.+3.*t+3.*t*t-3.*t*t*t)/6.;
				b2 = (4.-6.*t*t+3.*t*t*t)/6.;
				b3 = (1.-t)*(1.-t)*(1.-t)/6.;
				xp = b3*xd[i-2]+b2*xd[i-1]+b1*xd[i]+b0*xd[i+1];
				zp = b3*zd[i-2]+b2*zd[i-1]+b1*zd[i]+b0*zd[i+1];
				if (i==3 && j==0)
					pDC->MoveTo((int) xp,(int) zp);
				else
					pDC->LineTo((int) xp,(int) zp);
			}
		}


	}
*/
}