Test ob eine Punkt einer Ebene in einem Dreieck in dieser Ebene liegt.



  • Hi Leute,

    ich habe einen Test geschrieben, der mathematisch fehlerfrei ist.

    Dieser beruht auf Vorzeichenüberprüfung von drei Vektoren.

    Jedoch habe ich das Problem, dass mir c++ rundungsfehler einbaut und ich statt "0" eine "-0" oder Statt 0,00001 eine -0,0000001 bekomme. Wie kann ich dem Fehler entgegenwirken?

    Gruß

    Alex



  • Das Rechnen mit float und double bewirkt Rundungsfehler, da diese Zahlen periodisch dargestellt werden. Dagegen kannst du in C++ nichts machen außer eine Fehlertoleranz einzubauen.



  • naja die toleranz bringt mir ja nichts...weil ich dann die 0 als -1 bis 1 definieren müsste und da aber auch richtige Dinge wegfallen würden..

    Gibt es keine Möglichkeit irgendeine Exaktheit zu bekommen?



  • Rechne mit Vielfachen in ints.



  • Du kannst mit rationalen Zahlen rechnen. boost.Rational bietet dafür eine fertige Klasse an.



  • uff ich schreibe mal morgen wenn ich am Arbeitsplatz bin den code hin, dann könnt ihr mir helfen mit dem rational da 🙂

    Gruß

    Alex



  • Necrofighter schrieb:

    naja die toleranz bringt mir ja nichts...weil ich dann die 0 als -1 bis 1 definieren müsste und da aber auch richtige Dinge wegfallen würden..

    Gibt es keine Möglichkeit irgendeine Exaktheit zu bekommen?

    Wieso -1 bis 1? Du kannst doch sagen abs(abstand) < 1e-6 liegt für mich in der Ebene. Die Aufagbe klingt nicht so als sei es entscheidend, ob der Punkt exakt oder quasi exakt in der Ebene liegt.



  • Necrofighter schrieb:

    ich habe einen Test geschrieben, der mathematisch fehlerfrei ist.

    Ach ja?

    Necrofighter schrieb:

    Dieser beruht auf Vorzeichenüberprüfung von drei Vektoren.

    Meinst Du vielleicht Skalarprodukte? Wenn ich einen solchen Test schreiben müsste, wären das 3 Skalarprodukte und wenn alle drei positiv sind, ist der Punkt im Dreieck.



  • bool isPointInTriangle(const Point& a, const Point& b, const Point& c, const Point& p)
    
    {
    	Point ab = b - a;
        Point bc = c - b;
        Point ca = a - c;
    
    	Point nab = b - a;
        Point nbc = c - b;
        Point nca = a - c;
    
    	if(isZero(nab)){}
    
    	else
    			{nab.normalize();}
    
    	if(isZero(nbc)){}
    
    	else
    			{nbc.normalize();}
    
    	if(isZero(nca)){}
    
    	else
    			{nca.normalize();}
    
    	Point pa = (p - a) ;
        Point pb = (p - b) ;
        Point pc = (p - c) ;
    
    	Point pab = (p - a) ;// sqNorm(p - a);
        Point pbc = (p - b) ;// sqNorm(p - b);
        Point pca = (p - c) ;// sqNorm(p - c);
    
    	if (isZero(pab)){}
    
    	else 
    		{ pab.normalize(); }
    
    	if (isZero(pbc)){}
    
    	else 
    		{pbc.normalize(); }
    
    	if (isZero(pca)){}
    
    	else 
    		{pca.normalize(); }
    
    	Point cab = crossProduct(nab, pab);
        Point cbc = crossProduct(nbc, pbc);
        Point cca = crossProduct(nca, pca);
    
    	if (isZero(cab))
    
    	{ 
    	   double n = sqNorm(ab);
    	   return  ( n >= sqNorm(pa)  && n >= sqNorm(pb) );
        }
    
    	else if (isZero(cbc))
    
    	{
    	   double n = sqNorm(bc);
    	   return ( n >= sqNorm(pb) && n >= sqNorm(pc) );
        }
    
    	else if (isZero(cca))
    
    	{
    	   double n = sqNorm(ca);
    	   //cout << "n=" << n << endl;
    	   return ( n >= sqNorm(pc) && n >= sqNorm(pa) );
        }
    
    	else
    
    	{
    		return
    			( 
    
    				((cab.getX() <= 0.000001 && cbc.getX() <= 0.000001 && cca.getX() <=0.000001) ||  (cab.getX() >= -0.000001 && cbc.getX() >= -0.000001 && cca.getX() >=-0.000001))
    
    				&&
    
    				((cab.getY() <= 0.000001 && cbc.getY() <= 0.000001 && cca.getY() <=0.000001) ||  (cab.getY() >= -0.000001 && cbc.getY() >= -0.000001 && cca.getY() >=-0.000001))
    
    				&&
    
    				((cab.getZ() <= 0.000001 && cbc.getZ() <= 0.000001 && cca.getZ() <=0.000001) ||  (cab.getZ() >= -0.000001 && cbc.getZ() >= -0.000001 && cca.getZ() >=-0.000001))
    
    				);
    
    	}
    }
    

    Also: Ich übergebe die dreieckspunkte a,b,c und punkt p.

    Ob Punkt p in der Ebene liegt ist nicht entscheidend.

    ich bekomme am Ende wenn alle ifs da ok waren drei normale der Ebene raus. Diese haben die selbe orientierung wenn der Punkt im Dreieck liegt. Ansonsten existiert eben ein Vorzeichenwechsel.

    Das wird am Ende überprüft.

    Gruß

    Alex



  • Du könntest mal die Redundanz deines Codes entfernen. Übrigends gibt es ein Nicht-zeichen:

    if(isZero(nab)){}
    
        else
                {nab.normalize();}
    

    wird zu:

    if(!isZero(nab)){nab.normalize();}
    


  • Welche Redundanz?

    Ja das es ein häßlicher Code ist ist mir klar, ich bin absoluter Neuling was programmieren angeht. Ich verstehe nur nicht, wie ich nun das ganze Lösen könnte.



  • Vielleicht funktioniert diese Version ja... KA ob ichs richtig gemacht habe ist lange her hehe

    struct POINT3D
    {
    	int x;
    	int y;
    	int z;
    };
    
    POINT3D KreuzProdukt(POINT3D a, POINT3D b)
    {
    	POINT3D p;
    
    	p.x = a.y * b.z - a.z * b.y;
    	p.y = a.z * b.x - a.x * b.z;
    	p.z = a.x * b.y - a.y * b.x;
    
    	return p;
    }
    
    bool PunktInEbene(POINT3D a,POINT3D b,POINT3D c, POINT3D p) 
    { 
    	//Normalenform bestimmen!
    	POINT3D e = KreuzProdukt(a, b);
    
    	//Wenn 0 dann liegt in Ebene
    	return !(((p.x - c.x) * e.x) + ((p.y - c.y) * e.y) + ((p.z - c.z) * e.z));
    }
    

    Edit: Ich bin davon ausgegangen, dass a und b Richtungsvektoren sind... Und die Parameterform übergeben wird. Die Richtungsvektoren lassen sich ja einfach bei gegebenen Punkten berechnen.

    Edit2: Oooops Musste mal statt minus sein... Ich habe es jetzt sicherheitshalter auch getestet.



  • meine a,b,c sind Dreickspunkte.

    Es ist ja mir nicht wichtig, ob das ganze in der Ebene liegt, da kann durch Rundungsfehler auch eine kleine Abweichung vorliegen. Es soll aber in einem Dreieck der Ebene liegen und da läufts bei der allerletzten Abfrage aufgrund von rundungsfehlern schief.



  • Oh, naja das könnte dann funktionieren. Du kannst die Fehler/unnötige Sachen sicherlich selbst korrigieren. Ich habe es auch so geändert, das du Punkte angeben kannst.

    #include <iostream>
    #include <cmath>
    
    using namespace std;
    
    #define PI 3.14159265
    
    struct Point3D
    {
    	double x;
    	double y;
    	double z;
    
    	Point3D(double x, double y, double z): x(x), y(y), z(z){}
    	friend bool operator==(Point3D a, Point3D b);
    };
    
    bool operator==(Point3D a, Point3D b)
    {
    	return (a.x == b.x && a.y == b.y && a.z == b.z);
    }
    
    Point3D KreuzProdukt(Point3D a, Point3D b)
    {
    	return Point3D(a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x);
    }
    
    Point3D RichtungsVektor(Point3D anfang, Point3D ende)
    {
    	return Point3D(ende.x - anfang.x, ende.y - anfang.y, ende.z - anfang.z);
    }
    
    bool PunktInEbene(Point3D a, Point3D b, Point3D c, Point3D p)
    {
    	Point3D r1 = RichtungsVektor(c, a);
    	Point3D r2 = RichtungsVektor(c, b);
    
    	Point3D e = KreuzProdukt(r1, r2);
    
    	return !(((p.x - c.x) * e.x) + ((p.y - c.y) * e.y) + ((p.z - c.z) * e.z));
    }
    
    double VektorWinkel(Point3D a, Point3D b)
    {
    	double rad = (a.x * b.x + a.y * b.y + a.z * b.z)/(sqrt(pow(a.x,2)+pow(a.y,2)+pow(a.z,2))*sqrt(pow(b.x,2)+pow(b.y,2)+pow(b.z,2)));
    
    	if(rad < -1.0) rad = -1.0;
    	else if(rad > 1.0) rad = 1.0;
    
    	return 180/PI * acos(rad);
    }
    
    bool PunktInDreieck(Point3D a, Point3D b, Point3D c, Point3D p)
    {
    	if(!PunktInEbene(a,b,c,p))
    		return false;
    
    	Point3D ap = RichtungsVektor(a, p);
    	Point3D bp = RichtungsVektor(b, p);
    	Point3D cp = RichtungsVektor(c, p);
    
    	double winkel = VektorWinkel(ap, bp) + VektorWinkel(ap, cp) + VektorWinkel(bp, cp);
    
    	if((winkel <= 360.000001 && winkel >= 359.999999) || a == p || b == p || c == p)
    		return true;
    
    	return false;
    }
    
    int main () 
    { 
    	Point3D a(0,0,0);
    	Point3D b(2,0,0);
    	Point3D c(0,2,2);
    	Point3D p(1,1,1);
    
    	if(PunktInDreieck(a,b,c,p))
    		cout << "Punkt in Dreieck!" << endl;
    
    	return 0; 
    }
    

    EDIT: Der Vollständigkeit halber...



  • der Test ist mir neu, aber ich werde das mal testen.

    danke



  • Vielleicht gehts auch nur 2D -.-



  • Jedoch müsste meiner ja eigentlich auch normal gehen....



  • Ich habe nochmal die VektorWinkel Funktion überarbeitet und in 3D getestet funktioniert.

    double VektorWinkel(POINT3D a,POINT3D b)
    {
    	double rad = (a.x * b.x + a.y * b.y + a.z * b.z)/(sqrt(pow(a.x,2)+pow(a.y,2)+pow(a.z,2))*sqrt(pow(b.x,2)+pow(b.y,2)+pow(b.z,2)));
    
    	if(rad < -1.0)
    		rad = -1.0;
    	else if( rad > 1.0)
    		rad = 1.0;
    
    	return 180/PI * acos(rad);
    }
    

    Edit: s************** Rundungsfehler.



  • Hallo Necrofighter,

    ich glaube, dass es nur indirekt an Rundungsfehlern liegt. Dein Algorithmus prüft doch letztlich, ob die drei Kreuzproduktvektoren in identische Richtungen zeigen. Zeigen sie zufällig in Richtung einer der Hauptachsen, so reicht schon ein winziges Delta, das P außerhalb der Ebene des Dreiecks liegt, damit der Test fehlschlägt. Und das auch dann, wenn P augenscheinlich mitten im Dreieck liegt.

    Prüfe die Produktvektoren daraufhin ab, dass jeweils das Skalarprodukt eines Paares nicht negativ wird.

    Gruß
    Werner



  • ich verstehe nicht ganz was du meinst 🙂


Anmelden zum Antworten