Schnitttest(Gerde, Gerade) im R^3 und Schittpunktberechnung
Robert Hartmann
Der Test, ob sich zwei beliebige Geraden im R^3 schneiden,
ist recht einfach implementiert.
Gegeben zwei Geraden A und B.
A: A.origin + lambda* A.direction
B: B.origin + kappa * B.direction
A.origin, B.origin, A.direction, B.direction sind Vektoren des R^3,
lambda und kappa sind aus R.
Weiterhin gilt, dass weder A.direction noch B.direction
identisch mit dem Nullvektor sind.
Der Test auf Existenz eines Schnittpunktes mï¿œsste sein:
Falls
Richtungsvektoren der Geraden sind linear unabhï¿œngig, d.h.
crossProd(A.direction, B.direction)
!= (0,0,0)
und
Strahlen sind nicht windschief, d.h.
spatProd(A.direction, B.direction, (B.origin - A.origin))
==
dotProd(crossProd(A.direction, B.direction),(B.origin - A.origin))
== 0
dann
schneiden sich die Strahlen
sonst
nicht.
Implementation mit C++
bool Ray3D::isIntersecting(const Ray3D& A, const Ray3D& B)
{
assert((A.direction!=Vector3D(0,0,0))&&(B.direction!=Vector3D(0,0,0)));
Vector3D cP = crossProd(A.direction, B.direction);
Vector3D oD = B.origin - A.origin;
if (( cP != Vector3D(0,0,0) ) && ( dotProd(cP,oD) == 0 ))
return true;
else
return false;
}
Bei der Implementation einer Methode fï¿œr die
Schnittpunktberechnung komme ich irgendwie zu
keinem schï¿œn ï¿œbersichtlichen Ergebnis.
Wenn ich die Geradengleichungen gleichsetze
A.origin + lambda* A.direction
= B.origin + kappa * B.direction
dann erkenne ich in Komponentenschreibweise
A.origin.x+ lambda * A.direction.x
= B.origin.x + kappa * B.direction.x
A.origin.y+ lambda * A.direction.y
= B.origin.y + kappa * B.direction.y
A.origin.z+ lambda * A.direction.z
= B.origin.z + kappa * B.direction.z
dass man auf jeden Fall, um lambda
und kappa ausrechnen zu kï¿œnnen,
Fallunterscheidungen braucht;
denn durch Null zu dividieren
ist ja nicht gerade zulï¿œssig.
Betrachten wir den Fall
(A.direction.x!=0)&&(B.direction.y!=0)
etwas genauer:
Es gilt nach Umstellung der Komponentenschreibweise:
lambda
= (B.origin.x + kappa * B.direction.x - A.origin.x)
/ A.direction.x
= (B.origin.x - A.origin.x) / A.direction.x
+ (kappa * B.direction.x) / A.direction.x
und
kappa
= (A.origin.y+ lambda * A.direction.y - B.origin.y)
/ B.direction.y
= (A.origin.y - B.origin.y) / B.direction.y
+ lambda * A.direction.y / B.direction.y
Um den Wert fï¿œr lambda zu bestimmen,
setze ich nun kappa in lambda ein und lï¿œse nach lambda auf:
lambda
= (B.origin.x - A.origin.x) / A.direction.x
+ (kappa * B.direction.x) / A.direction.x
= (B.origin.x - A.origin.x) / A.direction.x
+ ((A.origin.y - B.origin.y) / B.direction.y
+ lambda * A.direction.y / B.direction.y)
* B.direction.x) / A.direction.x
damit ist, wenn ich mich nicht verrechnet habe:
lambda
= ( [(B.origin.x-A.origin.x)/A.direction.x]
+ ( [(A.origin.y-B.origin.y)/B.direction.y]
*[B.direction.x/A.direction.x ) )
/ ( 1
- [A.direction.y*B.direction.x]
/ [B.direction.y*A.direction.x] )
diesen Wert in die kappa-Berechnung einsetzen.
Wenn mit den so berechneten Werten fï¿œr
lambda und kappa die bisher ungenutzte
dritte Komponentengleichung zur wahren
Aussage wird
A.origin.z+ lambda * A.direction.z
== B.origin.z + kappa * B.direction.z ,
kann man endlich den Schnittpunkt zurï¿œckliefern.
Dies ist wie geschrieben nur der Fall
(A.direction.x!=0)&&(B.direction.y!=0) gewesen.
ï¿œhnliches mï¿œsste man nun fï¿œr die anderen Fï¿œlle,
dazu darf man dann aber auch keinen der Fï¿œlle
vergessen, auch machen.
Nun zu meiner eigentlichen Frage:
Kann man eine allgemeine Schnittpunktbestimmung
nicht irgendwie deutlich einfacher implementieren?
Dazu sei weiterhin angemerkt,
dass der Zugriff auf die Vektor-Komponenten
auch mit numerischen Indizes zugreifen kann:
So ist z.B.
A.origin.x == A.origin[0],
A.origin.y == A.origin[1],
A.origin.z == A.origin[2].
entsprechend ï¿œhnlich fï¿œr die anderen Vektoren
A.direction, B.origin, B.direction
Gruￜ Robert
Jakob Creutzig
> Der Test, ob sich zwei beliebige Geraden im R^3 schneiden,
> ist recht einfach implementiert.
[..]
> Bei der Implementation einer Methode f�r die
> Schnittpunktberechnung komme ich irgendwie zu
> keinem sch�n �bersichtlichen Ergebnis.
>
> Wenn ich die Geradengleichungen gleichsetze
>
> A.origin + lambda* A.direction
> = B.origin + kappa * B.direction
hast Du ein lineares Gleichungssystem mit zwei
Unbekannten in drei Dimensionen:
(lamba, -kappa) * (A.direction, B.direction)^T
= B.origin - A.origin
Quasi der Grossvater aller linearen Gleichungssysteme ;-).
fup2 dsm
Best,
Jakob
Markus Wichmann
Ich habe mal über das Problem meditiert: Erstmal, üm nachzusehen, ob es
überhaupt einen Schnittpunkt gibt, muss das LGS ja lösbar sein. Die
Lösbarkeitsbedingung ist, dass der Rang der Koeffizientenmatrix dem der
erweiterten Koeffizientenmatrix gleicht. Den Rang der
Koeffizientenmatrix bekommst du leicht raus: Der Einfachheit halber sei
a.d der Richtungsvektor der Geraden a. Dann ist die Koeffizientenmatrix
( a.d.x -b.d.x )
A = ( a.d.y -b.d.y )
( a.d.z -b.d.z )
Auf diese Matrix soll der Gaussche Algorithmus angewendet werden. Also
erstmal sehen, ob a.d.x==0. Denn wenn, dann muss man sich was einfallen
lassen, wie man die Null da weg bekommt. Da der Richtungsvektor niemals
der Nullvektor sein darf, reicht es aus, wenn man eine der beiden
anderen Zeilen zur ersten addiert. Das gleiche muss man gegebenenfalls
für b.d.y tun.
Dann kommen wir zur Rangermittlung. Da A21, A31 und A32 sowieso Null
werden, interessiert eigentlich nur A22: Wird A22 Null, wenn man die
Zeile so umformt, dass A21 Null wird? Wenn ja, ist der Rang der
Koeffizientenmatrix nämlich nicht 2, sondern 1. Kleiner kann er AFAIK
nur werden, wenn auch A11 und A12 Null werden, wenn man diese
Umformungen abgeschlossen hat.
Wir fassen in Quelltext zusammen:
#define EPSILON (1e-20)
bool inline is_null_vector(Vector3d v)
{
return ((abs(v.x) < EPSILON)
&& (abs(v.y) < EPSILON)
&& (abs(v.z) < EPSILON));
}
bool has_common_point(Line3d& a, Line3d& b)
{
double[3][2] A = {
{a.direction.x, -b.direction.x},
{a.direction.y, -b.direction.y},
{a.direction.z, -b.direction.z}
};
double[3][3] Aa = {
{a.direction.x, -b.direction.x, b.origin.x-a.origin.x},
{a.direction.y, -b.direction.y, b.origin.y-a.origin.y},
{a.direction.z, -b.direction.z, b.origin.z-a.origin.z}
};
int rank1, rank2;
assert(!is_null_vector(a.direction) && !is_null_vector(b.direction));
if (abs(A[0][0]) < EPSILON)
{
if (abs(A[1][0]) > EPSILON)
{
A[0][0] += A[1][0];
A[0][1] += A[1][1];
Aa[0][2] += Aa[1][2];
}
else
{
A[0][0] += A[2][0];
A[0][1] += A[2][1];
Aa[0][2] += Aa[2][2];
}
Aa[0][0] = A[0][0];
Aa[0][1] = A[0][1];
}
if (abs(A[1][1]) < EPSILON)
{
if (abs(A[0][1]) > EPSILON)
{
A[1][0] += A[0][0];
A[1][1] += A[0][1];
Aa[1][2] += Aa[0][2];
}
else
{
A[1][0] += A[2][0];
A[1][1] += A[2][1];
Aa[1][2] += Aa[2][2];
}
Aa[1][0] = A[1][0];
Aa[1][1] = A[1][1];
}
/**
* Jetzt kommt der witzige Teil:
* Um A[1][0] auf Null zu setzen, muss ich nach Gausschem Algo die
* A[0]-Zeile temporär mit -1/A[1][0] multiplizieren. Wenn A[1][0]
* schon 0 ist, kann man sich die Arbeit natürlich sparen. In dem
* Fall haben wir gerade dafür gesorgt, dass A[1][1] ungleich 0 ist,
* was natürlich heißt, dass der Rang dieser Matrix gleich 2 ist.
*/
if (abs(A[1][0]) < EPSILON)
rank1 = 2;
else
{
/* Wenn man A[1][1] auf die gleiche Weise behandelt wie A[1][0],
* dh. mit A[1][0]* -1/A[0][1] addiert,
* und dann 0 rauskommt, ist der Rang gleich 1, sonst 2
*/
if (abs(A[1][1] + A[1][0] * -1/A[0][1]) < EPSILON)
rank1 = 1;
else
rank1 = 2;
}
/* Nun zur erweiterten Koeffizientenmatrix. Hier wird es schon
* kniffliger, den Rang zu bestimmen. Aber glücklicherweise kann man
* den Schnittpunkt gleich mit bestimmen, wenn man den Gausschen
* Algorithmus angewendet hat.
*/
/* Erstmal muss ich nachsehen, ob alle drei Elemente der letzten
* Spalte 0 sind. In dem Fall kann ich mir den Rest sparen: Es gibt
* einen Schnittpunkt, und der ist gleich dem Ursprung der beiden
* Geraden. Denn auch vektoriell gilt: a - b = 0 <=> a = b.
*/
Vector3d v = {
.x = Aa[0][2];
.y = Aa[1][2];
.z = Aa[2][2];
};
if (is_null_vector(v))
return true;
for (int i = 1; i <= 2; i++)
{
for (int j = 0; j <= 2; j++)
{
Aa[i][j] -= Aa[i-1][j]/Aa[i-1][i-1];
}
}
/*Wenn ich mich nicht vertan habe, war das so ziemlich der Gaussche
* Algorithmus.
* Jetzt muss Aa[2][2] == 0, sonst ist das LGS unlösbar
*/
if (abs(Aa[2][2]) > EPSILON)
{
rank2 = 3;
}
else
/* Es gibt noch Hoffnung... zunächst den Rang fertig bestimmen.
*/
if (abs(Aa[1][1]) < EPSILON) && (abs(Aa[1][2]) < EPSILON)
rank2 = 1;
else
rank2 = 2;
if (rank2 == rank1)
return true;
else return false;
}
So, den Rest überlasse ich dir. Noch ein paar Tipps zum Abschluss: Wenn
der Rang kleiner ist als die Anzahl der unbekannten (Sprich, der Rang
ist 1), dann gibt es unendlich viele Lösungen, was bezogen auf das
Problem bedeuten dürfte, dass die Geraden identisch sind. Und wenn die
Ränge gleich 2 sind, dann sieht die erweiterte Koeffizientenmatrix
irgendwie so aus:
(a b c)
(0 d e)
(0 0 0)
In Gleichungen also
I. a*s + b*t = c
II. d*t = e
Da du schon weißt, dass das LGS lösbar ist, kannst du es auch einfach
bei der Bestimmung von t bewenden lassen und daraus den Schnittpunkt
ausrechnen.
HTH,
Markus
--
GUI - ein Hintergrundbild und zwölf XTerms
vim -c "exec \"norm iwHFG#NABGURE#IVZ#UNPXRE\"|%s/#/ /g|norm g??g~~"
Robert Hartmann
Robert Hartmann schrieb am 12.06.2009 im Beitrag mit
Message-ID: <h0tda5$2nho$1...@ariadne.rz.tu-clausthal.de>
[Followup an de.sci.mathematik gesetzt]
> Hallo zusammen,
>
> Der Test, ob sich zwei beliebige Geraden im R^3 schneiden,
> ist recht einfach implementiert.
[...]
>
> Bei der Implementation einer Methode fï¿œr die
> Schnittpunktberechnung komme ich irgendwie zu
> keinem schï¿œn ï¿œbersichtlichen Ergebnis.
Nur wegen der Vollstï¿œndigkeit will ich hier nun
meine C++ Implementation der Schnittberechnung zweier
beliebiger euklidischer Geraden im R^3 angeben.
Dabei ist "Vector3D" eine Klasse fï¿œr die Vektorrechnung im R^3,
und "euklGerade3D" ist eine Klasse, die eben euklidische Geraden
im R^3 implementiert.
Gruￜ Robert
//Test auf Schnitt:
bool isIntersecting
(const euklGerade3D& A, const euklGerade3D& B) {
//Die Richtungsvektoren der Geraden mï¿œssen _sinnvoll_ sein,
// d.h. sie dï¿œrfen nicht dem Nullvektor entsprechen.
assert( (A.direction!=Vector3D(0,0,0) )
&&( B.direction!=Vector3D(0,0,0)) );
/*Falls
Richtungsvektoren der euklid. Geraden sind linear unabhï¿œngig,
d.h. crossProd(A.direction, B.direction) != (0,0,0)
und
euklid. Geraden sind nicht windschief,
d.h. spatProd(A.direction, B.direction, (B.origin - A.origin))
== dotProd( crossProd(A.direction, B.direction),
(B.origin - A.origin))
== 0
dann
schneiden sich die euklidischen Geraden
sonst
nicht.
*/
Vector3D cP = crossProd(A.direction, B.direction);
Vector3D oD = B.origin - A.origin;
if ((cP != Vector3D(0, 0, 0)) && (dotProd(cP, oD) == 0))
return true;
else
return false;
}
//Schnittpunktberechnung
Vector3D getIntersection
(const euklGerade3D& A, const euklGerade3D& B) {
assert(isIntersecting(A,B));
/*
A.origin + lambda* A.direction
== B.origin + kappa * B.direction
*/
double lambda = 0;
double kappa = 0;
if ((A.direction.x != 0) && (B.direction.y != 0)) {
lambda = (B.direction.y * (B.origin.x - A.origin.x))
+ (B.direction.x
* (A.origin.y - B.origin.y));
lambda /= (A.direction.x * B.direction.y)
- (B.direction.x * A.direction.y);
kappa = ((A.origin.y - B.origin.y)
+ (lambda * A.direction.y))
/ B.direction.y;
assert( (A.origin.z+ lambda * A.direction.z)
== (B.origin.z + kappa * B.direction.z));
} else if ((A.direction.x != 0) && (B.direction.z != 0)) {
lambda = (B.direction.z * (B.origin.x - A.origin.x))
+ (B.direction.x
* (A.origin.z - B.origin.z));
lambda /= (A.direction.x * B.direction.z)
- (B.direction.x * A.direction.z);
kappa = ((A.origin.z - B.origin.z)
+ (lambda * A.direction.z))
/ B.direction.z;
assert( (A.origin.y+ lambda * A.direction.y)
== (B.origin.y + kappa * B.direction.y));
} else if ((A.direction.y != 0) && (B.direction.z != 0)) {
lambda = (B.direction.z * (B.origin.y - A.origin.y))
+ (B.direction.y
* (A.origin.z - B.origin.z));
lambda /= (A.direction.y * B.direction.z)
- (B.direction.y * A.direction.z);
kappa = ((A.origin.z - B.origin.z)
+ (lambda * A.direction.z))
/ B.direction.z;
assert( (A.origin.x+ lambda * A.direction.x)
== (B.origin.x + kappa * B.direction.x));
} else if ((A.direction.y != 0) && (B.direction.x != 0)) {
lambda = (B.direction.x * (B.origin.y - A.origin.y))
+ (B.direction.y
* (A.origin.x - B.origin.x));
lambda /= (A.direction.y * B.direction.x)
- (B.direction.y
* A.direction.x);
kappa = ((A.origin.x - B.origin.x)
+ (lambda * A.direction.x))
/ B.direction.x;
assert( (A.origin.z+ lambda * A.direction.z)
== (B.origin.z + kappa * B.direction.z));
} else if ((A.direction.z != 0) && (B.direction.x != 0)) {
lambda = (B.direction.x * (B.origin.z - A.origin.z))
+ (B.direction.z
* (A.origin.x - B.origin.x));
lambda /= (A.direction.z * B.direction.x)
- (B.direction.z
* A.direction.x);
kappa = ((A.origin.x - B.origin.x)
+ (lambda * A.direction.x))
/ B.direction.x;
assert( (A.origin.y+ lambda * A.direction.y)
== (B.origin.y + kappa * B.direction.y));
} else if ((A.direction.z != 0) && (B.direction.y != 0)) {
lambda = (B.direction.y * (B.origin.z - A.origin.z))
+ (B.direction.z
* (A.origin.y - B.origin.y));
lambda /= (A.direction.z * B.direction.y)
- (B.direction.z
* A.direction.y);
kappa = ((A.origin.y - B.origin.y)
+ (lambda * A.direction.y))
/ B.direction.y;
assert( (A.origin.x+ lambda * A.direction.x)
== (B.origin.x + kappa * B.direction.x));
}
return (A.origin + (lambda * A.direction));
}