While reading a post a while back about using spheres to approximate the area of a mesh (or at least I think that''s what the question was), I decided, "Hey, I feel like designing a function that finds the area of any polygon, concave or convex, big or small". Three days later, everything was finished. I drew the most complex polygon I could think of and plugged it into my system for testing. Result? 95% accuracy. Not good.
What it does is break up the polygon into a collection of triangles and find the area of the triangles using Heron''s Formula. It uses a recursive tree approach - we search for a good line to divide the polygon with. A "good" line is one which doesn''t intersect any sides of the polygon and doesn''t cross outside the polygon. This line creates two smaller polygons, each of which is passed back into the recursive function. When we finnaly reach a triangle, it''s added to a list, and we return (it''s basically a BSP-tree generation algorithm).
For starters, some pseudo-code so you know what''s happening:
void PolygonBSP(Polygon& poly)
{
if(poly.IsTriangle()
{
Add_To_Tri_List(poly);
return;
}
for(start = poly.Points.begin(); loop through point list)
{
for(end = start; loop throuh point list)
{
if(same point or adjacent point)
continue;
Line line(start, end);
if(''line'' is "bad")
continue;
// Both polygons include line points
Polygon front(all points on one side of line);
Polygon back(all points on other side of line);
PolygonBSP(front);
PolygonBSP(back);
goto out;
}
}
out:;
}
I''m having some trouble with the algorithm however. I''ve modified the program to print out the final triangle list before adding the areas, and it appears that sometimes the code choses a line that not only intersects the polygon somewhere else, but also passes outside the polygon. Like I mentioned above, it choses a "good" line 95% of the time, but the other 5% of the time it choses a "bad" one.
For now, I will show you the function that determines whether or not a line is "good", as well as any supporting functions. Hopefully you might be able to spot something. I wanted to switch my line intersection code over to parametric, but I''m so familiar with vectors that I decided to implement them for this. However, if they''re not working right, or parametrics work best for this situation, then I''ll gladly reimplement things:
// Return NULL if no intersection, or the point if an intersection
Vector3* Line::Intersect(const Line& line)
{
// HACK - Use a static intersection point
static Vector3 intersection;
intersection = line.Get_First_Point();
// Check points
if((Get_First_Point() == line.Get_First_Point() &&
Get_Second_Point() == line.Get_Second_Point()) ||
(Get_First_Point() == line.Get_Second_Point() &&
Get_Second_Point() == line.Get_First_Point()))
return &intersection
Vector3 v1 = DeltaVector();
Vector3 v2 = line.DeltaVector();
// Cross product produces vector normal to both v1 and v2
Vector3 cp = Math::CrossProduct(v1, v2);
// Create vector perpendicular to v1 (this vecotor)
Vector3 normal = Math::CrossProduct(cp, v1);
normal.Normalize();
// distance of plane in equation
float distance = Math::DotProduct(Get_First_Point(), normal);
// Get dot product of both points in vec to our line/plane
float dp1 = Math::DotProduct(line.Get_First_Point(), normal) - distance;
float dp2 = Math::DotProduct(line.Get_Second_Point(), normal) - distance;
// Same sign indicates same side
if(!Math::DiffSign(dp1, dp2))
return 0;
// Otherwise, find intersection!
float height = float(fabs(dp1) + fabs(dp2));
float ratio = fabs(dp1) / height;
Vector3 intersect_delta = v2 * ratio;
intersection = line.Get_First_Point() + intersect_delta;
// Now test if this is a false-positive
normal = Math::CrossProduct(cp, v2);
normal.Normalize();
distance = Math::DotProduct(line.Get_First_Point(), normal);
dp1 = Math::DotProduct(Get_First_Point(), normal) - distance;
dp2 = Math::DotProduct(Get_Second_Point(), normal) - distance;
return (Math::DiffSign(dp1, dp2)) ? &intersection : 0;
}
bool Polygon::Intersect(const Line& line)
{
// Return true if this line intersects our polygon at any point
LineListType line_list = Get_Line_List();
LineListIteratorType it(line_list.begin());
for(; it != line_list.end(); ++it)
{
if((*it).Intersect(line))
return true;
}
return false;
}
// Determines if we have a valid line - doesn''t intersect, inside polygon
bool PolygonBSP::IsValidLine(Polygon& poly, const Line& line, const Polygon::PointListIteratorType& start, const Polygon::PointListIteratorType& end)
{
if(poly.Intersect(line))
return false;
// Check to see if the line goes through the polygon
Polygon::PointListIteratorType it;
// Create right normal for line
Vector3 forward = line.DeltaVector();
Vector3 up = Vector3(0.0f,0.0f,1.0f);
Vector3 normal = Math::CrossProduct(forward, up);
normal.Normalize();
float distance = Math::DotProduct(normal, line.Get_First_Point());
for(it = start; it != end; ++it)
{
if((Math::DotProduct(normal, **it) - distance) > 0.0f)
return false;
}
return true;
}
The intersection code might have a few issues (i.e dealing with colinear lines that intersect on another line), but considering the information we are working with, I don''t know how much they matter.
If you''ve read this far, I thank you
Hopefully you will quickly spot a logical error on my part.
