Defining normal unit vector for arbitrary plane surface in 3D space

Dear All,

I'm trying to find unit vector which pointing perpendicularly outward from arbitrary shape of panel in 3D space. I found from other similar question, that it can be done by calculating the cross product of the points in the panel. I did the same thing, but there is some error which I can't really understand. A tried to do it by using these three methods below:

if true
  %
A = [-225.0000    2.7555e-14         0];
B = [-225.0000    2.7404e-14   -2.8802e-15];
C = [-223.6800    24.3270         0];
D = [-223.6800    24.1940      -2.5428];
center = (A+B+C+D)/4
% Method1: Calculating unit vector by seeing the coordinate of center point of the panel as a vector from axis of origin then normalized it
UnitVector1=center/norm(center)
CheckA1=(UnitVector1(1,1)^2+UnitVector1(1,2)^2+UnitVector1(1,3)^2)^0.5
CheckA2=cross(center,UnitVector1)
max(abs(CheckA2))
% Method2: Obtaining unit vector by calculating cross product of the vectors at panel's corner point
UnitVector2=cross((D-A),(B-A));
UnitVector2=UnitVector2/norm(UnitVector2)
CheckB1=(UnitVector2(1,1)^2+UnitVector2(1,2)^2+UnitVector2(1,3)^2)^0.5
CheckB2=cross(center,UnitVector2)
max(abs(CheckB2))
% Method3: Obtaining unit vector by calculating cross product of the panel's center point
Point1=(A+B)/2;
Point2=(B+C)/2;
UnitVector3=cross((Point1-center),(Point2-center));
UnitVector3=UnitVector3/norm(UnitVector3)
CheckC1=(UnitVector3(1,1)^2+UnitVector3(1,2)^2+UnitVector3(1,3)^2)^0.5
CheckC2=cross(center,UnitVector3)
max(abs(CheckC2))

end

These panel coordinates are taken from sphere-shaped 3D bodies. It means that the center point of the panel itself, when being normalized, is already become the unit vector which pointing outward of the panel (method 1). when I do cross product between the obtained unit vector and the center coordinate, the result is going to be zero (achieved with first method with error precision 10^-16) because there is no area generated between those vectors. But, since I'm dealing not only with sphere-shaped object, I need to find another way which more applicable for other case (method 2 and 3). But apparently, the resulting cross product between obtained unit vector and the center point coordinates is not zero, which means it's not perfectly perpendicuar outward the panel.

Any idea to improve/solve this issue will be highly appreciated.

Regards,

Fredo Ferdian

[yashodeep]

%program to find unit norma to the surface

Answer 1

[Matt J]

Here's a general way to fit a plane normal to any number of points,

V=[A;B;C;D];
 V=V-mean(V); %assumes R2016b or later, otherwise use bsxfun()

 [U,S,W]=svd(V,0);

 normal=W(:,end)

[fredo ferdian]

Dear Matt J,

Thank you very much for your reply. I seems that your solution is based on some kind of advance math theory. Apparently, it is not really applicable for any point coordinates. To be honest, I just blindly using your code without understanding the basic concept of orthonormal basis/singular value decomposition, and put these coordinate data,

if true
A = [-225.0000    2.7555e-14         0];
B = [-225.0000    2.7404e-14   -2.8802e-15];
C = [-223.6800    24.3270         0];
D = [-223.6800    24.1940      -2.5428];
end

and got empty matrix. Quickly browse "how" and "why", seems like it is related about inverse-ability of the matrix and kernel-thing which I never heard before. Does it really important to learn about those advance math theory? I just need a simple unit vector pointing outward of the panel in 3D space. Or is there any other simpler way to do it? Maybe something about increasing the precision of cross product? I just need to obtain this result from my coordinate data above

if true
  % -0.9985    0.0540   -0.0028
end

Thank you very much before,

Regards,

Fredo Ferdian

[Matt J]

I've modified my code. It should work now.

[fredo ferdian]

Dear Matt J,

Thank you very much for your response. It works like a charm. Unfortunately, the obtained unit vector is exactly the same with my original code at method 2 and 3. Which means, when I re-checked it again by doing the cross product with the center point, it generates some values (not zero -> means not exactly pointing outward like my first method isn't it?). I do it as your suggestion like this

if true
% V=[A;B;C;D];
Vnew = bsxfun(@minus, V, mean(V))
[U,S,W]=svd(Vnew,0);
UnitVector4=(W(:,end))'
CheckD1=(UnitVector4(1,1)^2+UnitVector4(1,2)^2+UnitVector4(1,3)^2)^0.5
CheckD2=cross(mean(V),UnitVector4)
end

And get this unit vector

if true
  %  -0.9985    0.0542   -0.0028
end

This is a bit weird. The values of 'CheckD2' is still quite large. I'm expecting something like the first method (CheckA2 -> order 10^-15)

Do you might have some idea what's happening here?

Thank you very much before.

Regards,

Fredo Ferdian

[LC]

Matt J,

Is the normal always pointing outward?

[Matt J]

@LC, the only guarantee is that the normal is perpendicular to the fitted plane. The code doesn't know which of the two possible directions you consider to be "outward".

[Diaa]

Is it possible to find from your code the basis vectors of the plane passing through the points without using the function pca?

[Matt J]

@Diaa Once you've calculated the normal, you can obtain basis vectors for the plane by doing

null(normal(:).')

[Diaa]

If you don't mind, would you please help me know how to calculate the direction cosines of the output normal vector?

Thanks in advance.

[Matt J]

That is what my posted answer gives you.

[Diaa]

In terms of perfromance, is it better to make it

V = V-mean(V);
[U,S,W] = svd(V,0);

or in one step using

[U,S,W] = svd(V-mean(V),0);

I am talking about tens of this operation, so is nesting the calculations always a good idea in favor of performance instead of overwriting the existent variables (e.g. V here)?

[Matt J]

A very minor improvement, I would guess, but maybe the profiler would tell you differently.

Answer 2

[Matt J]

Use planarFit() in this FEX download,

https://www.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-for-fitting-conics-and-quadrics