Polygon clipping - a little elaboration?

Question

I have been reading a lot about the sutherland hodgman polygon clipping algorithm and understand the general idea. However, when I see the actual implementation of it (like the one below), I get confused about the coordinate comparisons such as those in the intersection and inside methods. Therefore, I was wondering if someone could elaborate on what and why? I see a ton of videos and articles explaining the general concepts but I really have trouble finding some explanation of the actual details regarding the implementation.

  bool inside(b2Vec2 cp1, b2Vec2 cp2, b2Vec2 p) {
      return (cp2.x-cp1.x)*(p.y-cp1.y) > (cp2.y-cp1.y)*(p.x-cp1.x);
  }

  b2Vec2 intersection(b2Vec2 cp1, b2Vec2 cp2, b2Vec2 s, b2Vec2 e) {
      b2Vec2 dc( cp1.x - cp2.x, cp1.y - cp2.y );
      b2Vec2 dp( s.x - e.x, s.y - e.y );
      float n1 = cp1.x * cp2.y - cp1.y * cp2.x;
      float n2 = s.x * e.y - s.y * e.x;
      float n3 = 1.0 / (dc.x * dp.y - dc.y * dp.x);
      return b2Vec2( (n1*dp.x - n2*dc.x) * n3, (n1*dp.y - n2*dc.y) * n3);
  }

  //http://rosettacode.org/wiki/Sutherland-Hodgman_polygon_clipping#JavaScript
  //Note that this only works when fB is a convex polygon, but we know all 
  //fixtures in Box2D are convex, so that will not be a problem
  bool findIntersectionOfFixtures(b2Fixture* fA, b2Fixture* fB, vector<b2Vec2>& outputVertices)
  {
      //currently this only handles polygon vs polygon
      if ( fA->GetShape()->GetType() != b2Shape::e_polygon ||
           fB->GetShape()->GetType() != b2Shape::e_polygon )
          return false;

      b2PolygonShape* polyA = (b2PolygonShape*)fA->GetShape();
      b2PolygonShape* polyB = (b2PolygonShape*)fB->GetShape();

      //fill subject polygon from fixtureA polygon
      for (int i = 0; i < polyA->GetVertexCount(); i++)
          outputVertices.push_back( fA->GetBody()->GetWorldPoint( polyA->GetVertex(i) ) );

      //fill clip polygon from fixtureB polygon
      vector<b2Vec2> clipPolygon;
      for (int i = 0; i < polyB->GetVertexCount(); i++)
          clipPolygon.push_back( fB->GetBody()->GetWorldPoint( polyB->GetVertex(i) ) );

      b2Vec2 cp1 = clipPolygon[clipPolygon.size()-1];
      for (int j = 0; j < clipPolygon.size(); j++) {
          b2Vec2 cp2 = clipPolygon[j];
          if ( outputVertices.empty() )
              return false;
          vector<b2Vec2> inputList = outputVertices;
          outputVertices.clear();
          b2Vec2 s = inputList[inputList.size() - 1]; //last on the input list
          for (int i = 0; i < inputList.size(); i++) {
              b2Vec2 e = inputList[i];
              if (inside(cp1, cp2, e)) {
                  if (!inside(cp1, cp2, s)) {
                      outputVertices.push_back( intersection(cp1, cp2, s, e) );
                  }
                  outputVertices.push_back(e);
              }
              else if (inside(cp1, cp2, s)) {
                  outputVertices.push_back( intersection(cp1, cp2, s, e) );
              }
              s = e;
          }
          cp1 = cp2;
      }

      return !outputVertices.empty();
  }

(code stolen from iforce2d :) )

Answer 1

You say you understand the general idea, presumably from reading something like Sutherland Hodgman Algorithm. That explains at a high level exactly what inside and intersection do.

As to the details of how they achieve their objectives, that is all just straight up text book linear algebra.

inside is testing the sign of (cp2 - cp2) cross (p - cp1) and returning true iff the sign is strictly greater than zero. You could rewrite the return statement as:

return (cp2.x-cp1.x)*(p.y-cp1.y) - (cp2.y-cp1.y)*(p.x-cp1.x) > 0;

by moving the second term to the left of the > which gives you exactly the cross product on the left.

Note that a cross product is typically a vec3 cross vec3 operation and requires computation of all three terms. However, we're doing this in 2d meaning the vec3s have the form (x, y, 0). Therefore we only need to compute the z output term, since the cross must be perpendicular to the xy plane, and therefore be of the form (0, 0, value).

intersection finds the point at which two vectors intersect using exactly the algorithm listed here: Line Intersection from Two Points. In particular, we care about the formula immediately following the text "The determinants can be written out as:"

In the context of that formula n1 is (x1 y2 - y1 x2), n2 is (x3 y4 - y3 x4) and n3 is 1 / ((x1 - x2) (y3 - y4) - (y1 - y2) (x3 - x4))

-- EDIT --

To cover the issue raised in the comments, here is as full an explanation as I can give for why the return value from inside() is a test of the sign of the cross product.

I'm going to go off on a slight tangent, show my age, and note that the cross product formula has a very simple memory aid. You just need to remember the first magic word from Woods and Crowther's text adventure game Colossal Cave. xyzzy.

If you have two vectors in three dimensions: (x1, y1, z1) and (x2, y2, z2), their cross product (xc, yc, zc) is evaluated thus:

xc = y1 * z2 - z1 * y2;
yc = z1 * x2 - x1 * z2;
zc = x1 * y2 - y1 * x2;

Now, look at the first line, remove the c, 1 and 2 suffixes from the terms, all the spaces and operators and just look at the remaining letters. It's the magic word. Then you just go vertically down, replacing x with y, y with z and z with x as you go from line to line.

Getting back to business, both the terms on the right of the expansions of xc and yc contain either z1 or z2. But we know that both of these are zero since our input vectors are in the xy plane, and therefore have a zero z component. That's why we can completely elide computing those two terms, because we know they'll be zero.

This is 100% consistent with the definition of what the cross product does, the resulting vector is always perpendicular to both input vectors. Hence if both input vectors are in the xy plane, we know the output vector must be perpendicular to the xy plane, and therefore have the form (0, 0, z)

So, what do we have for the z term?

zc = x1 * y2 - y1 * x2;

in this case vector 1 is cp2-cp1 and vector 2 is p-cp1. So plugging that into the above we get:

zc = (cp2.x-cp1.x)*(p.y-cp1.y) - (cp2.y-cp1.y)*(p.x-cp1.x);

But as noted, we don't care about its value, only it's sign. We want to know if that is greater than zero. Hence:

return (cp2.x-cp1.x)*(p.y-cp1.y) - (cp2.y-cp1.y)*(p.x-cp1.x) > 0;

which is then rewritten as:

return (cp2.x-cp1.x)*(p.y-cp1.y) > (cp2.y-cp1.y)*(p.x-cp1.x);

Finally, what does the sign of that term have to do with whether the point p is inside or outside the clipping polygon? You are quite correct that all the clipping takes place in the 2d xy plane, so why are we involving 3d operations at all?

The important thing to realize is that the cross product formula in 3d is not commutative. The order of the two vector operands is significant, in terms of the angle between them. The first image on the Wikipedia page for Cross Product shows it perfectly. In that diagram, if you look down from above, when evaluating a cross b, the shortest angular direction from a to b is counter-clockwise. In that instance, that leads to a cross product with a positive z value, assuming positive z goes up the page. However if you evaluate b cross a, the shotest angular distance from b to a is clockwise, and the cross product has a negative z value.

Thinking back to the Wikipedia page for the algorithm itself, you've got the blue "clipping" line that works its way counter-clockwise around the clipping polygon. If you think of that vector as always having positive magnitude in the counter-clockwise direction it'll always be cp2 - cp1 for any pair of adjacent vertices in the clipping polygon.

Keeping this in mind, imagine what you'd see if you stood at cp1, with your nose pointed straight at cp2. The interior of the clipping polygon will be on your left, and the exterior on the right. Now consider two points p1 and p2. We'll say p1 is inside the clipping poly, and p2 is outside. That means that the quickest way to point yourt nose at p1 is to rotate counter-clockwise, and the quickest way to point at p2 is to rotate clockwise.

So by studying the sign of the cross product, we're really asking 'Do we rotate clockwise or counter-clockwise from the current edge to look at the point' which equates to asking if the point is inside the clipping polygon or outside.

I'll add one final suggestion. If you're at all interested in this sort of thing, or 3d rendering, or any programming that involves modelling the mathematical representation of the real world, taking a good solid course in linear algebra that covers the likes of cross products, dot products, vectors, matrices and the interactions between them all will be one of the best things you can ever do. It will provide a very strong foundation for a fair amount of what is done with computers.