Best fit of plane for polygon

Started by
4 comments, last by Dave Eberly 16 years, 3 months ago
Hi, I have some polygons wich I am checking to see that they are coplaner, I can determine the plane of the first 3 points but what I want to do is use all the points to determine the plane that best fits the polygon. Im a bit stuck as to how to do this, theres loads of info on how to do it for 3 points, and some sites just say its just simple analytical geometry, but i cant seem to find a good reference. I can think of some ways of doing this 1) solving simultaneous eq for all points so aXn+bYn+cZn+d=0 this is beyond me unfortunatly, but im sure some people here find this sort of thing easy. 2) summing the cross products of the normal of the line vectors at each vertex, im not sure if this would give me a best fit ? if it is concave would it still handle this? .... at the moment im just using the given normal and measuring the distance of each point from that plane, but I need to be able to tweak the points to make it coplaner, so I could do with taking an average. also rounding errors makes this a bit dubious for very large planes, or planes very far away from the origin. im using directx in c# wich limits me to float wich isnt ideal, unless I use another library wich has double vectors and convert back to draw with directx, any suggestions? the models are existing assets, some also seem to have some concave polygons. thanks.
Advertisement
this is what ive come up with,...
I normalise after summing so the length of the lines adds some weight to the plane. seems to agree with the supplied normal 99.9% of the time

        public Plane GetBestPlane()        {            Plane bestPlane=new Plane(0,0,0,0);            for (int i = 0; i < Points.Count; i++)            {                Vector3 p0 = Points;                Vector3 p1 = Points[(i + 1) % Points.Count];                Vector3 p2 = Points[(i + 2) % Points.Count];                Vector3 v1 = p1 - p0;                Vector3 v2 = p2 - p1;                Vector3 vc = Vector3.Cross(v1, v2);                bestPlane.Normal += vc;            }            bestPlane.Normal.Normalize();            foreach (Point3d p in Points)                bestPlane.D -= Vector3.Dot(p.Position, bestPlane.Normal);            bestPlane.D /= Points.Count;            return bestPlane;        }

I don't completely understand the goal. There are easier ways to test for planarity than to compute the plane of best fit. For example, separating the vertices into triplets, treating them as triangles and testing the normals for near-equality would do. However, if you need a plane-of-best-fit for an arbitrary set of points, you'll have to work a bit harder.

There is no single solution to this problem. For best accuracy you should reduce it to a least-squares problem under some appropriate set of criteria, but your method is fine for a quick-and-dirty approximation, given that the vertex density is roughly uniform. For better results, you could weigh the terms according to the edge-lengths, for example, but this may be overkill. Tell us more about the origin of the data and the degree of accuracy you require.

By the way, this probably fits better in the Maths and Physics forum.
Ring3 Circus - Diary of a programmer, journal of a hacker.
thanks, the goal is to validate model input data, I have a given normal wich I wish to make sure is close the given set of points, and that the points all lie close to a plane, then either adjust the points or adjust the normal.

the problem with using triplets is that they may be very short or even concave,
or a straight line, so I figure a better way is to measure the distance of each point from the plane, if there is an error then I need to know if the point needs adjusting or the normal.

the above function does in fact take into acount the length of the lines ast it sums the cross product before taking the normal so it is proportional to the lengths too, them then takes the normal of the sum, its proved quite an effective solution, in many places less points are in error with the generated normal then the given normal.

the problem with lines with nearly zero angle is taken care of by adding a central point and using a line from point 1 to the center then to point 2 etc..

a solution involving least squared erros would be the optimum but might be more involved and cpu intensive, I think this will do for now.
The matrix method I'm aware of for doing least-squares requires a function, so your plane equation aX + bY + cZ + d = 0 would have to be rewritten as X = -(b/a)Y - (c/a)Z - (d/a) or something like that where one variable is a function of the others. The problem with this approach is that either a, b, or c (depending on how you write the function) is assumed to be 1, so if your best-fit plane is one in which the wrong variable ('a' in the above example) would be zero or very close to it, things would blow up. It's the same problem you encounter using the equation "Y = mX + b" for a line instead of "aX + bY + c = 0". The first form assumes the 'Y' coefficient is 1, so you can't describe a line in which the normal has a zero 'Y' coefficient (which would make it vertical). And if you try, your slope goes to infinity, which presents numerous precision problems in and of itself.
For least-squares fitting a plane to a set of points, see Least Squares Fitting of Data. I also have an implementation at my web site. See the Foundation/Approximation files Wm4ApprPlaneFit3.{h,cpp}.

This topic is closed to new replies.

Advertisement