Distance from point to plane

Hi everywhere!

A question; How can I compute the distance between a plane (in ellipsoidal coordinates) to a point (in ellipsoidal coordinates too)?

The plane have 4 points (the borders points), and I need calculate the closest distance from this plane to a point.

For example;

          %Lon.            %Lat.         %Depth (km)
Plane= [-74.65180      -37.74230        1.62620
        -72.86330      -33.14430        1.62620
        -70.88940      -33.64880       62.19350
        -72.67800      -38.24670       62.19350
        -74.65180      -37.74230        1.62620]

Point=[-72.7060  -37.7950]

Thank you so much!

Answer 1

[Matt J]

Well, the fact that you have ellipsoidal coordinates shouldn't matter, because you can convert them to Cartesian coordinates PlaneCart and PointCart. Once you've done so

      c=mean(PlaneCart);
      N=null(bsxfun(@minus, PlaneCart,c));
            N=mean(N,2); 

      distance=abs(dot(N,c-PointCart));

[Nicolás]

Thank you for answered.

I use ell2xyz.m to convert ellipsoidal to cartesian coordinates.

When I use your script, I have this problem;

N =

     Empty matrix: 3-by-0

%--------------------------------------------------------

I tried with this code;

Plane= [-37.7423  -74.6518    1.6262
        -33.1443  -72.8633    1.6262
        -33.6488  -70.8894   62.1935
        -38.2467  -72.6780   62.1935];

Point=[ -33.3210  -71.4110];

for i=1:4
[x(i),y(i),z(i)]=ell2xyz(deg2rad(Plane(i,1)),deg2rad(Plane(i,2)),-Plane(i,3)*1000);
end

P1=[x(1),y(1),z(1)]/1000; %en km
P2=[x(2),y(2),z(2)]/1000; %en km
P3=[x(3),y(3),z(3)]/1000; %en km
normal = cross(P1-P2, P1-P3);
d=dot(normal, P1); d=-d;
%Equation form: a*x+b*y+c*z+d=0

%Traspaso a coordenadas cartesianas
[xPo,yPo,zPo]=ell2xyz(deg2rad(Point(1)),deg2rad(Point(2)),0);
xPo=xPo/1000; yPo=yPo/1000; zPo=zPo/1000; %en km

d1=norm(normal(1)*xPo+normal(2)*yPo+normal(3)*zPo+d);
d2=sqrt(normal(1)^2+normal(2)^2+normal(3)^2);

d=(d1/d2);

but, the value of "d" is incorrect. Any idea why?

Best Regards.

[Matt J]

Looks like it should be

     d=abs(dot(normal, P1-Point))

[Nicolás]

not work :/

Thank you anyway.

[Matt J]

One solution would be to discard one of the points in PlaneCart so that you have only 3 points. You shouldn't need more to define the plane.

My original solution worked fine (I tested it), when you have "PlaneCart" contains 3 points.

[Nicolás]

Yes, but obtain similar "d" values with my script. This "d" values are wrong.

[Matt J]

Let's do this together with the following sample data

    PlaneCart=eye(3); PointCart=[0,0,0];

When I run with the above inputs, I obtain

   distance =

    0.5774

Do you not obtain the same thing? And do you disagree that this is the correct distance from the origin to the plane containing the rows of eye(3)? If you disagree, then you have given us the wrong description of your problem in some way.

Answer 2

[Thorsten]

This is a math question rather than a Matlab question. First convert your plane from 3-point-form to Hessian normal form http://mathworld.wolfram.com/Plane.html. Second compute the point-plane distance http://mathworld.wolfram.com/HessianNormalForm.html

[Nicolás]

Thank you for answered.

The mathematical problem is clear, the use of Hessian normal is useful for define the signs if the point is on the same side of the plane as the normal vector and negative if it is on the opposite side ( refer ). Therefore the next script is similar to use the Hessian normal form (but not work);

Plane= [-37.7423  -74.6518    1.6262
        -33.1443  -72.8633    1.6262
        -33.6488  -70.8894   62.1935
        -38.2467  -72.6780   62.1935];
Point=[ -33.3210  -71.4110];
for i=1:4
[x(i),y(i),z(i)]=ell2xyz(deg2rad(Plane(i,1)),deg2rad(Plane(i,2)),-Plane(i,3)*1000);
end
P1=[x(1),y(1),z(1)]/1000; %en km
P2=[x(2),y(2),z(2)]/1000; %en km
P3=[x(3),y(3),z(3)]/1000; %en km
normal = cross(P1-P2, P1-P3);
d=dot(normal, P1); d=-d;
%Equation form: a*x+b*y+c*z+d=0
%Traspaso a coordenadas cartesianas
[xPo,yPo,zPo]=ell2xyz(deg2rad(Point(1)),deg2rad(Point(2)),0);
xPo=xPo/1000; yPo=yPo/1000; zPo=zPo/1000; %en km
d1=norm(normal(1)*xPo+normal(2)*yPo+normal(3)*zPo+d);
d2=sqrt(normal(1)^2+normal(2)^2+normal(3)^2);
d=(d1/d2);

The problem is to implement in Matlab, because, the result "d" is wrong.

Regards.

Answer 3

[Thorsten]

 Plane= [-37.7423  -74.6518    1.6262
        -33.1443  -72.8633    1.6262
        -33.6488  -70.8894   62.1935
        -38.2467  -72.6780   62.1935];
Point=[ -33.3210  -71.4110];
for i=1:4
  [x(i),y(i),z(i)]=ell2xyz(deg2rad(Plane(i,1)),deg2rad(Plane(i,2)),-Plane(i,3)*1000);
end
P1=[x(1),y(1),z(1)]/1000; %en km
P2=[x(2),y(2),z(2)]/1000; %en km
P3=[x(3),y(3),z(3)]/1000; %en km
P4=[x(4) y(4) z(4)]/1000;
n = cross(P1-P2, P1-P3);
d = dot(-P1, n)
nnorm = n/norm(n); 
p = d/norm(n);

 Distance = dot(nnorm, P4) + p

 Distance =

    2.8967

That means that Point 4 is not in the plane, but has a distance of 2.89 km.

 [xPo,yPo,zPo]=ell2xyz(deg2rad(Point(1)),deg2rad(Point(2)),0);
Po = [xPo yPo zPo];

 Distance = dot(nnorm, Po) + p

Distance =

 -6.0717e+006

So point Po lies 6.07 km low the plane. Does this make sense?