Gauss Algorithmus
-
Hallo zusammen,
mein Problem: Ich wollte den Schnitt von zwei Dreiecken im Raum berechnen. Dazu habe ich die entsprechenden Ebenengleichungen in Parameterdarstellung aufgestellt und anschließend das LGS aufgestellt und versucht mittels Gauss zu lösen.
Leider sind die Ergebnisse nicht wie gewünscht. Angeblich schneiden sich alle Dreiecke.
Hier mal mein Code:
#define NMAX 10 for(it=list_triangles->begin(); it!=list_triangles->end();it++) { for(iit=list_triangles->begin(); iit!=list_triangles->end();iit++) { if(*it != *iit) { double a[NMAX][NMAX]; const float* t1_v1 = (*it)->vertices[0]->floatData(); // A1 const float* t1_v2 = (*it)->vertices[1]->floatData(); // B1 const float* t1_v3 = (*it)->vertices[2]->floatData(); // C1 const float* t2_v1 = (*iit)->vertices[0]->floatData(); // A2 const float* t2_v2 = (*iit)->vertices[1]->floatData(); // B2 const float* t2_v3 = (*iit)->vertices[2]->floatData(); // C2 //lies Matrix // r2 a[0][0] = t2_v2[0] - t2_v1[0]; a[1][0] = t2_v2[1] - t2_v1[1]; a[2][0] = t2_v2[2] - t2_v1[2]; // s2 a[0][1] = t2_v3[0] - t2_v1[0]; a[1][1] = t2_v3[1] - t2_v1[1]; a[2][1] = t2_v3[2] - t2_v1[2]; // r1 a[0][2] = t1_v2[0] - t1_v1[0]; a[1][2] = t1_v2[1] - t1_v1[1]; a[2][2] = t1_v2[2] - t1_v1[2]; // s1 a[0][3] = t1_v3[0] - t1_v1[0]; a[1][3] = t1_v3[1] - t1_v1[1]; a[2][3] = t1_v3[2] - t1_v1[2]; //b a[0][4] = t1_v1[0] - t2_v1[0]; a[1][4] = t1_v1[1] - t2_v1[1]; a[2][4] = t1_v1[2] - t2_v1[2]; double x[NMAX]; if(Gaussalg(a,3, x)) { } } } } // Gauss'sches Eliminationsverfahren mit Zeilenpivotisierung // Argumente: // double a[N][N+1] erweiterte Koeffizientenmatrix Read/Write // int n Anzahl der Gleichungen Read // double x[N] Loesungen Write // Resultat: // int Fehlercode 0 fuer Fehler, 1 fuer Erfolg int Gaussalg (double a[][NMAX], int n, double x[]) { int i, j; // Zeile, Spalte int s; // Elimininationsschritt int pzeile; // Pivotzeile int fehler = 0; // Fehlerflag double f; // Multiplikationsfaktor const double Epsilon = 0.01; // Genauigkeit double Maximum; // Zeilenpivotisierung //extern FILE *fout; s = 0; do { // die einzelnen Eliminationsschritte //fprintf(fout, "Schritt %2i von %2i\n", s+1, n-1); //cout << ("Schritt %1 von %2\n").arg(s+1).arg(n-1) << endl; Maximum = fabs(a[s][s]); // groesstes Element pzeile = s ; // suchen for (i = s+1; i < n; i++) if (fabs(a[i][s]) > Maximum) { Maximum = fabs(a[i][s]) ; pzeile = i; } fehler = (Maximum < Epsilon); if (fehler) break; // nicht loesbar if (pzeile != s) // falls erforderlich, Zeilen tauschen { double h; for (j = s ; j <= n; j++) { h = a[s][j]; a[s][j] = a[pzeile][j]; a[pzeile][j]= h; } } // Elimination --> Nullen in Spalte s ab Zeile s+1 for (i = s + 1; i < n; i++ ) { f = -(a[i][s]/a[s][s]); // Multiplikationsfaktor a[i][s] = 0.0; for (j = s+1; j <= n ; j++) // die einzelnen Spalten a[i][j] += f*a[s][j]; // Addition der Zeilen i, s } s++; } while ( s < n-1 ) ; if (fehler) { cout << "gauss: Gleichungssystem nicht loesbar\n" << endl; //fprintf (fout, "gauss: Gleichungssystem nicht loesbar\n"); return 0; } else { // Berechnen der Loesungen aus der entstandenen Dreiecksmatrix // letzte Zeile x[n-1] = a[n-1][n] / a[n-1][n-1]; // restliche Zeilen for (i = n-2 ; i >= 0; i-- ) { for (j = n-1 ; j > i ; j-- ) { a[i][n] -= x[j]*a[i][j]; // rechte Seite berechnen } x[i] = a[i][n] / a[i][i]; // Loesung } return 1; } }Warum sind alle Schnittests positiv?
-
Ich hab mir das jetzt nicht alles durchgelesen (Quelltext).
Aber so wie ich dich verstehe stellst du Ebenengleichungen auf. Diese beschreiben aber Ebenen und keine Dreiecke.
Dass sich Ebenen irgendwann schneiden, ist wohl klar (außer die sind parallel).Vllt. solltest du prüfen, ob die entstehende Schnittgerade durch die beiden am Schnitt beteiligten Dreiecke läuft, und erst dann ein 1 zurück geben.
Grüße
Franz
-
Hallo,
ja ich dachte mir ich nehm die 3 Punkte der Dreiecke und stellen jeweils Ebenengleichungen auf um dann den Schnitt berechnen zu können.
Ich hab mittlerweile die Fälle isoliert in denen die Dreiecke die gleichen Punkte oder Kanten haben können. Dennoch zeigt er mir immer noch Dreiecke an, die eigentlich keinen Schnitt haben dürften.
-
Ich hab dieses Problem auch mal lösen sollen... Evtl bringt dich der Quelltext ja ein Stückchen weiter - mir ist damals zumindest kein besserer Weg eingefallen - ich guck mal, ob ich noch paar Kommentare hinterlasse ^^
template <typename Tpt> bool TriangleHits(const my::TTriangle<Tpt> &dreieck1, const my::TTriangle<Tpt> &dreieck2) //vorrausgesetzt, die beiden Dreiecke sind beides auch wirklich dreiecke und nicht nur eine Gerade - weiß zwar net mehr genau, wieso das so wichtig war, aber irgendwas war da ^^ { typedef typename Tpt::ValueType valtype; typedef Tpt point; const size_t dimension = point::dimension; typedef my::TVector <dimension, valtype> vector; typedef my::TLine_Segment<point, vector> line_segment; //strecke typedef my::TPlain <vector> plain; plain E1 = my::GetPlainFromPoints <point, vector> ( dreieck1.GetA(), dreieck1.GetB(), dreieck1.GetC() ); //übers Kreuzprodukt ausrechnen und dann exception werfen, wenns 0 ist (also die vectoren linear abhängig sind) line_segment g[3] = { line_segment ( dreieck1.GetA(), vector (dreieck1.GetB()-dreieck1.GetA()) ), line_segment ( dreieck1.GetB(), vector (dreieck1.GetC()-dreieck1.GetB()) ), line_segment ( dreieck1.GetC(), vector (dreieck1.GetA()-dreieck1.GetC()) ) }; plain E2 = my::GetPlainFromPoints <point, vector> ( dreieck2.GetA(), dreieck2.GetB(), dreieck2.GetC() ); line_segment g2[3] = { line_segment ( dreieck2.GetA(), vector (dreieck2.GetB()-dreieck2.GetA()) ), line_segment ( dreieck2.GetB(), vector (dreieck2.GetC()-dreieck2.GetB()) ), line_segment ( dreieck2.GetC(), vector (dreieck2.GetA()-dreieck2.GetC()) ) }; for(size_t i(0); i != 3; ++i) { math::intersection::result res = math::intersects <plain, line_segment, point> (E2, g[i]); if (res == math::intersection::yes || res == math::intersection::contains) //schneiden sich oder sind identisch { for(size_t j(0); j != 3; ++j) { res = math::intersects <plain, line_segment, point> (E1, g2[j]); if (res == math::intersection::yes || res == math::intersection::contains) return true; } } } return false; }Kann dir zwar nicht mehr alles ganz genau ausm Kopf erklären, aber wenn du Fragen hast, könnte ich noch mal überlegen und ggf. versuchen zu erklären ^^
Ums Zusammenzufassen:
Es reicht nicht, zu überprüfen, ob sich die beiden Ebenen schneiden...bb