davideberly / GeometricTools

A collection of source code for computing in the fields of mathematics, geometry, graphics, image analysis and physics.
Boost Software License 1.0
1.08k stars 202 forks source link

Plane to plane intersection looks wrong #57

Closed maquemo closed 1 year ago

maquemo commented 1 year ago

When trying to calculate a plane to plane intersection with the below input, the origin point for the resulting line looks wrong:

Input: plane1: (origin=(0, 0, -6), normal vector=(1, 4, -1)), plane2:(origin(0, 0, 4), normal vector(1, -2, 1))

The point p=(5, 0, -1) is a known point that it is on the resulted intersecting line of the two planes so the next assert should be true:

gte::Vector3 vA{ 1.0, 4.0, -1.0 }; gte::Normalize(vA); gte::Vector3 vB{ 1.0, -2.0, 1.0 }; gte::Normalize(vB);

gte::Plane3 A{ vA, 6.0 }; gte::Plane3 B{ vB, 4.0 };

gte::FIQuery<double, gte::Plane3, gte::Plane3> query;

const auto result = query(A, B);

assert(gte::Length(gte::Cross((gte::Vector3{ 5, 0, -1 } - result.line.origin), result.line.direction)) < 0.00000

image

but it fails. The direction of the intersection line is ok, it is only the origin line.

In the image the point B is the returned origin point for the intersection.

davideberly commented 1 year ago

Your assert statement computes the length of a vector and asserts it is negative. That does not seem right to me. Regardless, I believe the heart of the problem is the following.

The comments at the beginning of Hyperplane.h indicate the plane normal must be unit length. You are using the plane constructor Plane3(N,C), where N is a plane normal and C is a plane constant. You have normalized your inputs to satisfy the unit-length-normal constraint. The planes are of the form Dot(N,X-C)=0. My code computes the intersection line origin to be P=(13.889289794715246, 3.7376117522089771, 3.3838926808354231). For unit-length plane normals N1 and N2, Dot(N1,P)-C1 = 4.4408910^{-15} and Dot(N2,P) -C2 = 1.7763610^{-15}).

Let Q=(5,0,-1). If you use N1=(1,4,-1) and N2 - (1,-2,1), which are NOT unit length, you get Dot(N1,Q)-C1=0 and Dot(N2,Q)-C2=0. You need to use the unit-length plane normals. In particular, when N1 and N2 are unit-length, Dot(N1,Q)-C1=19.4558 and Dot(N2,Q)-C2=5.79796, so Q is NOT a point on the intersection line.

davideberly commented 1 year ago

I forgot to mention. If your application has the plane defined as Dot(N,X) = C for non-unit-length N, create the GT plane using: T length = Normalize(N); Plane3 plane(N,C/length);

maquemo commented 1 year ago

Hello, thank you for the help.

Regarding the weird assert comparation, it is a typo, in the original code is 0.000001.

So then, because the normal vector has been normalizaed then the constants should be divided by the length of the vector before be normalized, that is right?

Thank you

davideberly commented 1 year ago

Yes, that is right. See my previous comment "I forgot to mention...".