Archived

This topic is now archived and is closed to further replies.

Area of polygon - Problems

This topic is 5623 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

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.

Share this post


Link to post
Share on other sites
Ok, after some extreme debugging, I found that the problems lies in the code that determines whether the line is outside the polygon. I was only writing code with this case in mind ->

B .---.
/ /
A / /
.-------. /
| C /
| /
.--------.

The code generates a "right" normal. So if I had a vector that started at A and went to B, it would catch that it was outside the polygon, because C would be in front based on the "right" normal. However, that''s assuming that point A is the first one in the list. What if point B is the first one? When the code comes around to check to see if the line from point B to point A is valid, it passes, beacuse the right normal faces away from the polygon and point C is behind it, which counts as a pass.

Belive me, it took diligent file output to catch this one However, now I must rethink a robust solution to this problem. Hmmm...

Share this post


Link to post
Share on other sites
The first obvious answer that comes to mind is to simply perform the check in both directions (AB then BA) before considering that the line is acceptable.

The other way to deal with this problem is to represent the polygon as a sequence of line-angle pairs. You can instantly detect if a vertex-vertex line falls outside the polygon by checking the angles at all intervening vertices along the shortest path around the polygon. So, for example, the polygon you drew might be given by the sequence:

2,90,3,60,2+sqrt(5),150,1,60,-2,-60,3,90

where I went anticlockwise from A and used anticlockwise as the positive direction for angles. Also, I assigned arbitrary lengths. AC was length 3 and the height of the base rectangle at the bottom left of the poly was 2.

The negative angle corresponds to the angle ACB. Clearly the line AB lies outside the polygon in this situation, as does the line AD, where I'm making D the top-right most vertex. If the next vertex around is E, then AE lies inside the polygon since there is a shorter path going from A to E in the anticlockwise direction that has no negative angles.

Furthermore, you can quickly check using this encoding that your polygon is convex by checking the sign of each angle. If they're all positive then you can implement a quicker area algorithm for convex polygons.

Have fun,

Timkin


[edited by - Timkin on July 25, 2002 4:35:08 AM]

Share this post


Link to post
Share on other sites
Here's another way to find a triangle area. Just find the crossproduct between 2 of its sides and take 1/2 of its length. IMHO, it's easier:


    
float TriangleArea(D3DVECTOR& a,D3DVECTOR& b,D3DVECTOR& c)
{
D3DVECTOR ab = b -a;
D3DVECTOR ac = c-a;
D3DVECTOR cross(ab.y * ac.z - ab.z * ac.y,
ab.z * ac.x - ab.x * ac.z,
ab.x * ac.y - ab.y * ac.x);

return 0.5f * sqrtf(cross.x*cross.x+cross.y*cross.y+cross.z*cross.z);
}


Oh. while at it, the formula can be readily generalized to find an area of an arbitrary (convex or concave) polygon. The only condition is that the polygon has to be flat.


  
float PolyArea(D3DVECTOR* vertex, int vertexCount)
{
D3DVECTOR v0 = vertex[0];
D3DVECTOR acc(0,0,0);
for ( int i = 1; i < vertexCount-1; i++ )
{
D3DVECTOR r1 = vertex[i]-v[0];
D3DVECTOR r2 = vertex[i+1]-v[0];
D3DVECTOR cross(r1.y *r2.z - r1.z*r2.y,
r1.z *r2.x - r1.x*r2.z,
r1.x *r2.y - r1.y*r2.x);
acc+=cross;
}
return 0.5f * sqrtf(acc.x*acc.x+acc.y*acc.y+acc.z*acc.z);
}


Best,
M


[edited by - microsha on July 25, 2002 11:26:24 AM]

[edited by - microsha on July 25, 2002 11:26:46 AM]

Share this post


Link to post
Share on other sites
quote:
Oh. while at it, the formula can be readily generalized to find an area of an arbitrary (convex or concave) polygon. The only condition is that the polygon has to be flat.

That works great, and it trashes my whole exercise Oh well, learn something everyday. I completely forgot that the length of the crossproduct was the area of the parallelagram. Doh

Share this post


Link to post
Share on other sites