Kreis aus Kontrollpunkten berechnen

KXII

wenn die punkte nicht auf einer ebene liegen, dann kann man auch keinen kreis mit ALLEN bilden. das sollte doch klar sein.

das ist das selbe als ob man eine ebene mit vier oder mehr punkten bilden möchte und es möglich ist, das ein paar der punkte nicht in der eben liegen.

ein kreis hat auch die eigenschaften einer ebene und deshalb kann mann keinen kreis aus punkten bilden, die nicht auf einer ebene liegen.

was du bilden möchtest würde ich als band bezeichnen. also eine verkettung von mehreren gekrümmten bahnen zu einem durchgehenden objekt.

ich würde dann folgendermaßen vorgehen:

einen aglorithmus besorgen der 3d-bezierkurven berechnet. eine variante des 3d-bezier-algorithmus dürfte wahrscheinlich so arbeiten, das du ihm insgesamt 4 punkte geben musst, und er dann die krümmung zwischen den beiden inneren punkten berrechnet. dann gehst du hin und durchläufst alle deine 8 teilbahnen zwischen den punkten und berechnest ihren verlauf, indem du jeweils die punkte nimmst, die um die teilbahn liegen.

p3
                     *
              p2 *       * p4
Teilbahn b1 ->  /
            p1 *           * p5

              p8 *       * p6
                     *
                     p7

für die teilbahn b1 benutzt du dann p8/p1/p2/p3. der alogithmus berrechnet dann für die teilbahn von p1 nach p2 die bezierkurve. die erbgnisse des algorithmuss musst du dann entsprechend der gewünschten genauigkeit durch die interpolationsschritte anpassen. aus den ergebnissen des algorithmuss hast du den verlauf der bahn in einzelschritte zerlegt und kannst dann weiter machen indem du z.b "16-eckige-kreise" senkrecht zur bahn bildets und die eckpunkte der "kreise" verbindest, dann hast du einen schönen torusähnlichen körper. das ganze in echtzeitveränderbar in 3d dargestellt sieht bestimmt geil aus (so änhlich wie in 3dsmax die grundkörper)

[ Dieser Beitrag wurde am 21.06.2003 um 14:58 Uhr von KXII editiert. ]

Eho

Gibt es auch Kreise die nicht auf einer Ebene liegen???
Wie sind die denn definiert?

KXII

@kreisforscher

das verstehe ich nicht.

wenn drei punkte nicht ausreichen, das muss man doch einen kreis, der nur durch drei punkte beschrieben ist, verändern können ohne die drei punkte ändern zu müssen. wenn ich mir aber einen kreis in 3d vorstelle der 3 punkte hat, dann führt jede veränderung des kreises zu einer veränderung von mindestens einem punkt.

ein kreis im 3d-raum kann entweder gedreht, skaliert oder verschoben werden.

beim verschieben ändern sich die drei punkte.

beim drehen ändern sich die drei punkte auch.

beim skalieren ändern sich die drei punkte auch.

erklär mal bitte wofür man die anderen fünf punkte braucht.

Ein Kreis Kann um Seine eigene Achse Gedreht Sein im 3D-Raum, wobei Dies im 2D-Raum Nicht Möglich ist. Um Dies zu Kompensieren Braucht man Mindestens n+2 Andere Punkte. Im 3D-Raum Ist n=3, Also Braucht man Mindestens 5 Weitere Punkte.

Um Es Klarer zu Formulieren. Im 2D-Raum Ist ein Kreis Bestimmt Durch 3 Punkte. Jedoch Nicht im 3D-Raum. Denn ein Kreis Ist Eigentlich ein 2D-Objekt. Wenn Dieses in den 3D-Raum eintaucht, Muss Die Ebene, in Der es Liegt, Mit Angegeben sein. Eine Ebene Braucht 4 Angaben. Die 8. Angabe Steckt In der Rotation Des Kreises um Seine eigene Achse.

KXII

was ist die eigene achse eines 3d-kreises? die senkrechte zum mittelpunkt des kreises bezogen auf die ebene auf die der kreis liegt? wenn ja, dann ändern sich die drei punkte auch hier, weil das nichts anderes als eine rotation ist.

das wäre dann so, als ob man für eine linie im 3d-raum mehr als 2 punkte benötigt um sie genau zu beschreiben, weil eine linie sich ja auch um ihre eigene achse drehen kann.

ich glaube nicht das das von polo benötigt wird.

[ Dieser Beitrag wurde am 21.06.2003 um 15:46 Uhr von KXII editiert. ]

KXII

wer definiert, das ein kreis ein 2d-objekt ist?

ich meine, wir behandeln das thama hier nicht nach der korrekten SPRACHLICHEN mathematischen definition. wenn das so wäre, dann ist fast alles falsch was wir hier von uns lassen. aber wenn interessiert das schon?

D1BAKEL

Der Kreisforscher ist ein Troll, oder ein Verwirrter mit enormem Selbstvertrauen: 3 Punkte im Raum definieren einen Kreis eindeutig:
Der Kreis im Raum ist definiert als Menge aller Punkte einer Ebene, die von einem bestimmten Punkt dieselbe Entfernung haben. Der Anschauligkeit liegt dieser Punkt auch auf dieser Ebene, dass es auch ohne geht, sieht man daran, dass auch ein Krei entsteht, wenn man von einer Kugel eine Kappe abschnedet. Die entsprechende Ebene kann man ohne Weiteres aus den drei gegebenen Punkten berechnen. Der Rest geschieht auf genau dieser Ebene und wird konstruiert,wie im 2-dimensionalen... Jeder weiter gegebene Punkt würde zur Mehrdeutigkeit führen, man könnte höchstens testen, ob er auf dem Kreis liegt, den die übrigen drei beschreiben. Weitere Punkte bringen aber keine weiteren Informationen.

PS: Es gibt keine Kreise, die nicht in einer Ebene liegen, jedenfalls heißen sie dann nicht mehr so!

mfG D1B

@D1BAKEL:
Wer Nicht viel Ahnung Hat Sollte Lieber still Sein, Also möchte ich Von Dir hierzu in Dieser Sache nichts Mehr außer Einer Entschuldigung Hören. Nur Weil sich Diese Sache Deinen Horizont Überschreitet musst Du nicht Gleich die Fassung Verlieren. Kopf Hoch! Sieh Es positiv: Nun Hast Du etwas Neues, Was Du Lernen Kannst.

KXII

@kreisforscher

ehrlich gesagt: ich glaube das es hier niemanden gibt der sich für deine 2d-definition des kreises interessiert.

EIN KREIS KANN AUCH EIN 3D-GEBILDE SEIN, FEDDISCH!!!

wenn du dir darunter nichts vorstellen kannst dann bist du der mit dem kleineren horizont!

TGGC

Ach mir fällt wieder ein, warum ich nicht mehr Beiträge von Unregistrierten Usern antworten wollte. Sollte mich in Zukunft wirklich an den Vorsatz halten.

Schön' Schrank noch.

Sgt. Nukem

Original erstellt von KXII:
wenn du dir darunter nichts vorstellen kannst dann bist du der mit dem kleineren horizont!

Wo wir g'rad' dabei sind: Wie viele Dimensionen hat eigentlich so ein Horizont...!??

Original erstellt von KXII:
@kreisforscher
EIN KREIS KANN AUCH EIN 3D-GEBILDE SEIN, FEDDISCH!!!:-D

Gerade Darum Sind ja Auch 8 Punkte Notwendig! Du Hast wohl Meine Ausführungen nicht Gelesen! :o

KXII

ich würde sagen ein horizont hat "(1+1/x)^x für x->inf" dimensionen. kein bischen mehr und kein bischen weniger, alles andere würde zu fehler im raumzeitkontinuum führen!

@kreisforscher

danke für den 100%igen beweis das du es dir nicht vorstellen kannst.

[ Dieser Beitrag wurde am 21.06.2003 um 18:11 Uhr von KXII editiert. ]

kreisforscher hat wirklich recht
ich habs grad mal durchgerechnet
ES PASST!!!

Sgt. Nukem

KXII hat wirklich recht! Ein Horizont hat (1+1/x)^x für x->inf Dimensionen!!
Ich hab's g'rad' mal durchgerechnet...
ES PASST!!!

NEIN!!!

Original erstellt von <geraldo>:
NEIN!!!

DOCH!!!!

Achromat

Alle Punkte wirken auf den optimalen kreis, mit mehr punkten kann dieser noch feiner vorhergesagt werden, der Rechenaufwand ist mit irregulären Matritzen- Fehlern dann recht haarig.

Das ganze basiert auch auf die Vorhersage eines Kreises wenn n aneinandergesetzte Punkte die Abstarhlrichtung vorgeben. Z.b erzeugt ein Sensor
eine Punktewolke mit Ausreißern, es ist nun sehr schwer die Ausreißer auf der nicht genau bestimmten Kreisbahn zu lokalisieren, man kann dann die Streuung der letzten Steigung über statistisch Funktionen finden Sigma-Rauschen.

Oftmals werden auch synthetisch so viele Modelkreise über die Punkte Wolke gelegt bis einer am besten passt..

Alle Lösungen sind immer nicht hinreichend zufriedenstellend.

Erste kunde: Die Kreisbahn

//teilt einen gedachten kreis(Radius) in seqmente und erzeugt eine 
//Position für jedes seqment, kreispunktbahn 
void Circle::Filter(CRect rc,int DotCnt,double Radius)
{
	#define  FC (2.0f*PI/(double)SEGMENTS)
	CPoint ptc(rc.CenterPoint());
	double r(rc.Width()/2);
	for(register short nr=0; nr<SEGMENTS; nr++)//calc dotring
	{
		int x = ptc.x + cos(FC*nr) * r;
		int y = ptc.y + sin(FC*nr) * r;
	}
}

Zweite Kunde: Bresenhalm Kreis !

oid CGdiObj::Circle(CDC *pDC, WORD x0, WORD y0, WORD radius)
{//breshalm :-) 
   short f = 1 - radius;
   WORD ddFx = 0,ddFy = -2 * radius;
   WORD x = 0, y = radius;

   while(x < y) 
   {
     if(f >= 0) 
       y--, ddFy+=2, f+=ddFy;

     x++, ddFx+=2, f+=ddFx+1;

     pDC->SetPixel(x0 + x, y0 + y,RGB(0,255,0));
     pDC->SetPixel(x0 - x, y0 + y,RGB(0,255,0));
     pDC->SetPixel(x0 + x, y0 - y,RGB(0,255,0));
     pDC->SetPixel(x0 - x, y0 - y,RGB(0,255,0));
     pDC->SetPixel(x0 + y, y0 + x,RGB(0,255,0));
     pDC->SetPixel(x0 - y, y0 + x,RGB(0,255,0));
     pDC->SetPixel(x0 + y, y0 - x,RGB(0,255,0));
     pDC->SetPixel(x0 - y, y0 - x,RGB(0,255,0));
   }
}

Dritte Kunde Kreis aus drei punkten

bool Circle::Calc(double ScaleX, double ScaleY,double &PosX, double &PosY, double &Rad)
{
	// Die Unbekannten: X, Y und Radius mit Naeherungswerten
	double X=0,Y=0,R=1;
	unsigned int pointCnt = m_Points.size();
	for(unsigned int iteration=0; iteration < pointCnt; iteration++)
	{
	   GdiMatrix A(pointCnt, 3);
	   GdiMatrix b(pointCnt, 1);

		for(unsigned int i=0; i<pointCnt; i++)
		{
			GdiPoint &p = m_Points[i];

			p.x *= ScaleX;
			p.y *= ScaleY;

			// X
			X+=    DX; double px=GdiPoint(X,Y).Cross(p)-R;
			X-=2.0*DX; double mx=GdiPoint(X,Y).Cross(p)-R;
			X+=    DX;
			// Y
			Y+=    DX; double py=GdiPoint(X,Y).Cross(p)-R;
			Y-=2.0*DX; double my=GdiPoint(X,Y).Cross(p)-R;
			Y+=    DX;

			A.SetElement(i,0,(px-mx)/(2.0*DX));
			A.SetElement(i,1,(py-my)/(2.0*DX));
			A.SetElement(i,2,-1);

			double val = GdiPoint(X,Y).Cross(p)-R;
			b.SetElement(i,0,-val);
		}

		GdiMatrix tmpA = A.GetTransposed()*A;
		GdiMatrix tmpB = A.GetTransposed()*b;

		GdiMatrix dx = tmpA.GetInverted() * tmpB;

		X = X+dx.GetElement(0,0);
		Y = Y+dx.GetElement(1,0);
		R = R+dx.GetElement(2,0); 
	} 

	PosX = X, PosY = Y , Rad = R;

    return true;	
}

Vierte Kunde: Matrix Interpolation mit irregulärer Matrix (Kreis aus drei punkten Zweite Version)

bool Circle::Info(GdiPoint a,GdiPoint b,GdiPoint c,int &PosX, int &PosY, int &Rad)
{
	#define COLS 4
    #define ROWS 3

	double *pM,*pMat = new double[COLS*ROWS];
	if(!(pM=pMat))
	 return false;

	pM[0] = 1,pM[1] = -2*a.x,
	pM[2] = -2*a.y,
	pM[3] = -a.x*a.x-a.y*a.y;
	pM+=COLS;

	pM[0] = 1,pM[1] = -2*b.x,
	pM[2] = -2*b.y,
	pM[3] = -b.x*b.x-b.y*b.y;
	pM+=COLS;

	pM[0] = 1,pM[1] = -2*c.x,
	pM[2] = -2*c.y,
	pM[3] = -c.x*c.x-c.y*c.y;
	pM+=COLS;

	for(register int j = 0; j < ROWS; j++)
	{   // Diagonalenfeld normalisieren
		register double q = pMat[j * COLS + j];

		if(q == 0)//Gewährleisten, daß keine 0 in der Diagonale steht
		  for(register int i = j + 1; i < ROWS; i++)
		   if(pMat[i * COLS + j] != 0)// Suche Reihe mit Feld <> 0 und addiere dazu
		   {
			  for(register int k = 0 ; k < COLS; k++)
				pMat[j * COLS + k] += pMat[i * COLS + k];

			  q = pMat[j * COLS + j];
			  break;
		   }

		if(q != 0)// Diagonalen auf 1 bringen
		  for(register int k = 0; k < COLS; k++)
			pMat[j * COLS + k] = pMat[j * COLS + k] / q;

		for(register int i = 0 ; i < ROWS; i++)
		 if(i != j )// Spalten außerhalb der Diagonalen auf 0 bringen
		 {
			q = pMat[i * COLS + j];

			for(register int k = 0; k < COLS; k++)
			  pMat[i * COLS + k] -=  q * pMat[j * COLS + k];    
		  } 
	}

	bool stat(pMat[0] != 1 || pMat[5] != 1 || pMat[10] != 1);

	if (stat)//irrationale matrix
	{
		PosX = pMat[7];
		PosY = pMat[11];
		Rad = (sqrt(PosX*PosX + PosY*PosY - pMat[3]));
	}

	delete pMat;

	return stat;
}

Fünfte Kunde: Die spline Line
bringt einen auch leicht in die klapse.

#ifndef _SPLINE_H
#define _SPLINE_H

#define DIV_FACTOR 4 //adjust this factor to adjust the curve smoothness
class Curve
{
public:
	float  Ax,Ay;
	float  Bx,By;
	float  Cx,Cy;
	int    Ndiv;

	Curve(float ax,float ay,float bx,float by,float cx,float cy,int ndiv) 
	{
		Ax = ax;
		Ay = ay;
		Bx = bx;
		By = by;
		Cx = cx;
		Cy = cy;
		Ndiv = ndiv;
	}

	Curve(float ax,float ay,float bx,float by,float cx,float cy)
	{
		Ax = ax; 
		Ay = ay;
		Bx = bx; 
		By = by;
		Cx = cx; 
		Cy = cy;
		Ndiv = (int)(max(abs((int)Ax),abs((int)Ay))/DIV_FACTOR);
	}
	Curve() 
	{	
	};

	void PutCurve(float ax,float ay,float bx,float by,float cx,float cy) 
	{
		Ax = ax;
		Ay = ay;
		Bx = bx; 
		By = by;
		Cx = cx; 
		Cy = cy;
		Ndiv = (int)(max(abs((int)Ax),abs((int)Ay))/DIV_FACTOR);
	}

	void draw(HDC hdc,float x,float y) 
	{
		int OrigX,OrigY,NewX,NewY;
		float  t,f,g,h;
		if (Ndiv==0)
			Ndiv=1;

		OrigX = (int)x; 
		OrigY= (int)y;
		for(int i=1; i<=Ndiv ; i++)
		{
			t = 1.0f / (float)Ndiv * (float)i;
			f = t*t*(3.0f-2.0f*t);
			g = t*(t-1.0f)*(t-1.0f);
			h = t*t*(t-1.0f);
			NewX = (int)(x + Ax*f + Bx*g + Cx*h);
			NewY = (int)(y + Ay*f + By*g + Cy*h);
			MoveToEx(hdc, OrigX, OrigY, NULL);
			LineTo(hdc, NewX, NewY);
			OrigX = NewX;  
			OrigY=NewY;
		}
	}
	int GetCount()
	{
		if (Ndiv==0)
			Ndiv=1;
		int PointCount = 1;

		for(int i=1; i<=Ndiv ; i++)
		{
			PointCount++;
		}
		return PointCount;
	}
	void GetCurve(float x,float y, POINT points[], int& PointCount)
	{
		int X,Y;
		float  t,f,g,h;
		if (Ndiv==0)
			Ndiv=1;

		X = (int)x; 
		Y= (int)y;
		points[PointCount].x = X;
		points[PointCount].y = Y;
		PointCount++;

		for(int i=1; i<=Ndiv ; i++)
		{
			t = 1.0f / (float)Ndiv * (float)i;
			f = t*t*(3.0f-2.0f*t);
			g = t*(t-1.0f)*(t-1.0f);
			h = t*t*(t-1.0f);
			X = (int)(x + Ax*f + Bx*g + Cx*h);
			Y = (int)(y + Ay*f + By*g + Cy*h);
			points[PointCount].x = X;
			points[PointCount].y = Y;
			PointCount++;
		}
	}

};

class Spline 
{
public:
	float* Px;
	float* Py;
	float* Ax;
	float* Ay;
	float* Bx;
	float* By;
	float* Cx;
	float* Cy;
	float*  k;
	float*  Mat[3];

	int  NP;

	// constructor
	Spline(POINT pt[], int np)
	{
		NP = np;
		Px = new float[NP];
		Py = new float[NP];
		Ax = new float[NP];
		Ay = new float[NP];
		Bx = new float[NP];
		By = new float[NP];
		Cx = new float[NP];
		Cy = new float[NP];
		k = new float[NP];
		Mat[0] = new float[NP];
		Mat[1] = new float[NP];
		Mat[2] = new float[NP];

		for(int i=0;i<NP ;i++) 
		{
			Px[i] = (float)pt[i].x;  
			Py[i] = (float)pt[i].y;
		}
	}

	Spline(float px[] , float py[] , int np) 
	{
		NP = np;
		Px = new float[NP];
		Py = new float[NP];
		Ax = new float[NP];
		Ay = new float[NP];
		Bx = new float[NP];
		By = new float[NP];
		Cx = new float[NP];
		Cy = new float[NP];
		k = new float[NP];
		Mat[0] = new float[NP];
		Mat[1] = new float[NP];
		Mat[2] = new float[NP];

		for(int i=0;i<NP ;i++) 
		{
			Px[i] = px[i];  
			Py[i] = py[i];
		}
	}

	~Spline()
	{
		delete[] Px;
		delete[] Py;
		delete[] Ax;
		delete[] Ay;
		delete[] Bx;
		delete[] By;
		delete[] Cx;
		delete[] Cy;
		delete[] k;
		delete[] Mat[0];
		delete[] Mat[1];
		delete[] Mat[2];
	}

	void Generate() 
	{
		float AMag , AMagOld;
    	// vector A
		for(int i= 0 ; i<=NP-2 ; i++ ) 
		{
			Ax[i] = Px[i+1] - Px[i];
			Ay[i] = Py[i+1] - Py[i];
		}
		// k
		AMagOld = (float)sqrt(Ax[0]*Ax[0] + Ay[0]*Ay[0]);
		for(register int i=0 ; i<=NP-3 ; i++) 
		{
			AMag = (float)sqrt(Ax[i+1]*Ax[i+1] + Ay[i+1]*Ay[i+1]);
			k[i] = AMagOld / AMag;
			AMagOld = AMag;
		}
		k[NP-2] = 1.0f;

		// Matrix
		for(register int i=1; i<=NP-2;i++) 
		{
			Mat[0][i] = 1.0f;
			Mat[1][i] = 2.0f*k[i-1]*(1.0f + k[i-1]);
			Mat[2][i] = k[i-1]*k[i-1]*k[i];
		}
		Mat[1][0] = 2.0f;
		Mat[2][0] = k[0];
		Mat[0][NP-1] = 1.0f;
		Mat[1][NP-1] = 2.0f*k[NP-2];

		// 
		for(register int i=1; i<=NP-2;i++) 
		{
			Bx[i] = 3.0f*(Ax[i-1] + k[i-1]*k[i-1]*Ax[i]);
			By[i] = 3.0f*(Ay[i-1] + k[i-1]*k[i-1]*Ay[i]);
		}
		Bx[0] = 3.0f*Ax[0];
		By[0] = 3.0f*Ay[0];
		Bx[NP-1] = 3.0f*Ax[NP-2];
		By[NP-1] = 3.0f*Ay[NP-2];

		//
		MatrixSolve(Bx);
		MatrixSolve(By);

		for(register int i=0 ; i<=NP-2 ; i++ ) 
		{
			Cx[i] = k[i]*Bx[i+1];
			Cy[i] = k[i]*By[i+1];
		}
	}

	void MatrixSolve(float B[]) 
	{
		float* Work = new float[NP];
		float* WorkB = new float[NP];
		for(int i=0;i<=NP-1;i++) 
		{
			Work[i] = B[i] / Mat[1][i];
			WorkB[i] = Work[i];
		}

		for(int j=0 ; j<10 ; j++) 
		{ ///  need convergence judge
			Work[0] = (B[0] - Mat[2][0]*WorkB[1])/Mat[1][0];
			for(int i=1; i<NP-1 ; i++ ) 
			{
				Work[i] = (B[i]-Mat[0][i]*WorkB[i-1]-Mat[2][i]*WorkB[i+1])
							/Mat[1][i];
			}
			Work[NP-1] = (B[NP-1] - Mat[0][NP-1]*WorkB[NP-2])/Mat[1][NP-1];

			for(register int i=0 ; i<=NP-1 ; i++ ) 
			{
				WorkB[i] = Work[i];
			}
		}
		for(register int i=0 ; i<=NP-1 ; i++ ) 
		{
			B[i] = Work[i];
		}
		delete[] Work;
		delete[] WorkB;
	}

	void draw(HDC hdc) 
	{
		Curve c;
		for(int i=0; i<NP-1 ; i++) 
		{
			c.PutCurve(Ax[i],Ay[i],Bx[i],By[i],Cx[i],Cy[i]);
			c.draw(hdc,Px[i],Py[i]);
		}

	}
	int GetCurveCount()
	{
		Curve c;
		int count = 0;
		for(int i=0; i<NP-1 ; i++) 
		{
			c.PutCurve(Ax[i],Ay[i],Bx[i],By[i],Cx[i],Cy[i]);
			count += c.GetCount();
		}
		return count;
	}
	void GetCurve(POINT points[], int& PointCount)
	{
		Curve c;
		for(int i=0; i<NP-1 ; i++) 
		{
			c.PutCurve(Ax[i],Ay[i],Bx[i],By[i],Cx[i],Cy[i]);
			c.GetCurve(Px[i],Py[i], points, PointCount);
		}
	}
  //////////// closed cubic spline ////////////////////
	void GenClosed() 
	{
		float AMag , AMagOld , AMag0;
        // vector A
		for(int i= 0 ; i<=NP-2 ; i++ ) 
		{
			Ax[i] = Px[i+1] - Px[i];
			Ay[i] = Py[i+1] - Py[i];
		}
		Ax[NP-1] = Px[0] - Px[NP-1];
		Ay[NP-1] = Py[0] - Py[NP-1];

		// k
		AMag0 = AMagOld = (float)sqrt(Ax[0]*Ax[0] + Ay[0]*Ay[0]);
		for(register int i=0 ; i<=NP-2 ; i++) 
		{
			AMag = (float)sqrt(Ax[i+1]*Ax[i+1] + Ay[i+1]*Ay[i+1]);
			k[i] = AMagOld / AMag;
			AMagOld = AMag;
		}
		k[NP-1]=AMagOld/AMag0; 

		// Matrix
		for(register int i=1; i<=NP-1;i++) 
		{
			Mat[0][i] = 1.0f;
			Mat[1][1] = 2.0f*k[i-1]*(1.0f + k[i-1]);
			Mat[2][i] = k[i-1]*k[i-1]*k[i];
		}
		Mat[0][0] = 1.0f;
		Mat[1][0] = 2.0f*k[NP-1]*(1.0f + k[NP-1]);
		Mat[2][0] = k[NP-1]*k[NP-1]*k[0];

		// 
		for(register int i=1; i<=NP-1;i++) 
		{
			Bx[i] = 3.0f*(Ax[i-1] + k[i-1]*k[i-1]*Ax[i]);
			By[i] = 3.0f*(Ay[i-1] + k[i-1]*k[i-1]*Ay[i]);
		}
		Bx[0] = 3.0f*(Ax[NP-1] + k[NP-1]*k[NP-1]*Ax[0]);
		By[0] = 3.0f*(Ay[NP-1] + k[NP-1]*k[NP-1]*Ay[0]);

		//
		MatrixSolveEX(Bx);
		MatrixSolveEX(By);

		for(register int i=0 ; i<=NP-2 ; i++ ) 
		{
			Cx[i] = k[i]*Bx[i+1];
			Cy[i] = k[i]*By[i+1];
		}
		Cx[NP-1] = k[NP-1]*Bx[0];
		Cy[NP-1] = k[NP-1]*By[0];
	}

  ///// tridiagonal matrix + elements of [0][0], [N-1][N-1] //// 
	void MatrixSolveEX(float B[]) 
	{
		float* Work = new float[NP];
		float* WorkB = new float[NP];

		for(int i=0;i<=NP-1;i++) 
		{
			Work[i] = B[i] / Mat[1][i];
			WorkB[i] = Work[i];
		}

		for(int j=0 ; j<10 ; j++) 
		{  // need judge of convergence
			Work[0] = (B[0]-Mat[0][0]*WorkB[NP-1]-Mat[2][0]*WorkB[1])
					/Mat[1][0];
			for(int i=1; i<NP-1 ; i++ ) 
			{
				Work[i] = (B[i]-Mat[0][i]*WorkB[i-1]-Mat[2][i]*WorkB[i+1])
							/Mat[1][i];
			}
			Work[NP-1] = (B[NP-1]-Mat[0][NP-1]*WorkB[NP-2]-Mat[2][NP-1]*WorkB[0])
							/Mat[1][NP-1];

			for(register int i=0 ; i<=NP-1 ; i++ ) 
			{
				WorkB[i] = Work[i];
			}
		}

		for(register int i=0 ; i<=NP-1 ; i++ ) 
		{
			B[i] = Work[i];
		}

		delete[] Work;
		delete[] WorkB;
	}

	void drawClosed(HDC hdc) 
	{
		Curve c;
		for(int i=0; i<NP ; i++) 
		{
			c.PutCurve(Ax[i],Ay[i],Bx[i],By[i],Cx[i],Cy[i]);
			c.draw(hdc ,Px[i],Py[i]);
		}
	}
};

#endif

Hier noch der benötigte Matrix Stuff:

Header für Matrix

class GdiMatrix
{
public:

public:
	friend class	MatrixHelper ;						// used for operator[][]
	// construction and destruction
					GdiMatrix() ;							// default constructor
					GdiMatrix(const GdiMatrix &other) ;		// copy constructor
					GdiMatrix(int nCols, int nRows) ;		// constructs an empty matrix of this size
					GdiMatrix(int size, bool set_diagonal = true) ;	// creates a square matrix
//					GdiMatrix(VARIANT& var) ;				// from a SAFEARRAY variant
	virtual			~GdiMatrix();							// destructor

	// matrix mathematical operations
	GdiMatrix&		operator=(const GdiMatrix &other) ;
	GdiMatrix			operator+(const GdiMatrix &other) const ;
	GdiMatrix			operator-(const GdiMatrix &other) const ;
	GdiMatrix			operator*(const GdiMatrix &other) const ;
	void			operator+=(const GdiMatrix &other) ;
	void			operator-=(const GdiMatrix &other) ;
	void			operator*=(const GdiMatrix &other) ;
	void			operator*=(double value) ;
	friend GdiMatrix	operator*(const GdiMatrix &other, double value) ;
	bool			operator==(const GdiMatrix &other) const ;
	const MatrixHelper	operator[](int nCol) const ;			// reading version
	MatrixHelper	operator[](int nCol) ;					// writing version

	// element access
	bool			SetElement(int nCol, int nRow, double value) ;
#ifdef _DEBUG
	double			GetElement(int nCol, int nRow) const ;
#else
	inline double	GetElement(int nCol, int nRow) const { return m_pData[nCol + nRow * m_NumColumns] ; } ;
#endif
	inline int		GetNumColumns() const { return m_NumColumns ; } ;
	inline int		GetNumRows() const { return m_NumRows ; } ;
	double			SumColumn(int col) const ;
	double			SumRow(int row) const ;
	double			SumColumnSquared(int col) const ;
	double			SumRowSquared(int row) const ;
	double			GetRowMin(int row) const ;
	double			GetRowMax(int row) const ;
	double			GetColumnMin(int col) const ;
	double			GetColumnMax(int col) const ;
	// matrix transposition
	GdiMatrix			GetTransposed() const ;
	void			Transpose() ;
	// matrix inversion
	GdiMatrix			GetInverted() const ;
	void			Invert() ;
	// covariant (A' * A)
	GdiMatrix			Covariant() const ;
	// normalisation
	GdiMatrix			GetNormalised(double min, double max) const ;
	void			Normalise(double min, double max) ;
	// ranges functions
	void			GetNumericRange(double &min, double &max) const ;
	// matrix concatenation
	GdiMatrix			GetConcatinatedColumns(const GdiMatrix& other) const ;
	void			ConcatinateColumns(const GdiMatrix &other) ;
	GdiMatrix			GetConcatinatedRows(const GdiMatrix& other) const ;
	void			ConcatinateRows(const GdiMatrix &other) ;
	// adds an new row / column to the matrix
	void			AddColumn(const double *pData) ;
	void			AddRow(const double *pData) ;
	// sub matrix extraction, setting
	GdiMatrix			ExtractSubMatrix(int col_start, int row_start, int col_size, int row_size) const ;
	void			SetSubMatrix(int col_start, int row_start, const GdiMatrix &other) ;
	GdiMatrix			ExtractDiagonal() const ;
	// squaring the matrix functions
	GdiMatrix			GetSquareMatrix() const ;
	void			MakeSquare() ;
	// export functions/import
	void			CopyToClipboard() const ;

private:
	// internal variables
	int				m_NumColumns ;			// number of columns in matrix
	int				m_NumRows ;				// number of rows in matrix
	double*			m_pData ;				// pointer to data, may be shared among objects
#ifdef _DEBUG
	// variables used in debug for obejct counting
	static int		m_NextObjectNumber ;
	int				m_ObjectNumber ;
#endif
private:
	// private internal functions
	double*			AllocateMemory(int nCols, int nROws) ;
	// reference counting functions
	void			IncrementReferenceCount() ;	// increments the m_pData reference count
	void			DecrementReferenceCount() ;	// decrements the m_pData reference count
	void			DecrementAndRelease() ;		// decrements the count and releases the memory if required
	int				GetReferenceCount() const ;	// returns the m_pData's reference count
	// helper functions
	CString			GetRowAsText(int row) const ;
	static CString	ReadLine(CFile &file) ;		// reads a \r\n delimited line of text from a file
	static int		GetStringToken(CString source, CString &destination, int start, char ch) ;
};

// this class is used to help the operator[][] on a matrix object work correctly
// it only provides the operator[]
class MatrixHelper
{
public:
	friend class	GdiMatrix ;
protected:
	// protected constructor so only friend class can construct
	// a MatrixHelper object
					MatrixHelper(GdiMatrix* pMatrix, int col) : m_pMatrix(pMatrix), m_pMatrixConst(NULL)
						{
						m_Col = col ;
						} ;
					MatrixHelper(const GdiMatrix* const pMatrix, int col) : m_pMatrixConst(pMatrix), m_pMatrix(NULL)
						{
						m_Col = col ;
						} ;
					MatrixHelper(MatrixHelper& other) : m_pMatrix(other.m_pMatrix), m_pMatrixConst(other.m_pMatrixConst)
						{
						m_Col = other.m_Col ;
						} ;
public:
					~MatrixHelper()
						{
						} ;
	double			operator[](int row) const
						{
						ASSERT(row >= 0) ;									// array bounds error
						ASSERT(row < m_pMatrixConst->m_NumRows) ;				// array bounds error
						return m_pMatrixConst->m_pData[row * m_pMatrixConst->m_NumColumns + m_Col] ;
						} ;
	double&			operator[](int row)
						{
						// first check the reference count on our data object to see whether we need to create a copy
						if (m_pMatrix->GetReferenceCount() > 1)
							{
							// we need to make a copy
							double	*pData = m_pMatrix->m_pData ;			// take a copy of the pointer
							m_pMatrix->DecrementReferenceCount() ;			// decrement the current reference count
							m_pMatrix->m_pData = m_pMatrix->AllocateMemory(m_pMatrix->m_NumColumns, m_pMatrix->m_NumRows) ;
							memcpy(m_pMatrix->m_pData, pData, sizeof(double) * m_pMatrix->m_NumColumns * m_pMatrix->m_NumRows) ;
							m_pMatrix->IncrementReferenceCount() ;			// increment the new data's reference count
							}
						ASSERT(row >= 0) ;									// array bounds error
						ASSERT(row < m_pMatrix->m_NumRows) ;					// array bounds error
						ASSERT(m_pMatrix->m_pData) ;
						return m_pMatrix->m_pData[m_Col + row * m_pMatrix->m_NumColumns] ;
						} ;
private:

	GdiMatrix*		m_pMatrix ;
	const GdiMatrix*	m_pMatrixConst ;
	int				m_Col ;						// column index for operator[][]
} ;

Matrix Code

#include "stdafx.h"
#include "Matrix.h"

GdiMatrix::GdiMatrix()
{
	// default constructor, create a 1 * 1 array
	m_NumColumns = 1 ;
	m_NumRows = 1 ;
	m_pData = NULL ;
	m_pData = AllocateMemory(m_NumColumns, m_NumRows) ;
	IncrementReferenceCount() ;			// count the reference to this memory
}

GdiMatrix::GdiMatrix(const GdiMatrix &other)
{
	// copy constructor
	m_pData = NULL ;
	// use the other objects data pointer
	m_NumColumns = other.m_NumColumns ;
	m_NumRows = other.m_NumRows ;
	m_pData = other.m_pData ;				// copy the pointer
	IncrementReferenceCount() ;				// this thread can get the mutex multiple times without blocking
}

GdiMatrix::GdiMatrix(int nCols, int nRows)
{
	// size constructor
	ASSERT(nCols > 0) ;					// matrix size error
	ASSERT(nRows > 0) ;					// matrix size error
	m_pData = NULL ;
	m_NumColumns = nCols ;
	m_NumRows = nRows ;
	m_pData = AllocateMemory(m_NumColumns, m_NumRows) ;
	IncrementReferenceCount() ;				// count the reference to this memory
}

GdiMatrix::GdiMatrix(int size, bool set_diagonal)
{
	// size constructor
	ASSERT(size > 0) ;						// matrix size error
	m_pData = NULL ;
	m_NumColumns = size ;
	m_NumRows = size ;
	m_pData = AllocateMemory(m_NumColumns, m_NumRows) ;
	IncrementReferenceCount() ;				// count the reference to this memory
	// set the dialognal if required
	if (set_diagonal)
		{
		for (int i = 0 ; i < size ; ++i)
			SetElement(i, i, 1.0) ;
		}
}

GdiMatrix::~GdiMatrix()
{
#ifdef _DEBUG
	TRACE("Destroying GdiMatrix object %1d\n", m_ObjectNumber) ;
#endif
	DecrementAndRelease() ;			// free's m_pData if no other references
}

double* GdiMatrix::AllocateMemory(int nCols, int nRows)
{
	ASSERT(nCols > 0) ;							// size error
	ASSERT(nRows > 0) ;							// size error
	// allocates heap memory for an array
	double	*pData = NULL ;

	pData = new double[nCols * nRows + 1] ;			// all in one allocation (+1 for reference count)
	ASSERT(pData != NULL) ;					// allocation error
	ASSERT(FALSE == IsBadReadPtr(pData, sizeof(double) * (nCols * nRows + 1))) ;
	// empty the memory
	memset(pData, 0, sizeof(double) * (nCols * nRows + 1)) ;		// starts with a 0 reference count
	return pData ;
}

GdiMatrix& GdiMatrix::operator=(const GdiMatrix &other)
{
	if (&other == this)
		return *this ;
	// this does the same job as a copy constructor except we have to de-allocate any
	// memory we may have already allocated
	DecrementAndRelease() ;			// free's m_pData if no other references
	// now copy the other matrix into ourselves
	// use the other objects data pointer
	m_NumColumns = other.m_NumColumns ;
	m_NumRows = other.m_NumRows ;
	m_pData = other.m_pData ;				// copy the pointer
	IncrementReferenceCount() ;				// this thread can get the mutex multiple times without blocking
	// finally return a reference to ourselves
	return *this ;
}

bool GdiMatrix::operator==(const GdiMatrix &other) const
{
	// only return true if the matrices are exactly the same
	if (&other == this)
		return true ;		// comparing to ourselves
	if (m_pData == other.m_pData)
		return true ;		// both pointing to same data, must be same
	if (m_NumColumns != other.m_NumColumns || m_NumRows != other.m_NumRows)
		return false ;		// different dimensions
	if (memcmp(m_pData, other.m_pData, sizeof(double) * m_NumColumns * m_NumRows) == 0)
		return true ;		// buffers are the same
	return false ;			// must be different
}

GdiMatrix GdiMatrix::operator+(const GdiMatrix &other) const
{
	// first check for a valid addition operation
	if (m_NumColumns != other.m_NumColumns)
		throw "Invalid operation" ;
	if (m_NumRows != other.m_NumRows)
		throw "Invalid operation" ;
	// now that we know that the operation is possible
	ASSERT(FALSE == IsBadReadPtr(other.m_pData, sizeof(double) * other.m_NumColumns * other.m_NumRows)) ;
	// construct the object we are going to return
	GdiMatrix		result(*this) ;		// copy ourselves
	// now add in the other matrix
	for (int i = 0 ; i < m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < m_NumRows ; ++j)
			result.SetElement(i, j, result.GetElement(i, j) + other.GetElement(i, j)) ;
		}
	return result ;
}

GdiMatrix GdiMatrix::operator-(const GdiMatrix &other) const
{
	// first check for a valid subtraction operation
	if (m_NumColumns != other.m_NumColumns)
		throw "Invalid operation" ;
	if (m_NumRows != other.m_NumRows)
		throw "Invalid operation" ;
	// now that we know that the operation is possible
	ASSERT(FALSE == IsBadReadPtr(other.m_pData, sizeof(double) * other.m_NumColumns * other.m_NumRows)) ;
	// construct the object we are going to return
	GdiMatrix		result(*this) ;		// copy ourselves
	// now subtract the other matrix
	for (int i = 0 ; i < m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < m_NumRows ; ++j)
			result.SetElement(i, j, result.GetElement(i, j) - other.GetElement(i, j)) ;
		}
	return result ;
}

GdiMatrix GdiMatrix::operator*(const GdiMatrix &other) const
{
	// first check for a valid multiplication operation
	if (m_NumRows != other.m_NumColumns)
		throw "Matrices do not have common size" ;
	// now that we know that the operation is possible
	ASSERT(FALSE == IsBadReadPtr(other.m_pData, sizeof(double) * other.m_NumColumns * other.m_NumRows)) ;
	// construct the object we are going to return
	GdiMatrix		result(m_NumColumns, other.m_NumRows) ;

	// e.g.
	// [A][B][C]   [G][H]     [A*G + B*I + C*K][A*H + B*J + C*L]
	// [D][E][F] * [I][J] =   [D*G + E*I + F*K][D*H + E*J + F*L]
	//             [K][L]
	//
	double	 value ;
	for (int i = 0 ; i < result.m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < result.m_NumRows ; ++j)
			{
			value = 0.0 ;
			for (int k = 0 ; k < m_NumRows ; ++k)
				{
				value += GetElement(i, k) * other.GetElement(k, j) ;
				}
			result.SetElement(i, j, value) ;
			}
		}
	return result ;
}

void GdiMatrix::operator+=(const GdiMatrix &other)
{
	// first check for a valid addition operation
	if (m_NumColumns != other.m_NumColumns)
		throw "Invalid operation" ;
	if (m_NumRows != other.m_NumRows)
		throw "Invalid operation" ;
	// now that we know that the operation is possible
	ASSERT(FALSE == IsBadReadPtr(other.m_pData, sizeof(double) * other.m_NumColumns * other.m_NumRows)) ;
	// now add in the other matrix
	for (int i = 0 ; i < m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < m_NumRows ; ++j)
			SetElement(i, j, GetElement(i, j) + other.GetElement(i, j)) ;
		}
}

void GdiMatrix::operator-=(const GdiMatrix &other)
{
	// first check for a valid subtraction operation
	if (m_NumColumns != other.m_NumColumns)
		throw "Invalid operation" ;
	if (m_NumRows != other.m_NumRows)
		throw "Invalid operation" ;
	// now that we know that the operation is possible
	ASSERT(FALSE == IsBadReadPtr(other.m_pData, sizeof(double) * other.m_NumColumns * other.m_NumRows)) ;
	// now subtract the other matrix
	for (int i = 0 ; i < m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < m_NumRows ; ++j)
			SetElement(i, j, GetElement(i, j) - other.GetElement(i, j)) ;
		}
}

void GdiMatrix::operator*=(const GdiMatrix &other)
{
	// first check for a valid multiplication operation
	if (m_NumRows != other.m_NumColumns)
		throw "Matrices do not have common size" ;

	*this = *this * other ;
}

void GdiMatrix::operator*=(double value)
{
	// just multiply the elements by the value
	for (int i = 0 ; i < m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < m_NumRows ; ++j)
			{
			SetElement(i, j, GetElement(i, j) * value) ;
			}
		}
}

// MatrixHelper is only used for this to simulate a GdiMatrix::operator[][]
const MatrixHelper GdiMatrix::operator[](int nCol) const
{
	ASSERT(nCol >= 0) ;							// array bounds error
	ASSERT(nCol < m_NumColumns) ;				// array bounds error
	// construc the MatrixHelper object to allow operator[][] to work
	MatrixHelper	mh(this, nCol) ;
	return mh ;
}

MatrixHelper GdiMatrix::operator[](int nCol)
{
	ASSERT(nCol >= 0) ;							// array bounds error
	ASSERT(nCol < m_NumColumns) ;				// array bounds error
	// construc the MatrixHelper object to allow operator[][] to work
	MatrixHelper	mh(this, nCol) ;
	return mh ;
}

bool GdiMatrix::SetElement(int nCol, int nRow, double value)
{
	// first check the reference count on our data object to see whether we need to create a copy
	if (GetReferenceCount() > 1)
		{
		// we need to make a copy
		double	*pData = m_pData ;				// take a copy of the pointer
		DecrementReferenceCount() ;				// decrement the current reference count
		m_pData = AllocateMemory(m_NumColumns, m_NumRows) ;
		memcpy(m_pData, pData, sizeof(double) * m_NumColumns * m_NumRows) ;
		IncrementReferenceCount() ;				// increment the new data's reference count
		}
	ASSERT(nCol >= 0) ;							// array bounds error
	ASSERT(nCol < m_NumColumns) ;						// array bounds error
	ASSERT(nRow >= 0) ;							// array bounds error
	ASSERT(nRow < m_NumRows) ;						// array bounds error
	ASSERT(m_pData) ;							// bad pointer error
	m_pData[nCol + nRow * m_NumColumns] = value ;
	return true ;
}

#ifdef _DEBUG
// release version is in-line
double GdiMatrix::GetElement(int nCol, int nRow) const
{
	ASSERT(nCol >= 0) ;							// array bounds error
	ASSERT(nCol < m_NumColumns) ;						// array bounds error
	ASSERT(nRow >= 0) ;							// array bounds error
	ASSERT(nRow < m_NumRows) ;						// array bounds error
	ASSERT(m_pData) ;							// bad pointer error
	return m_pData[nCol + nRow * m_NumColumns] ;
}
#endif

// 
// To avoid big hits when constructing and assigning GdiMatrix objects, multiple GdiMatrix's can reference
// the same m_pData member. Only when a matrix becomes different from the other does a new version of the array
// get created and worked with.
//
void GdiMatrix::IncrementReferenceCount()
{
	// get a pointer to the end of the m_pData object where the reference count resides
	int*	pReference = (int*)&m_pData[m_NumColumns * m_NumRows] ;
	++(*pReference) ;				// increment the reference count
	// done!
}

void GdiMatrix::DecrementReferenceCount()
{
	// get a pointer to the end of the m_pData object where the reference count resides
	int*	pReference = (int*)&m_pData[m_NumColumns * m_NumRows] ;
	--(*pReference) ;				// decrement the reference count
	// done!
}

void GdiMatrix::DecrementAndRelease()
{
	// get a pointer to the end of the m_pData object where the reference count resides
	int*	pReference = (int*)&m_pData[m_NumColumns * m_NumRows] ;
	--(*pReference) ;				// decrement the reference count
	if (*pReference == 0)
		{
		// the memory is no longer needed, release it
		delete []m_pData ;
		m_pData = NULL ;
		}
	// done!
}

int GdiMatrix::GetReferenceCount() const
{
	// get a pointer to the end of the m_pData object where the reference count resides
	int*	pReference = (int*)&m_pData[m_NumColumns * m_NumRows] ;
	return *pReference ;
}

GdiMatrix GdiMatrix::GetTransposed() const
{
	GdiMatrix	transposed(*this) ;		// make a copy of ourselves

	transposed.Transpose() ;
	return transposed ;
}

void GdiMatrix::Transpose()
{
	// first check the reference count on our data object to see whether we need to create a copy
	GdiMatrix	mcopy(*this) ;
	// swap the x/y values
	int	copy = m_NumColumns ;
	m_NumColumns = m_NumRows ;
	m_NumRows = copy ;
	// copy across the transposed data
	for (int i = 0 ; i < m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < m_NumRows ; ++j)
			SetElement(i, j, mcopy.GetElement(j, i)) ;
		}
}

GdiMatrix GdiMatrix::GetInverted() const
{
	// matrix inversion will only work on square matrices
	if (m_NumColumns != m_NumRows)
		throw "GdiMatrix must be square." ;
	// return this matrix inverted
	GdiMatrix	copy(*this) ;
	copy.Invert() ;

	return copy ;
}

void GdiMatrix::Invert()
{
	// matrix inversion will only work on square matrices
	// invert ourselves
	if (m_NumColumns != m_NumRows)
		throw "GdiMatrix must be square." ;
	double	e ;

	for (int k = 0 ; k < m_NumColumns ; ++k)
		{
		e = GetElement(k, k) ;
		SetElement(k, k, 1.0) ;
		if (e == 0.0)
			break;//throw "GdiMatrix inversion error" ;

		for (int j = 0 ; j < m_NumColumns ; ++j)
			SetElement(k, j, GetElement(k, j) / e) ;
		for (int i = 0 ; i < m_NumColumns ; ++i)
			{
			if (i != k)
				{
				e = GetElement(i, k) ;
				SetElement(i, k, 0.0) ;
				for (int j = 0 ; j < m_NumColumns ; ++j)
					SetElement(i, j, GetElement(i, j) - e * GetElement(k, j)) ;
				}
			}
		}
}

// A' * A
GdiMatrix GdiMatrix::Covariant() const
{
	GdiMatrix	result ;
	GdiMatrix trans(GetTransposed()) ;

	result = *this * trans ;
	return result ;
}

GdiMatrix GdiMatrix::ExtractSubMatrix(int col_start, int row_start, int col_size, int row_size) const
{
	ASSERT(col_start >= 0) ;						// bad start index
	ASSERT(row_start >= 0) ;						// bad start index
	ASSERT(col_size > 0) ;						// bad size
	ASSERT(row_size > 0) ;						// bad size
	// make sure the requested sub matrix is in the current matrix
	if (col_start + col_size > m_NumColumns)
		throw "Sub matrix is not contained in source" ;
	if (row_start + row_size > m_NumRows)
		throw "Sub matrix is not contained in source" ;

	GdiMatrix sub(col_size, row_size) ;

	for (int i = 0 ; i < col_size ; ++i)
		{
		for (int j = 0 ; j < row_size ; ++j)
			{
			sub.SetElement(i, j, GetElement(col_start + i, row_start + j)) ;
			}
		}
	return sub ;
}

void GdiMatrix::SetSubMatrix(int col_start, int row_start, const GdiMatrix &other)
{
	ASSERT(col_start >= 0) ;						// bad start index
	ASSERT(row_start >= 0) ;						// bad start index
	ASSERT(col_start + other.m_NumColumns <= m_NumColumns) ;	// bad size
	ASSERT(row_start + other.m_NumRows <= m_NumRows) ;	// bad size
	for (int i = 0 ; i < other.m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < other.m_NumRows ; ++j)
			{
			SetElement(col_start + i, row_start + j, other.GetElement(i, j)) ;
			}
		}
}

GdiMatrix GdiMatrix::ExtractDiagonal() const
{
	if (m_NumColumns != m_NumRows)
		throw "Can only extract diagonal from square matrix" ;
	GdiMatrix	diagonal(m_NumColumns, 1) ;

	for (int i = 0 ; i < m_NumColumns ; ++i)
		diagonal.SetElement(i, 0, GetElement(i, i)) ;
	return diagonal ;
}

GdiMatrix GdiMatrix::GetConcatinatedColumns(const GdiMatrix& other) const
{
	if (m_NumRows != other.m_NumRows)
		throw "Cannot concatenate matrices, not same size" ;
	// copy ourselves and then return the concatenated result
	GdiMatrix		result(*this) ;

	result.ConcatinateColumns(other) ;
	return result ;
}

// concatinate the other matrix to ourselves
void GdiMatrix::ConcatinateColumns(const GdiMatrix &other)
{
	if (m_NumRows != other.m_NumRows)
		throw "Cannot concatenate matrices, not same size" ;
	// create a matrix big enough to hold both
	GdiMatrix		result(m_NumColumns + other.m_NumColumns, m_NumRows) ;

	// now populate it
	for (int i = 0 ; i < m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < m_NumRows ; ++j)
			{
			result.SetElement(i, j, GetElement(i, j)) ;
			}
		}
	// now add the other matrix
	for (int i = 0 ; i < other.m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < m_NumRows ; ++j)
			{
			result.SetElement(i + m_NumColumns, j, other.GetElement(i, j)) ;
			}
		}
	*this = result ;					// assign it to us
}

GdiMatrix GdiMatrix::GetConcatinatedRows(const GdiMatrix& other) const
{
	if (m_NumColumns != other.m_NumColumns)
		throw "Cannot concatenate matrices, not same size" ;
	// copy ourselves and then return the concatenated result
	GdiMatrix		result(*this) ;

	result.ConcatinateRows(other) ;
	return result ;
}

void GdiMatrix::ConcatinateRows(const GdiMatrix &other)
{
	if (m_NumColumns != other.m_NumColumns)
		throw "Cannot concatenate matrices, not same size" ;
	// create a matrix big enough to hold both
	GdiMatrix		result(m_NumColumns, m_NumRows + other.m_NumRows) ;

	// now populate it
	for (int i = 0 ; i < m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < m_NumRows ; ++j)
			{
			result.SetElement(i, j, GetElement(i, j)) ;
			}
		}
	// now add the other matrix
	for (int i = 0 ; i < other.m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < m_NumRows ; ++j)
			{
			result.SetElement(i, j + m_NumRows, other.GetElement(i, j)) ;
			}
		}
	*this = result ;					// assign it to us
}

void GdiMatrix::AddColumn(const double *pData)
{
	ASSERT(FALSE == IsBadReadPtr(pData, sizeof(double) * m_NumRows)) ;

	GdiMatrix	result(m_NumColumns + 1, m_NumRows) ;			// costruct the result

	result.SetSubMatrix(0, 0, *this) ;				// copy ouselves across
	// now add the new row
	for (int i = 0 ; i < m_NumRows ; ++i)
		{
		result.SetElement(m_NumColumns, i, pData[i]) ;
		}
	*this = result ;								// assign result to us
}

void GdiMatrix::AddRow(const double *pData)
{
	ASSERT(FALSE == IsBadReadPtr(pData, sizeof(double) * m_NumColumns)) ;

	GdiMatrix	result(m_NumColumns, m_NumRows + 1) ;			// costruct the result

	result.SetSubMatrix(0, 0, *this) ;				// copy ouselves across
	// now add the new row
	for (int i = 0 ; i < m_NumColumns ; ++i)
		{
		result.SetElement(i, m_NumRows, pData[i]) ;
		}
	*this = result ;								// assign result to us
}

GdiMatrix	operator*(const GdiMatrix &other, double value)
{
	GdiMatrix copy(other) ;

	// just multiply the elements by the value
	for (int i = 0 ; i < copy.m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < copy.m_NumRows ; ++j)
			{
			copy.SetElement(i, j, copy.GetElement(i, j) * value) ;
			}
		}
	return copy ;
}

GdiMatrix GdiMatrix::GetSquareMatrix() const
{
	GdiMatrix	copy(*this) ;

	copy.MakeSquare() ;
	return copy ;
}

void GdiMatrix::MakeSquare()
{
	// make the current matrix square by either stepping in the x or y directions
	// square to the smallest side
	int size = m_NumColumns ;
	if (size > m_NumRows)
		size = m_NumRows ;

	GdiMatrix	work(size, size) ;				// construct result
	double	x_step = m_NumColumns / size ;
	double	y_step = m_NumRows / size ;

	for (int i = 0 ; i < size ; ++i)
		{
		for (int j = 0 ; j < size ; ++j)
			work.SetElement(i, j, GetElement((int)(i * x_step), (int)(j * y_step))) ;
		}
	*this = work ;				// copy the result to ourselves
}

GdiMatrix GdiMatrix::GetNormalised(double min, double max) const
{
	GdiMatrix copy(*this) ;
	copy.Normalise(min, max) ;
	return copy ;
}

void GdiMatrix::Normalise(double min, double max)
{
	// get the lower and upper limit values in the matrix
	// we use the range to normalise
	double	e_min ;
	double	e_max ;

	GetNumericRange(e_min, e_max) ;

	double	range = e_max - e_min ;
	double	r_range = max - min ;			// required range
	double	value ;
	for (int i = 0 ; i < m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < m_NumRows ; ++j)
			{
			value = GetElement(i, j) ;
			value -= e_min ;			// 0 - range
			value /= range ;
			value *= r_range ;
			value += min ;
			SetElement(i, j, value) ;
			}
		}
}

// gets the lowest and highest values in the matrix
void GdiMatrix::GetNumericRange(double &min, double &max) const
{
	double	e_min = GetElement(0, 0) ;
	double	e_max = e_min ;
	double	value ;

	for (int i = 0 ; i < m_NumColumns ; ++i)
		{
		for (int j = 0 ; j < m_NumRows ; ++j)
			{
			value = GetElement(i, j) ;
			if (value < e_min)
				e_min = value ;
			else if (value > e_max)
				e_max = value ;
			}
		}
	min = e_min ;
	max = e_max ;
}

double GdiMatrix::SumColumn(int column) const
{
	ASSERT(column >= 0) ;					// bad column
	ASSERT(column < m_NumColumns) ;				// bad column
	double	sum = 0.0 ;

	for (int i = 0 ; i < m_NumRows ; ++i)
		sum += GetElement(column, i) ;
	return sum ;
}

double GdiMatrix::SumRow(int row) const
{
	ASSERT(row >= 0) ;						// bad row
	ASSERT(row < m_NumRows) ;					// bad row
	double	sum = 0.0 ;

	for (int i = 0 ; i < m_NumColumns ; ++i)
		sum += GetElement(i, row) ;
	return sum ;
}

double GdiMatrix::SumColumnSquared(int column) const
{
	double value = SumColumn(column) ;
	return (value * value) ;
}

double GdiMatrix::SumRowSquared(int row) const
{
	double value = SumRow(row) ;
	return (value * value) ;
}

// returns the minimum value in a row of the matrix
double GdiMatrix::GetRowMin(int row) const
{
	ASSERT(row >= 0) ;
	ASSERT(row < m_NumRows) ;
	double	value = GetElement(0, row) ;
	for (int i = 1 ; i < m_NumColumns ; ++i)
		{
		if (GetElement(i, row) < value)
			value = GetElement(i, row) ;
		}
	return value ;
}

// returns the maximum value in a row of the matrix
double GdiMatrix::GetRowMax(int row) const
{
	ASSERT(row >= 0) ;
	ASSERT(row < m_NumRows) ;
	double	value = GetElement(0, row) ;
	for (int i = 1 ; i < m_NumColumns ; ++i)
		{
		if (GetElement(i, row) > value)
			value = GetElement(i, row) ;
		}
	return value ;
}

// returns the minimum value in a column of the matrix
double GdiMatrix::GetColumnMin(int column) const
{
	ASSERT(column >= 0) ;
	ASSERT(column < m_NumColumns) ;
	double	value = GetElement(column, 0) ;
	for (int i = 1 ; i < m_NumRows ; ++i)
		{
		if (GetElement(column, i) < value)
			value = GetElement(column, i) ;
		}
	return value ;
}

// returns the maximum value in a column of the matrix
double GdiMatrix::GetColumnMax(int column) const
{
	ASSERT(column >= 0) ;
	ASSERT(column < m_NumColumns) ;
	double	value = GetElement(column, 0) ;
	for (int i = 1 ; i < m_NumRows ; ++i)
		{
		if (GetElement(column, i) > value)
			value = GetElement(column, i) ;
		}
	return value ;
}

// returns a row as , separated values
CString GdiMatrix::GetRowAsText(int row) const
{
	ASSERT(row >= 0) ;						// bad row
	ASSERT(row < m_NumRows) ;					// bad row
	CString	text ;
	CString	token ;

	for (int i = 0 ; i < m_NumColumns ; ++i)
		{
		if (i > 0)
			text += "," ;
		token.Format("%e", GetElement(i, row)) ;
		text += token ;
		}
	return text ;
}

// read a line of text from the current file
CString GdiMatrix::ReadLine(CFile& file)
{
	CString	line("") ;
	char	ch[3] ;
	DWORD	hBytesRead	= 0 ;

	while (true)
		{
		hBytesRead = file.Read(&ch[0], 2) ;
		if (hBytesRead == 2)
			file.Seek(-1L, CFile::current) ;		
		if (hBytesRead == 0)
			break ;					// end of file
		if (ch[0] == '\n')
			{
			if (ch[1] == '\r')
				file.Seek(1L, CFile::current) ;
			break ;
			}
		if (ch[0] == '\r')
			{
			ch[0] = '\n' ;
			if (ch[1] == '\n')
				file.Seek(1L, CFile::current) ;
			break ;
			}
		if (ch[0] != '\015')
			{
			// ignore LF characters
			ch[1] = '\0' ;
			line += ch ;
			}
		}
	return line ;
}

int GdiMatrix::GetStringToken(CString source, CString &destination, int start, char ch)
{
	ASSERT(start >= 0) ;

	// at the end of the source string ?
	if (start >= source.GetLength())
		{
		destination = "" ;					// no token available
		return source.GetLength() ;			// return @ end of string
		}

	// skip past any termination characters at the start position
	while (start < source.GetLength())
		{
		if (ch == source[start])
			start++ ;
		else
			break ;
		}
	// find the next occurance of the terminating character
	int pos = source.Find(ch, start) ;			// find termination character
	if (pos < 0)
		{
		// terminator not found, just return the remainder of the string
		destination = source.Right(source.GetLength() - start) ;
		return source.GetLength() ;
		}
	// found terminator, get sub string
	destination = source.Mid(start, pos - start) ;
	return pos ;
}

Fazit.. Du musst aus allen Verfahren das beste für Dich nutzen, keines der
Maßnahmen erreicht ein ahha Erlebnis.

Ich Rechteck-Forscher habe festgestellt, das es keine Kreise gibt. Das ist das eigentlich Problem ^^

Viel Spass
Karsten aus Preußen (www.FlexxVision.de)