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

krümelkacker

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.

Necrofighter

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

HighLigerBiMBam

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();}

Necrofighter

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.

HighLigerBiMBam

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.

Necrofighter

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.

HighLigerBiMBam

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...

Necrofighter

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

danke

HighLigerBiMBam

Vielleicht gehts auch nur 2D -.-

Necrofighter

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

HighLigerBiMBam

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.

Werner Salomon

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

Necrofighter

ich verstehe nicht ganz was du meinst

Werner Salomon

Necrofighter schrieb:

ich verstehe nicht ganz was du meinst

nehme die Vektoren cab, cbc und cca und bereche von zwei Paaren (z.B. cab*cbc und cbc*cca) das Skalarprodukt. Das Ergebnis darf nie kleiner 0 sein, sonst liegt P außerhalb des Dreiecks.

Necrofighter

Das Problem ist ja, dass die nicht vollkommen identisch sind

Es sollsten ja vom Betrag 3 gleiche Vektoren da stehen, aber durch die Schwankung wird auch dein Test nicht gehen.

Werner Salomon

Necrofighter schrieb:

Das Problem ist ja, dass die nicht vollkommen identisch sind

Es sollsten ja vom Betrag 3 gleiche Vektoren da stehen, aber durch die Schwankung wird auch dein Test nicht gehen.

natürlich sind die nicht identisch; sie müssen nur in eine 'ähnliche' Richtung zeigen - unabhängig von ihrere Länge. Ist das Skalarprodukt zweier belieber Vektoren >0 so ist ihr Winkel untereinander <90°.
Und es geht doch nur darum, ob alle Kreuzproduktvektoren zur selben Seite des Dreiecks hin zeigen, oder eben vom Betrag 0 sind. In letzteren Fall liegt P auf einer der Seiten.

Wenn Du alle drei Paare untersuchst, so kannst Du auch die if( isZero(..) )-Kaskade einfach weglassen.

Gruß
Werner

Necrofighter

Alles klar ich teste das mal.

Danke für die gute Hilfe.

Kannst mir ja vielleicht in meinem anderen Thread mit der map of arrays versuchen zu helfen

Necrofighter

Mal Mathematik:

Wenn ich die Determinante von (ab,ap,111) und (bc,bp,111) und (ca,cp,111)
ausrechne, dann ist beim Vorzeichenwechsel der Punkt ausserhalb und wenn nicht sitzt er drin.

Als weitere Lösung für das Problem

geom

Hallo!

Ich mache beruflich Geometrie-Software. Geometrische Tests mit Double-
Koordinaten sind nicht zuverlässig. Aber es gibt einige Möglichkeiten:

Verwende einen exakten Zahlentyp. Das ist vermutlich am einfachsten,
aber es ist auch sehr langsam. Also bei einer Million Tests würde ich
es nicht mehr so machen. http://gmplib.org/
Unter www.cgal.org findest Du eine hochwertige Geometrie-Library. Die
Basiskomponenten, die Du brauchst, stehen unter der LGPL. Der Aufwand,
das aufzusetzen und sich einzuarbeiten, ist allerdings nicht gering.
Schau Dir http://www.cs.cmu.edu/~quake/robust.html an. Der Source
Code ist in C geschrieben, aber das sollte sich verwenden lassen.

lg

krümelkacker

Kann jemand mal das Problem genauer definieren? Ich verstehe gar nicht, was der OP hier berechnen will.