Reply to topic

Intersecting two hyperplanes in 3D space

antimatter
Registered Member
Posts
5
Karma
0
Hello,

I've got two Hyperplane<double, 3> instances h1 and h2. They are each constructed from 3 points with the static function Hyperplane<...>::Through(). But when I call
Code: Select all
h1.intersection(h2);

I get an assertion error with the macro THIS_METHOD_IS_ONLY_FOR_VECTORS_OF_A_SPECIFIC_SIZE at compile-time.
Is this intersection() method capable, of intersecting planes in a 3D space at all? From the API description, I doubt that. It returns a VectorType. But a line in 3D space needs at least two components (e.g. support vector and direction vector).

Is there an elegant way to calculate the intersection line of two intersecting planes in 3D space with Eigen?

Best regards,
Andreas
User avatar ggael
Moderator
Posts
2740
Karma
17
OS
As explained in the doc, this method is to intersect two lines.
antimatter
Registered Member
Posts
5
Karma
0
You are right, now I see it too. Sorry for being blind.

So, in order to find the support vector of the intersection line, I have implemented a common algorithm, I already used a while ago. There, both planes share a common support vector which will be the support vector of the intersection line. And they each have their own normal vectors.

Image

Assuming, that the z-component of the support vector will be 0, I use the x- and y-components of both normal vectors to build the equation system:

Image

Now, given 3 points for each plane, I construct a plane with the Hyperplane::Through() function and get their normal vectors by calling the normal() method. According to doc, this returns a unit normal vector of that plane. When I use these for calculating the support vector, I get wrong results.
Another way is the plane equation ax+by+cz=d. Setting d to 1, and putting those points in a 3x3 coefficient matrix A, I can solve Ax=b and get the normal vector (a,b,c) for that plane. Using these normal vectors for calculating the support vector, I get the correct results.

Why? What is the difference between the normal vectors from Hyperplane::normal() and my Ax=b method ?
User avatar ggael
Moderator
Posts
2740
Karma
17
OS
hm, the support vector is simply the cross product of the two normals.
antimatter
Registered Member
Posts
5
Karma
0
I'm pretty sure, that the cross product of the normals is a vector that is parallel to the line of intersection and is therefore the directional vector.
Imagine a plane, that is identical with the y-z-plane of a coordinate system (z pointing up). It's normal is parallel with the x-axis.
Imagine a nother plane identical with the x-z-plane, which normal is parallel with the y-axis.
Then it's easy to imagine, that the line of intersection is identical to the z-axis, which is normal to both x- and y-axis.
User avatar ggael
Moderator
Posts
2740
Karma
17
OS
sorry, I did not really read your post and interpreted "support vector" as the direction vector. If you're looking for an arbitrary point on the intersection line, then fixing z=0 is a pretty bad idea as the other components might tends to infinite values in some cases. I suggest you to use a SVD to pick the minimal norm solution, which is the orthogonal projection of the origin on the intersection line. Note that the normal of a Eigen::HyperPLane is a unit vector, and you thus have to consider the offset() values (http://eigen.tuxfamily.org/dox/classEig ... plane.html). Again, fixing it to 1 would be a rather bad idea since this does not allow to represent planes passing through the origin.
antimatter
Registered Member
Posts
5
Karma
0
Hi,
ok, SVD is unknown to me. I researched for it, and found something useful for my problem. This is without fixing any coefficients to predefined numbers. http://www.math.usm.edu/lambers/mat419/lecture15.pdf

So, after creating the Hyperplanes from 3 points, I get the normal(), which is the coefficients a, b, c and a negated offset() is d in ax+by+cz=d. That way I can define the underdetermined equation system Ax=b. Following the linked PDF, this will calculate the minimum norm solution:

Image

In Eigen:
Code: Select all
Vector3d x = A.transpose() * ((A * A.transpose()).inverse() * b);

So far, it delivers the results I expect. Is that what you meant? If not, maybe you can show some short example or a better way, how to do that in Eigen.

Merci and best regards,
Andreas
User avatar ggael
Moderator
Posts
2740
Karma
17
OS
yes, that's one way. The drawback is that method requires that rank(A)==rows(A). A more general approach is to use SVD:

x = A.jacobiSvd(ComputeFullU|ComputeFullV).solve(b);
antimatter
Registered Member
Posts
5
Karma
0
Very nice. Thank you. :-)

 
Reply to topic

Bookmarks



Who is online

Registered users: askinner, Baidu [Spider], Bing [Bot], davidemme, Google [Bot], google01103, Ian_Fin, jribeiro, linuxfluesterer, Mamarok, wolfi323, Yahoo [Bot]