📄 dlgplotting.cpp

📁 Hi guys, I have a B-spline curve algorithm by using Genetic algorithm. Any interested?
💻 CPP
字号:
// dlgPlotting.cpp : implementation file
//GA

#include "stdafx.h"
#include "GA.h"
#include "dlgPlotting.h"
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include "matrix.h"
#include "knotc.h"
#include <time.h>
#include "dlgFittingProcess.h"

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

#define data 400
#define Draw 10000

/////////////////////////////////////////////////////////////////////////////
// CdlgPlotting dialog
/////////////////////////////////////////////////////////////////////////////
// global variable

	extern int number;
	extern int space;
	extern float x[data],y[data],z[data];
	extern int IndividualLoop, times, area;
	extern int countY,countYd,begincountY;
	extern int countX,countXd,begincountX;
	extern int scale;
	extern int countKnot[102];

	extern float pocx0[data],pocy0[data],pocx1[data],pocy1[data];
	extern int countPoc;	

	float C1[4000],CT[4000],MultiMatrix[4000],Dx[data],Dy[data],Dz[data],B1[4000];
	float Px[data],Py[data],Pz[data],B2x[data],B2y[data],B2z[data];

	 float nbasis2[100],nbasis2b[100];
	 float Nmatrix2[4000];
	 float Nmatrix[4000],Nmatrixb[4000];
	 float C[data][data];
	 float RMS[101],pm[101],tknot[100],tknotb[100];

	float random[101][100],randomb[101][100],IndiCrossBest[2][100],IndiCrossBestb[2][100],IndiCrossBestCurve[50][100];
	int count, stepTCount;
	float valueT,valueTb,Pxfix[data],Pyfix[data],Pzfix[data];
	float RMSMin,RMSMinCurve[100];
	int countBest,countBestb;
	int interiorKnotX,interiorKnotY;
	float DrawPxfix[Draw],DrawPyfix[Draw],DrawPzfix[Draw],DrawB2x[Draw],DrawB2y[Draw],DrawB2z[Draw];	
/////////////////////////////////////////////////////////////////////////////


CdlgPlotting::CdlgPlotting(CWnd* pParent /*=NULL*/)
	: CDialog(CdlgPlotting::IDD, pParent)
{
	//{{AFX_DATA_INIT(CdlgPlotting)
		// NOTE: the ClassWizard will add member initialization here
	//}}AFX_DATA_INIT
}


void CdlgPlotting::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CdlgPlotting)
	DDX_Control(pDX, IDC_MSHFLEXGRID1, m_grid1);
	//}}AFX_DATA_MAP
}


BEGIN_MESSAGE_MAP(CdlgPlotting, CDialog)
	//{{AFX_MSG_MAP(CdlgPlotting)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CdlgPlotting message handlers


BOOL CdlgPlotting::OnInitDialog()
{

	CDialog::OnInitDialog();
	int i,j,k,n;
	char str2[5];
	float check;
	m_grid1.SetRows(101);
	m_grid1.SetCols(8);
	area = 2;
	Refresh();
	srand( (unsigned)time( NULL ) );	

	countX = 0;	begincountX = 0;	
	countY = 0;	begincountY = 0;
	
	n=1;

	countX = space;
//	interiorKnotX = countX+4+1-8;
	interiorKnotX = countX*0.5;				//L =  Lamda*N from paper
	for (i = 1; i <= IndividualLoop; i++) {
		for (j = 5; j <= interiorKnotX+4; j++) {
			random[i][j] = rand();
			random[i][j] = random[i][j]/RAND_MAX;
		}
	}

	for (i = 1; i <= IndividualLoop; i++) {
		for (j = 1; j <= 4; j++) {
			random[i][j] = 0;
		}
				
		for (j = interiorKnotX+5; j <= interiorKnotX+8; j++) {
			random[i][j] = 1;
		}
	}

	for (k = 1; k <= IndividualLoop; k++) {	
		for (i = 5; i <= interiorKnotX+3; i++) {
			for (j = i+1; j <= interiorKnotX+4; j++) {
				if(random[k][i] > random[k][j]){
					check = random[k][i];
					random[k][i] = random[k][j];
					random[k][j] = check;
				}
			}
		}
	}
/*
	for (k = 1; k <= IndividualLoop; k++) {
		for (i = 5; i <= interiorKnotX+4; i++) {
			sprintf(str2,"%7.4f",random[k][i]);
			m_grid1.SetTextMatrix(k,i-3,str2);
		}
	}*/

//	MessageBox("Message","Title",MB_OK);
	GA();
	InitialPopulationData();

	UpdateData(FALSE);
	return TRUE;
}

void CdlgPlotting::Refresh()
{
	char strNum[11];
	m_grid1.Clear();
	m_grid1.SetTextMatrix(0,0,"Individuals");
	for (int i = 1; i < 7; i++) {
		sprintf(strNum,"%d",i);		
		m_grid1.SetTextMatrix(0,i,strNum);
	}

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

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

	}

	for (i = 1; i < IndividualLoop; i++) {
		sprintf(strNum,"%d",1);		
		m_grid1.SetTextMatrix(i,6,strNum);

	}
}

void CdlgPlotting::bbasis(int c, float t, int npts, float xv[], float *nbasis)
{
	
	float temp[200];
   	float b1,b2;
	int nplusc; 
	int i,k;

	float d,e;

	nplusc=npts+c;
// -------------------------------------------------------------------*
//  zero the temporary arrays
// -------------------------------------------------------------------*
/*	for(i=1;i<=nplusc;i++){
	 temp[i-1]=0.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]))
			temp[i] = 1;
	    else
			temp[i] = 0;
	}
// calculate basis fuction
//--------------------------------------------------------------------*
	for (k = 2; k <= c; k++){
    	for (i = 1; i <= nplusc-k; i++){
        	if (temp[i] != 0)    /* if the lower order basis function is zero skip the calculation */
           		d = ((t-xv[i])*temp[i])/(xv[i+k-1]-xv[i]);
	        else
				d = 0;

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

    	    temp[i] = d + e;
		}
	}

	if (t == xv[nplusc]){		/*    pick up last point	*/
 		temp[npts] = 1;
	}

/* put in n array	*/
	for (i = 1; i <= npts; i++) {
    	nbasis[i] = temp[i];
	}

	return ;
}

int CdlgPlotting::getArea()
{
	return area;
}

void CdlgPlotting::draw(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;
	
	//----------------------------------------------------------------------*
	//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/3;zp = pocy1[k]*scalez/3;
		pDC->LineTo((int) xp,(int) zp);
	}
}

void CdlgPlotting::InitialPopulationData()
{
	FILE *f1, *f2;
	f1 = fopen("test.db","wb");
	f2 = fopen("IntialData.txt","w+t");
	
	char title[18];int j;
	strcpy(title,"Initial Generation");

	fwrite(&title,1,18,f1);
	fprintf(f2,"%s\n",title);
//Knots vector in 1st generation
	fprintf(f2,"Knots vector in 1st generation for curve: ");
	fprintf(f2,"\n");
	for (int 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 above: ");fprintf(f2,"\n");
	for (i=1;i<=space;i++){
		fwrite(&DrawB2x[i],8,1,f1);
		fprintf(f2,"%10.4lf",DrawB2x[i]);
		fwrite(&DrawB2z[i],8,1,f1);
		fprintf(f2,"%10.4lf",DrawB2z[i]);
		if (fmod(i,space) == 0)	fprintf(f2,"\n");
	}
	fprintf(f2,"\n");

	fclose(f1);
	fclose(f2);
}

void CdlgPlotting::GA()
{
	float RMSAverage,check,checkb;
	int i,j,k,m,n,l,beginnumber,temp1,q,s;
	float stepT,stepTb;
    float u,w;

	RMSMin = 10000;
	countBest = 0;
	RMSAverage = 0;
	int counttimes;
	int counttest = 0;
//----------------------------------------------------------------------*
//Loop for initial populations (loop: individuals in initial population)
//----------------------------------------------------------------------*
for (int 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++){
		bbasis(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;
	}	

	//----------------------------------------------------------------------*
	//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;
	}
	for (i=1;i<=countX;i++){
		Dx[i] = x[i];
		Dz[i] = z[i];
	}
	
	m=(countKnot[recursive]);
	n=countX;
	MultimatrixF(B1,Dx,B2x,m,n);
	MultimatrixF(B1,Dz,B2z,m,n);
	for (i=1; i<=countX; i++){
		Px[i] = 0;			Py[i] = 0;			Pz[i] = 0;
	}
	//----------------------------------------------------------------------*
	//Calculate a Cartesian product B-spline surface P[]
	//----------------------------------------------------------------------*
	temp1 = 1;
	for (i=1; i<=countX; i++){
		for (j=1; j<=countX; j++){	
			m=j+countX*(i-1);
			Px[temp1] = Px[temp1] + Nmatrix[m]*B2x[j];
			Pz[temp1] = Pz[temp1] + Nmatrix[m]*B2z[j];
		}
		temp1 = temp1 + 1;
	}
	//----------------------------------------------------------------------*
	//Fitness function BIC: Bazes Information Criterion/
	//----------------------------------------------------------------------*
	check = 0;
	for (i = 1;i<=countX;i++){
		k = i;
		check=check +((Px[k] - Dx[k]) * (Px[k] - Dx[k]) 
					+ (Pz[k] - Dz[k]) * (Pz[k] - Dz[k]));
	}
	RMS[recursive] = countX*log(check)+log(countX)*(2*(countX+4-8)+4);
	
	if (RMS[recursive] < 0) RMS[recursive] = RMS[recursive] * (-1);
	if(RMS[recursive] < RMSMin){
		for (i=1; i<=countX; i++){
			Pxfix[i] = Px[i];
			Pzfix[i] = Pz[i];
		}
		countBest = 0;
		for (j=1; j<=(countKnot[recursive]-8+4)+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];
	}
}
	CreateOpPocupine();
}
void CdlgPlotting::MultimatrixF(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];
		}
	}
}
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -