Punkt in Dreieck (3D)

[TDM]

Grüße.

Ich versuche momentan rauszufinden ob ein Punkt P sich in einem Dreieck ABC befindet.

Zweidimensional klappt das schon ganz gut, allerdings hab ich Probleme bei der dritten Dimension.

Die Idee ist wie folgt:

Durch P wird eine Gerade gezogen die das Dreieck an zwei Seiten schneidet. Die Schnittpunkte seien s und t.

Folglich gilt mit den Vektoren a = Strecke[A,B] und b = Strecke[A,C] für den Ortsvektor p = Strecke[o, P]:

p = Strecke[0, A] + s * a + t * b

Der Punkt liegt drin, wenn: 0<= s,t <= 1 UND 0<= s+ t <= 1

(die Gleichheitszeichen gelten für den Rand)

Zweidimensional ist das alles kein Problem, denn da gilt:

I. px = ax + s*(bx -ax ) +t*(cx - ax)

II. py = ay + s*(by -ay) + t*(cy - ay)

Das wird wundertoll gelößt:

I'

II''.

II'' in I'.

Folglich sieht meine Funktion (C#) dann so aus:

private static bool isPointinTriangle(Point3D p, Point3D a, Point3D b, Point3D c)

        {

            double s =  -((c.X - a.X) * p.Z + (a.Z - c.Z) * p.X + a.X * c.Z - a.Z * c.X) /

                        ((b.X - a.X) * c.Z + (a.Z - b.Z) * c.X + a.X * b.Z - a.Z * b.X);


            double t = (-((b.X - a.X) * s - p.X + a.X) / (c.X - a.X));


            return ((0 < s) && (s <= 1) && (0 < t) && (t <= 1) && (0 < (s + t)) && ((s + t) <= 1));


        }

[TDM]

Das geht auch solange ich nicht die Y-Koordinate beachte.

So, nun soll das ganze aber im dreidimensionalen Raum gelten.

Die Idee einfach den Punkt zu erweitern, dass gilt:

III. pz = az + s * (bz - az) + t * (cz - az)

machte es bisschen unübersichtlich, aber naja.

Folglich gilt nun:

III' in II''.

und III'. in II'' in I'

[TDM]

Dann gilt für die Funktion:

        private static bool isPointinTriangle(Point3D p, Point3D a, Point3D b, Point3D c)

        {

            double s =  -(((b.X - a.X) * c.Y - a.Y * b.X + a.X * a.Y) * p.Z + ((a.X - b.X) * c.Z + (b.Z - a.Z) * c.X - a.X * b.Z + a.Z * b.X) * p.Y + ((a.Z - b.Z) * c.Y + a.Y * b.Z - a.Y * a.Z) * p.X + (a.Y * b.X - a.X * a.Y) * c.Z + (a.X * b.Z - a.Z * b.X) * c.Y + (a.Y * a.Z - a.Y * b.Z) * c.X) /

                        (((b.X - a.X) * b.Y - a.Y * b.X + a.X * a.Y) * c.Z + ((a.Y - b.Y) * b.Z + a.Z * b.Y - a.Y * a.Z) * c.X + (a.X * b.Y - a.X * a.Y) * b.Z - a.Z * b.X * b.Y + a.Y * a.Z * b.X);


            double t = (-((b.X - a.X) * s - p.X + a.X) / (c.X - a.X));


            return ((0 < s) && (s <= 1) && (0 < t) && (t <= 1) && (0 < (s + t)) && ((s + t) <= 1));


        }

Soweit zur Theorie. Getestet wird momentan mit folgendem Aufruf:

            Point3D a = new Point3D(-1, 1, 1);

            Point3D b = new Point3D(1, 1, 1);

            Point3D c = new Point3D(1, 1, -1);


            isPointinTriangle(new Point3D(0.25, 1, 0.25), a, b, c);

Problem an der ganzen Sache ist, dass mein s bei dieser Konstellation leider immer außerhalb des Wertebereichs ist (NaN).

Deswegen meine Frage: Kann jemand den Algorithmus nachvollziehen und mir sagen, wo mein (Denk-/Tipp-) Fehler liegt?

Memo@Mods:

Links zu lang -> Überschreitet immer 10.000 Zeichen, folglich 3 Posts.

[Klotzkopp]

Du hast nicht bedacht, dass in R3 der Fall auftreten kann, dass P gar nicht in der Ebene des Dreiecks liegt. In diesem Fall gibt es keine Lösung für s und t, da du ja durch keine Linearkombination aus dieser Ebene herauskommst.

Das müsstest du also vorher prüfen.

Das hier könnte interessant für dich sein:

Point in triangle test

Bedenke bitte auch, dass dir die Ungenauigkeit der Fließkommatypen hier ordentlich in die Suppe spucken kann.

[TDM]

Erstmal danke für den Link, les ich mir gleich mal durch.

Wegen der Ebene:

Von welcher reden wir ? x, y oder z ?

Das Problem an der Sache ist nämlich, dass ein Dreieck auch gekippt liegen kann.

Z.B. {(-1, 2, 3);(1, 5, 5);(0, -4, -5)}

Böse Sache das.

Aber ich hatte ja die Punkte erstmal so gewählt, dass der Punkt theoretisch drinnen liegen würde.

(Dreieck liegt komplett auf der Y-1er-Ebene)

Geht dennoch nicht.

Das mit der Ungenauigkeit hatte ich vorhin schon mit ein paar Kollegen besprochen. Ich werd das am Besten auf 1000tel Stellen runden, da stört das eh keinen.

[flashpixx]

Bestimme die Baryzentrischen Koordinaten und erzeuge damit eine Projektion des Dreiecks in eine Ebene. Damit kannst Du dann prüfen ob sich der Punkt innerhalb der Fläche befindet bzw. Du kannst auch damit Schnitttests zwischen Strahl und Dreieck bestimmen

Phil

[TDM]

Wenn ich mir Klotzkopp's Link so anschaue, war es wahrscheinlich genau das, was ich machen wollte. O.o'

Ich teste mal ein bisschen rum, mal sehen obs klappt.

[Klotzkopp]

Wegen der Ebene:

Von welcher reden wir ? x, y oder z ?

Von der, die durch das Dreieck selbst aufgespannt wird. Das Dreieck definiert eine Ebene im R3, und wenn dein Punkt gar nicht erst in dieser Ebene liegt, hat dein Gleichungssystem gar keine Lösung.

Aber ich hatte ja die Punkte erstmal so gewählt, dass der Punkt theoretisch drinnen liegen würde.

Auch dann ist deine Formel nicht "sicher". Schau dir die Formel für t an: Egal wo P liegt, wenn c.X == a.X ist, dividierst du da durch 0.

[TDM]

Auch dann ist deine Formel nicht "sicher". Schau dir die Formel für t an: Egal wo P liegt, wenn c.X == a.X ist, dividierst du da durch 0.

Das ist ja egal, weil in der Formel für t ja jeweils nur eine Achse beachtet wird.

Erst durch das einsetzen von t in die s-Gleichung werden mehrere Achsen beachtet.

Und dann muss ja folglich x oder y unterschiedlich sein, weil sonst würde C auf A liegen und es wäre kein Dreieck mehr.

[Klotzkopp]

Das ist ja egal, weil in der Formel für t ja jeweils nur eine Achse beachtet wird.

Das ist nicht egal, und die Begründung ist Unsinn.

Erst durch das einsetzen von t in die s-Gleichung werden mehrere Achsen beachtet.

Ich rede nicht von deinen Umformungen. Ich rede vom Endergebnis. Schau dir deinen Code doch an.

Du dividierst durch (c.X - a.X). Das bedeutet, dass du durch Null dividierst, wenn c.X gleich a.X ist. Daran gibt's nichts zu rütteln.

Deine Ausgangsformel für s hat ein ähnliches Problem. Beim Auflösen nach s hast du durch (b.X - a.X) dividiert. Das darfst du aber gar nicht, wenn b.X gleich a.X ist. Du hast diverse Umformungen gemacht, dabei aber anscheinend völlig außer Acht gelassen, dass dabei Ausnahmen auftreten, die du vorher ausschließen musst.

[TDM]

Ach, jetzt weiß ich wo du meinst.

Jup. War ein Fehler, danke.

Eigentlich wollte ich Vektoren der Strecken nehmen und nicht die Punkte selber. O.o