Projektionsmatrix und Modelbetrachtungsmatrix in OpenGL

[megaweber]

MisterX schrieb:

...Und bitte nicht nach Spiele-/Grafikprogrammierung verschieben, da in diesem Forum (Mathematik) hier mehr Mathematik Kompetenz vorliegt.

Und du nicht registriert bist

Noch ne kleine zusätzliche Erklärung:
In homogenen Koordinaten kann man die Translationen und Rotationen geschickt in eine Matrix packen. In 3x3 wären das zwei verschiedene Operationen.

[MisterX]

Danke, jetzt ist mir einiges klarer. Aber eine Frage habe ich noch, denn ich versuche folgenden code zu verstehen:

class TFrustum{

	float projM[16];  //Projektionsmatrix
	float modM[16];   //Modelbetrachtungsmatrix
	float clip[16];
	float Frustum[6][4];
public:
	void NormalizePlane(int pPlane){

		float laenge;

		laenge = sqrt(Frustum[pPlane][0]*Frustum[pPlane][0]+Frustum[pPlane][1]*Frustum[pPlane][1]+Frustum[pPlane][2]*Frustum[pPlane][2]);
		laenge=laenge/1.0;
		Frustum[pPlane][0] = Frustum[pPlane][0]*laenge;
		Frustum[pPlane][1] = Frustum[pPlane][1]*laenge;
		Frustum[pPlane][2] = Frustum[pPlane][2]*laenge;
		Frustum[pPlane][3] = Frustum[pPlane][3]*laenge;
	}

	void berechne_frustum() {

		glGetFloatv(GL_PROJECTION_MATRIX, projM);
		glGetFloatv(GL_MODELVIEW_MATRIX, modM);

		clip[ 0] = modM[ 0]*projM[ 0] + modM[ 1]*projM[ 4] + modM[ 2]*projM[ 8] + modM[ 3]*projM[12];
		clip[ 1] = modM[ 0]*projM[ 1] + modM[ 1]*projM[ 5] + modM[ 2]*projM[ 9] + modM[ 3]*projM[13];
		clip[ 2] = modM[ 0]*projM[ 2] + modM[ 1]*projM[ 6] + modM[ 2]*projM[10] + modM[ 3]*projM[14];
		clip[ 3] = modM[ 0]*projM[ 3] + modM[ 1]*projM[ 7] + modM[ 2]*projM[11] + modM[ 3]*projM[15];
		clip[ 4] = modM[ 4]*projM[ 0] + modM[ 5]*projM[ 4] + modM[ 6]*projM[ 8] + modM[ 7]*projM[12];
		clip[ 5] = modM[ 4]*projM[ 1] + modM[ 5]*projM[ 5] + modM[ 6]*projM[ 9] + modM[ 7]*projM[13];
		clip[ 6] = modM[ 4]*projM[ 2] + modM[ 5]*projM[ 6] + modM[ 6]*projM[10] + modM[ 7]*projM[14];
		clip[ 7] = modM[ 4]*projM[ 3] + modM[ 5]*projM[ 7] + modM[ 6]*projM[11] + modM[ 7]*projM[15];
		clip[ 8] = modM[ 8]*projM[ 0] + modM[ 9]*projM[ 4] + modM[10]*projM[ 8] + modM[11]*projM[12];
		clip[ 9] = modM[ 8]*projM[ 1] + modM[ 9]*projM[ 5] + modM[10]*projM[ 9] + modM[11]*projM[13];
		clip[10] = modM[ 8]*projM[ 2] + modM[ 9]*projM[ 6] + modM[10]*projM[10] + modM[11]*projM[14];
		clip[11] = modM[ 8]*projM[ 3] + modM[ 9]*projM[ 7] + modM[10]*projM[11] + modM[11]*projM[15];
		clip[12] = modM[12]*projM[ 0] + modM[13]*projM[ 4] + modM[14]*projM[ 8] + modM[15]*projM[12];
		clip[13] = modM[12]*projM[ 1] + modM[13]*projM[ 5] + modM[14]*projM[ 9] + modM[15]*projM[13];
		clip[14] = modM[12]*projM[ 2] + modM[13]*projM[ 6] + modM[14]*projM[10] + modM[15]*projM[14];
		clip[15] = modM[12]*projM[ 3] + modM[13]*projM[ 7] + modM[14]*projM[11] + modM[15]*projM[15];
		/

		//rechts
		Frustum[0][0] = clip[ 3] - clip[ 0];  
		Frustum[0][1] = clip[ 7] - clip[ 4];  
		Frustum[0][2] = clip[11] - clip[ 8];  
		Frustum[0][3] = clip[15] - clip[12];  
		NormalizePlane( 0);

		//links
		Frustum[1][0] = clip[ 3] + clip[ 0];
		Frustum[1][1] = clip[ 7] + clip[ 4];
		Frustum[1][2] = clip[11] + clip[ 8];
		Frustum[1][3] = clip[15] + clip[12];
		NormalizePlane( 1);

		//unten
		Frustum[2][0] = clip[ 3] + clip[ 1];
		Frustum[2][1] = clip[ 7] + clip[ 5];
		Frustum[2][2] = clip[11] + clip[ 9];
		Frustum[2][3] = clip[15] + clip[13];
		NormalizePlane( 2);

		//oben
		Frustum[3][0] = clip[ 3] - clip[ 1];
		Frustum[3][1] = clip[ 7] - clip[ 5];
		Frustum[3][2] = clip[11] - clip[ 9];
		Frustum[3][3] = clip[15] - clip[13];
		NormalizePlane( 3);

		//hinten
		Frustum[4][0] = clip[ 3] - clip[ 2];
		Frustum[4][1] = clip[ 7] - clip[ 6];
		Frustum[4][2] = clip[11] - clip[10];
		Frustum[4][3] = clip[15] - clip[14];
		NormalizePlane( 4);

		//vorne
		Frustum[5][0] = clip[ 3] + clip[ 2];
		Frustum[5][1] = clip[ 7] + clip[ 6];
		Frustum[5][2] = clip[11] + clip[10];
		Frustum[5][3] = clip[15] + clip[14];
		NormalizePlane( 5);
	}

	bool Punkt_im_Frustum(float pX,float pY,float pZ){

		for (int i= 0; i < 6;i++) 
		{
			if (Frustum[i][0]*pX + Frustum[i][1]*pY + Frustum[i][2]*pZ + Frustum[i][3] <= 0)
			{

				return false;
			}

		}
		return true;
	}

};

Ich habe es folgendermaßen verstanden.

Es werden die einzelnen Werte der Projektionsmatrix und der Modelbetrachtungsmatrix geladen. Dann wird Modelbetrachtungsmatrix * Projektionsmatrix mittels Matrixmultiplikation berechnet.
Dann werden die Ebenen in Normalvektorform angegeben. Wobei:

Frustum[a][0] ist x der Ebenennormale
Frustum[a][1] ist y der Ebenennormale
Frustum[a][2] ist z der Ebenennormale
Frustum[a][3] gibt die Entfernung zum Durchstoßpunkt von Ebene und normale an.

Habe ich das bis hierhin richtig verstanden?

Ich verstehe den Sinn der folgenden Zeile nicht:

if (Frustum[i][0]*pX + Frustum[i][1]*pY + Frustum[i][2]*pZ + Frustum[i][3] <= 0)

Dies scheint zu berechnen, ob ein Punkt vor oder hinter der Ebene liegt, aber welche Formel / Idee steckt dahinter?
(Sieht etwas aus wie ein Skalarprodukt, aber nicht ganz)

[megaweber]

MisterX schrieb:

Ich verstehe den Sinn der folgenden Zeile nicht:

if (Frustum[i][0]*pX + Frustum[i][1]*pY + Frustum[i][2]*pZ + Frustum[i][3] <= 0)

Dies scheint zu berechnen, ob ein Punkt vor oder hinter der Ebene liegt, aber welche Formel / Idee steckt dahinter?
(Sieht etwas aus wie ein Skalarprodukt, aber nicht ganz)

Schau dir mal Normalenform, Hessesche Normalform an...

Mit einem Vorzeichentest kann man bestimmen auf welcher seite(in Abhängigkeit vom Normalenvektor) der Punkt sich befindet

[MisterX]

Hallo nochmal,

also bei Wikipedia steht zur Berechnung des Abstands von einer Ebene in Hessesche Normalform:

Allgemein erhält man den Abstand s eines beliebigen Punktes P von der Ebene E, indem man den Ortsvektor p von P für in die hessesche Normalform einsetzt:

s= p * n0 - d

Ist s < 0, so liegt P in demselben Halbraum von E wie der Ursprung, bei positivem Vorzeichen von s hingegen im anderen Halbraum.

Aber in dem Programmcode wird in dr Zeile

Frustum[i][0]*pX + Frustum[i][1]*pY + Frustum[i][2]*pZ + Frustum[i][3] <= 0

Frustum[i][3] addiert. Oder ist Frustum[i][3] nicht der Abstand der Ebene vom Ursprung?

[megaweber]

genau!

Es handelt sich um den Abstand der Ebene zum Ursprung. Den musst du natürlich noch draufaddieren. (Das ist dein d in der HNF)

[MisterX]

Das ganze funktioniert für Punkte ganz gut.

Nun möchte ich einen Test für beliebige Rechtecke.

Zuerst dachte ich: Teste einfach, ob die 4 Eckpunkte in Frustum sind.
Aber das 4 Eck kann ja sichtbar sein, ohne das man ein Eckpunkt sieht.

Hat jemand eine Idee wie man so einen Test machen könnte?
(Eine direkte Berechnung, denn mit Aufteilen des Rechtecks in viele kleine habe ich das schon hinbekommen.Dies ist nur entweder zu ungenau oder bei feiner Aufteilung zu langsamm)

[megaweber]

Stichwort:

Cohen-Sutherland Algorithmus!

Bei Fragen: Frag !

[rapso]

MisterX schrieb:

Das ganze funktioniert für Punkte ganz gut.

Nun möchte ich einen Test für beliebige Rechtecke.

Zuerst dachte ich: Teste einfach, ob die 4 Eckpunkte in Frustum sind.
Aber das 4 Eck kann ja sichtbar sein, ohne das man ein Eckpunkt sieht.

Hat jemand eine Idee wie man so einen Test machen könnte?
(Eine direkte Berechnung, denn mit Aufteilen des Rechtecks in viele kleine habe ich das schon hinbekommen.Dies ist nur entweder zu ungenau oder bei feiner Aufteilung zu langsamm)

mathematisch gesehen solltest du die schnittmenge der beiden koerper bilden... auch wenn das fast schon ein wortspiel ist.

[MisterX]

megaweber schrieb:

Stichwort:

Cohen-Sutherland Algorithmus!

Bei Fragen: Frag !

Hab ich mir angesehen. Aber dieser Algorythmus passt nicht ganz, weil ich im Frustum ja Seiten habe die ein 4 Eck bilden, welches kein rechteck ist. Da die 2 gegenüberliegenden Seiten (die nahe seite und die ferne Seite) verschieden lang sind. (z.B bei der rechten Seite)
Außerdem habe ich geehen, dass dieser Algorythmus nur uaf 2D arbeitet.
(Weil immer nur x und y Koordinaten verlang erden)

Also in etwa so:

----------
     \        /
      \      /
       \    /
        ----

[megaweber]

MisterX schrieb:

.... Aber dieser Algorythmus passt nicht ganz, weil ich im Frustum ja Seiten habe die ein 4 Eck bilden, welches kein rechteck ist...
... dass dieser Algorythmus nur uaf 2D arbeitet....

http://de.wikipedia.org/wiki/Cohen-Sutherland-Algorithmus

Dann sei mal ein bisschen kreativ

Benutze die Outcodes (unter Schnitte) um festzulegen ob ein Element ganz, teilweise oder gar nicht innerhalb des Frustums liegt... die 3D Version bekommst du sicher hin

Und als nächsten Schritt musst du alle Geraden(die Seiten von deinem Rechteck) gegen die Frustumebenen clippen, bei denen der Outcode ein Schnitt impliziert...

[MisterX]

Hallo, ich habe mir jetzt eine Methode zusammengebastelt. Diese scheint bis auf Rundungsfehler zu funktionieren. Ist die Linie im Frustum, aber ihre ecken nicht
wird es mal als drinen angezeigt und wenn ich die linie etwas verschiebe mal als drinnen und mals als draußen, schnell wchselnd.

Ich weiß nicht genau wie ich die Rundungsfehler vermeiden kann.

Die Idee ist, dass ich zu den 6 Seiten die den Frustum begrenzen Ebenen berechne. Dann teste ich, ob die Linie eine dieser Ebenen schneidet und berechne den Schnittpunkt. Nun teste ich, ob der Schnittpunkt im Frustum liegt.
(Hierbei scheinen die Rundungsfehler aufzutreten)

class TFrustum{

	Vector3D kollision;

	float projM[16];  //Projektionsmatrix
	float modM[16];   //Modelbetrachtungsmatrix
	float clip[16];
	float Frustum[6][4];
public:
	void NormalizePlane(int pPlane){

		float laenge;

		laenge = sqrt(Frustum[pPlane][0]*Frustum[pPlane][0]+Frustum[pPlane][1]*Frustum[pPlane][1]+Frustum[pPlane][2]*Frustum[pPlane][2]);
		laenge=laenge/1.0;
		Frustum[pPlane][0] = Frustum[pPlane][0]*laenge;
		Frustum[pPlane][1] = Frustum[pPlane][1]*laenge;
		Frustum[pPlane][2] = Frustum[pPlane][2]*laenge;
		Frustum[pPlane][3] = Frustum[pPlane][3]*laenge;
	}

	void berechne_frustum() {
		//Die Projektionsmatrix hat nun diese Form//
		// | 0   1   2   3|
        // | 4   5   6   7|
        // | 8   9  10  11|
		// | 12 13  14  15|
		/////////////////////////////////
		glGetFloatv(GL_PROJECTION_MATRIX, projM);
		//Die Modelbetrachtungsmatrix hat nun diese Form//
		// | 0   1   2   3|
        // | 4   5   6   7|
        // | 8   9  10  11|
		// | 12 13  14  15|
		/////////////////////////////////
		glGetFloatv(GL_MODELVIEW_MATRIX, modM);

		// Normale Matrixmultiplikation Modelbetrachtungsmatrix * Projektionsmatrix
		clip[ 0] = modM[ 0]*projM[ 0] + modM[ 1]*projM[ 4] + modM[ 2]*projM[ 8] + modM[ 3]*projM[12];
		clip[ 1] = modM[ 0]*projM[ 1] + modM[ 1]*projM[ 5] + modM[ 2]*projM[ 9] + modM[ 3]*projM[13];
		clip[ 2] = modM[ 0]*projM[ 2] + modM[ 1]*projM[ 6] + modM[ 2]*projM[10] + modM[ 3]*projM[14];
		clip[ 3] = modM[ 0]*projM[ 3] + modM[ 1]*projM[ 7] + modM[ 2]*projM[11] + modM[ 3]*projM[15];
		clip[ 4] = modM[ 4]*projM[ 0] + modM[ 5]*projM[ 4] + modM[ 6]*projM[ 8] + modM[ 7]*projM[12];
		clip[ 5] = modM[ 4]*projM[ 1] + modM[ 5]*projM[ 5] + modM[ 6]*projM[ 9] + modM[ 7]*projM[13];
		clip[ 6] = modM[ 4]*projM[ 2] + modM[ 5]*projM[ 6] + modM[ 6]*projM[10] + modM[ 7]*projM[14];
		clip[ 7] = modM[ 4]*projM[ 3] + modM[ 5]*projM[ 7] + modM[ 6]*projM[11] + modM[ 7]*projM[15];
		clip[ 8] = modM[ 8]*projM[ 0] + modM[ 9]*projM[ 4] + modM[10]*projM[ 8] + modM[11]*projM[12];
		clip[ 9] = modM[ 8]*projM[ 1] + modM[ 9]*projM[ 5] + modM[10]*projM[ 9] + modM[11]*projM[13];
		clip[10] = modM[ 8]*projM[ 2] + modM[ 9]*projM[ 6] + modM[10]*projM[10] + modM[11]*projM[14];
		clip[11] = modM[ 8]*projM[ 3] + modM[ 9]*projM[ 7] + modM[10]*projM[11] + modM[11]*projM[15];
		clip[12] = modM[12]*projM[ 0] + modM[13]*projM[ 4] + modM[14]*projM[ 8] + modM[15]*projM[12];
		clip[13] = modM[12]*projM[ 1] + modM[13]*projM[ 5] + modM[14]*projM[ 9] + modM[15]*projM[13];
		clip[14] = modM[12]*projM[ 2] + modM[13]*projM[ 6] + modM[14]*projM[10] + modM[15]*projM[14];
		clip[15] = modM[12]*projM[ 3] + modM[13]*projM[ 7] + modM[14]*projM[11] + modM[15]*projM[15];
		//Das Ergebnis ist wieder eine 4*4 Matrix, die clip Matrix
		// | 0   1   2   3|
        // | 4   5   6   7|
        // | 8   9  10  11|
		// | 12 13  14  15|
		/////////////////////////////////

		//Ebenen in "Hessesche Normalform": normale und Abstand zum Ursprung
		//rechts
		Frustum[0][0] = clip[ 3] - clip[ 0];  //x der normalen
		Frustum[0][1] = clip[ 7] - clip[ 4];  //y der normalen
		Frustum[0][2] = clip[11] - clip[ 8];  //z der normalen
		Frustum[0][3] = clip[15] - clip[12];  //Abstand der Ebene zum Ursprung
		NormalizePlane( 0);

		//links
		Frustum[1][0] = clip[ 3] + clip[ 0];
		Frustum[1][1] = clip[ 7] + clip[ 4];
		Frustum[1][2] = clip[11] + clip[ 8];
		Frustum[1][3] = clip[15] + clip[12];
		NormalizePlane( 1);

		//unten
		Frustum[2][0] = clip[ 3] + clip[ 1];
		Frustum[2][1] = clip[ 7] + clip[ 5];
		Frustum[2][2] = clip[11] + clip[ 9];
		Frustum[2][3] = clip[15] + clip[13];
		NormalizePlane( 2);

		//oben
		Frustum[3][0] = clip[ 3] - clip[ 1];
		Frustum[3][1] = clip[ 7] - clip[ 5];
		Frustum[3][2] = clip[11] - clip[ 9];
		Frustum[3][3] = clip[15] - clip[13];
		NormalizePlane( 3);

		//hinten
		Frustum[4][0] = clip[ 3] - clip[ 2];
		Frustum[4][1] = clip[ 7] - clip[ 6];
		Frustum[4][2] = clip[11] - clip[10];
		Frustum[4][3] = clip[15] - clip[14];
		NormalizePlane( 4);

		//vorne
		Frustum[5][0] = clip[ 3] + clip[ 2];
		Frustum[5][1] = clip[ 7] + clip[ 6];
		Frustum[5][2] = clip[11] + clip[10];
		Frustum[5][3] = clip[15] + clip[14];
		NormalizePlane( 5);
	}

	//Liegt ein Punkt im Frustum
	bool Punkt_im_Frustum(float pX,float pY,float pZ) 
	{
		return Punkt_im_Frustum(Vector3D(pX,pY,pZ));
	}

	//Liegt ein Punkt im Frustum
	bool Punkt_im_Frustum(Vector3D v)
	{

		for (int i = 0; i < 6;i++) 
		{
			if (Frustum[i][0]*v.get_x() + Frustum[i][1]*v.get_y() + Frustum[i][2]*v.get_z() + Frustum[i][3] <= 0)
			{

				return false;
			}

		}
		return true;
	}

	//Testen, ob Teile einer Linie gegeben duch die Koordinaten der beiden Endpunkte im Frustum liegen 
	bool Linie_im_Frustum(float p1X,float p1Y,float p1Z,float p2X,float p2Y,float p2Z  ) 
	{
		return Linie_im_Frustum(Vector3D(p1X,p1Y,p1Z),Vector3D(p2X,p2Y,p2Z));
	}

	//Testen, ob Teile einer Linie gegeben duch die Koordinaten der beiden Endpunkte im Frustum liegen 
	bool Linie_im_Frustum(Vector3D v1, Vector3D v2)
	{
		//Liegt ein Eckpunkt der Linie um Frustum, so muß die Linie gemalt werden
		if (Punkt_im_Frustum(v1)  && !Punkt_im_Frustum(v2)) {return true;}
		if (!Punkt_im_Frustum(v1) &&  Punkt_im_Frustum(v2)) {return true;}
		//Liegen beide Eckpunkte im Frustum, so muß die Linie gemalt werden
		if (Punkt_im_Frustum(v1)  &&  Punkt_im_Frustum(v2)) {return true;}

		//Wenn keine Ecke im Frustum liegt, teste, ob ein Schnittpunkt mit irgend einer Seite im Frustum liegt
		for (int i=0;i<6;i++) 
		{
			if ( linie_schneidet_Ebene(v1,v2,Vector3D(Frustum[i][0],Frustum[i][1],Frustum[i][2]),Frustum[i][3]))
			{
				if (Punkt_im_Frustum(kollision)) {return true;}
			}

		}
		//Linie liegt vollkommen außerhalb des Frustums
		return false;
	}

	bool linie_schneidet_Ebene(Vector3D linie_PunktA,Vector3D linie_PunktB,Vector3D normale, float d) 
	{
		//Richtungsvector der Linie berechnen
		Vector3D linie_Richtungsvector=(linie_PunktA-linie_PunktB);

		//Normale und Richtungsvector der Linie sind senkrecht => Ebene und Linie sind parallel => kein Schnitt
		float denominator=skalar_Produkt(normale,linie_Richtungsvector);
		if (abs(denominator)<0.0001)
		{

			return false;
		}

		//Streckenabschnitt des Strahl berechnen
		float s = -((skalar_Produkt(normale,linie_PunktA)+d))/denominator;

		//Liegt der Streckenabschnitt außerhalb von [0,1], so gibt es keinen Schnitt, denn der Schnitt läge außerhalb der Linie
		if(s<0.0f || s> 1.0f) {return false;}

		//Es gibt einen Schnittpunkt
		kollision=linie_PunktA+(s*linie_Richtungsvector);
		return true;
	}

     //Betrag berechnen
     private: float abs(float x) {if (x <0.0f) {return -x;} return x;}

};

Wie kann ich den Rundungsfehler abstellen?

[CStoll]

Einen exakteren Datentyp verwenden (z.B. double anstelle von float) - allerdings wirst du Rundungsfehler nur ganz vermeiden können, wenn du mit unbegrenzt genauen Zahlen rechnest (GMP oder ähnliches).

[MisterX]

Es funktioniert!!!!!!!!
Danke, ihr habt mich auf neue Ideen gebracht!

Es gibt keinen Rundungsfehler. Die Zeile ist falsch

Vector3D linie_Richtungsvector=(linie_PunktA-linie_PunktB);

Ich habe die beiden Punkte der Linie vertauscht(damit zeigt der Richtungsvector in die entgegengesetzte und somit falsche Richtung)

Wenn es Jemanden interessiert: Es muß so aussehen dann funzt es.

Vector3D linie_Richtungsvector=(linie_PunktB-linie_PunktA);