Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Expand fill of convex polygon

I have a convex polygon P1 of N points. This polygon could be any shape or proportion (as long as it is still convex).

I need to compute another polygon P2 using the original polygons geometry, but "expanded" by a given number of units. What might the algorithm be for expanding a convex polygon?

like image 956
Adam Harte Avatar asked Sep 20 '10 08:09

Adam Harte


People also ask

What is the use of convex hull?

Convex hulls have wide applications in mathematics, statistics, combinatorial optimization, economics, geometric modeling, and ethology. Related structures include the orthogonal convex hull, convex layers, Delaunay triangulation and Voronoi diagram, and convex skull.

What is convex hull problem?

Computing the convex hull is a problem in computational geometry. The indices of the points specifying the convex hull of a set of points in two dimensions is given by the command ConvexHull[pts] in the Wolfram Language package ComputationalGeometry` .

Is a convex hull a polygon?

The convex hull of a simple polygon is itself a convex polygon. Overlaying the original simple polygon onto its convex hull partitions this convex polygon into regions, one of which is the original polygon. The remaining regions are called pockets.


2 Answers

Shapes with their inflated equivalents

To expand a convex polygon, draw a line parallel to each edge and the given number of units away. Then use the intersection points of the new lines as the vertices of the expanded polygon. The javascript/canvas at the end follows this functional breakdown:

Step 1: Figure out which side is "out"

The order of the vertices (points) matters. In a convex polygon, they can be listed in a clockwise (CW), or a counter-clockwise (CCW) order. In a CW polygon, turn one of the edges 90 degrees CCW to obtain an outward-facing normal. On a CCW polygon, turn it CW instead.

CW and CCW polygons

If the turn direction of the vertices is not known in advance, examine how the second edge turns from the first. In a convex polygon, the remaining edges will keep turning in the same direction:

  1. Find the CW normal of the first edge. We don't know yet whether it's facing inward or outward.

  2. Compute the dot product of the second edge with the normal we computed. If the second edge turns CW, the dot product will be positive. It will be negative otherwise.

Dot product to find turn direction

Math:

// in vector terms:
v01 = p1 - p0                      // first edge, as a vector
v12 = p2 - p1                      // second edge, as a vector
n01 = (v01.y, -v01.x)              // CW normal of first edge
d = v12 * n01                      // dot product

// and in x,y terms:
v01 = (p1.x-p0.x, p1.y-p0.y)       // first edge, as a vector
v12 = (p2.x-p1.x, p2.y-p1.y)       // second edge, as a vector
n01 = (v01.y, -v01.x)              // CW normal of first edge
d = v12.x * n01.x + v12.y * n01.y; // dot product: v12 * n01

if (d > 0) {
    // the polygon is CW
} else {
    // the polygon is CCW
}

// and what if d==0 ?
// -- that means the second edge continues in the same
// direction as a first.  keep looking for an edge that
// actually turns either CW or CCW.

Code:

function vecDot(v1, v2) {
    return v1.x * v2.x + v1.y * v2.y;
}

function vecRot90CW(v) {
    return { x: v.y, y: -v.x };
}

function vecRot90CCW(v) {
    return { x: -v.y, y: v.x };
}

function polyIsCw(p) {
    return vecDot(
        vecRot90CW({ x: p[1].x - p[0].x, y: p[1].y - p[0].y }),
        { x: p[2].x - p[1].x, y: p[2].y - p[1].y }) >= 0;
}

var rot = polyIsCw(p) ? vecRot90CCW : vecRot90CW;

Step 2: Find lines parallel to the polygon edges

Now that we know which side is out, we can compute lines parallel to each polygon edge, at exactly the required distance. Here's our strategy:

  1. For each edge, compute its outward-facing normal

  2. Normalize the normal, such that its length becomes one unit

  3. Multiply the normal by the distance we want the expanded polygon to be from the original

  4. Add the multiplied normal to both ends of the edge. That will give us two points on the parallel line. Those two points are enough to define the parallel line.

Parallel line by adding a weighted normal vector

Code:

// given two vertices pt0 and pt1, a desired distance, and a function rot()
// that turns a vector 90 degrees outward:

function vecUnit(v) {
    var len = Math.sqrt(v.x * v.x + v.y * v.y);
    return { x: v.x / len, y: v.y / len };
}

function vecMul(v, s) {
    return { x: v.x * s, y: v.y * s };
}

var v01 = { x: pt1.x - pt0.x, y: pt1.y - pt0.y };  // edge vector
var d01 = vecMul(vecUnit(rot(v01)), distance);     // multiplied unit normal
var ptx0 = { x: pt0.x + d01.x, y: pt0.y + d01.y }; // two points on the
var ptx1 = { x: pt1.x + d01.x, y: pt1.y + d01.y }; //     parallel line

Step 3: Compute the intersections of the parallel lines

--these will be the vertices of the expanded polygon.

intersection of expanded polygon edges

Math:

A line going through two points P1, P2 can be described as:

P = P1 + t * (P2 - P1)

Two lines can be described as

P = P1 + t * (P2 - P1)
P = P3 + u * (P4 - P3)

And their intersection has to be on both lines:

P = P1 + t * (P2 - P1) = P3 + u * (P4 - P3)

This can be massaged to look like:

(P2 - P1) * t + (P3 - P4) * u = P3 - P1

Which in x,y terms is:

(P2.x - P1.x) * t + (P3.x - P4.x) * u = P3.x - P1.x
(P2.y - P1.y) * t + (P3.y - P4.y) * u = P3.y - P1.y

As the points P1, P2, P3 and P4 are known, so are the following values:

a1 = P2.x - P1.x    a2 = P2.y - P1.y
b1 = P3.x - P4.x    b2 = P3.y - P4.y
c1 = P3.x - P1.x    c2 = P3.y - P1.y

This shortens our equations to:

a1*t + b1*u = c1
a2*t + b2*u = c2    

Solving for t gets us:

t = (b1*c2 - b2*c1)/(a2*b1 - a1*b2)

Which lets us find the intersection at P = P1 + t * (P2 - P1).

Code:

function intersect(line1, line2) {
    var a1 = line1[1].x - line1[0].x;
    var b1 = line2[0].x - line2[1].x;
    var c1 = line2[0].x - line1[0].x;

    var a2 = line1[1].y - line1[0].y;
    var b2 = line2[0].y - line2[1].y;
    var c2 = line2[0].y - line1[0].y;

    var t = (b1*c2 - b2*c1) / (a2*b1 - a1*b2);

    return {
        x: line1[0].x + t * (line1[1].x - line1[0].x),
        y: line1[0].y + t * (line1[1].y - line1[0].y)
    };
}

Step 4: Deal with special cases

There is a number of special cases that merit attention. Left as an exercise to the reader...

  1. When there's a very sharp angle between two edges, the expanded vertex can be very far from the original one. You might want to consider clipping the expanded edge if it goes beyond some threshold. At the extreme case, the angle is zero, which suggests that the expanded vertex is at infinity, causing division by zero in the arithmetic. Watch out.

  2. When the first two edges are on the same line, you can't tell if it's a CW or a CCW polygon by looking just at them. Look at more edges.

  3. Non convex polygons are much more interesting... and are not tackled here.

Full sample code

Drop this in a canvas-capable browser. I used Chrome 6 on Windows. The triangle and its expanded version should animate.

browser snapshot

canvas { border: 1px solid #ccc; }
$(function() {
      var canvas = document.getElementById('canvas');  
      if (canvas.getContext) {  
          var context = canvas.getContext('2d');  

          // math for expanding a polygon

          function vecUnit(v) {
              var len = Math.sqrt(v.x * v.x + v.y * v.y);
              return { x: v.x / len, y: v.y / len };
          }

          function vecMul(v, s) {
              return { x: v.x * s, y: v.y * s };
          }

          function vecDot(v1, v2) {
              return v1.x * v2.x + v1.y * v2.y;
          }

          function vecRot90CW(v) {
              return { x: v.y, y: -v.x };
          }

          function vecRot90CCW(v) {
              return { x: -v.y, y: v.x };
          }

          function intersect(line1, line2) {
              var a1 = line1[1].x - line1[0].x;
              var b1 = line2[0].x - line2[1].x;
              var c1 = line2[0].x - line1[0].x;

              var a2 = line1[1].y - line1[0].y;
              var b2 = line2[0].y - line2[1].y;
              var c2 = line2[0].y - line1[0].y;

              var t = (b1*c2 - b2*c1) / (a2*b1 - a1*b2);

              return {
                  x: line1[0].x + t * (line1[1].x - line1[0].x),
                  y: line1[0].y + t * (line1[1].y - line1[0].y)
              };
          }

          function polyIsCw(p) {
              return vecDot(
                  vecRot90CW({ x: p[1].x - p[0].x, y: p[1].y - p[0].y }),
                  { x: p[2].x - p[1].x, y: p[2].y - p[1].y }) >= 0;
          }

          function expandPoly(p, distance) {
              var expanded = [];
              var rot = polyIsCw(p) ? vecRot90CCW : vecRot90CW;

              for (var i = 0; i < p.length; ++i) {

                  // get this point (pt1), the point before it
                  // (pt0) and the point that follows it (pt2)
                  var pt0 = p[(i > 0) ? i - 1 : p.length - 1];
                  var pt1 = p[i];
                  var pt2 = p[(i < p.length - 1) ? i + 1 : 0];

                  // find the line vectors of the lines going
                  // into the current point
                  var v01 = { x: pt1.x - pt0.x, y: pt1.y - pt0.y };
                  var v12 = { x: pt2.x - pt1.x, y: pt2.y - pt1.y };

                  // find the normals of the two lines, multiplied
                  // to the distance that polygon should inflate
                  var d01 = vecMul(vecUnit(rot(v01)), distance);
                  var d12 = vecMul(vecUnit(rot(v12)), distance);

                  // use the normals to find two points on the
                  // lines parallel to the polygon lines
                  var ptx0  = { x: pt0.x + d01.x, y: pt0.y + d01.y };
                  var ptx10 = { x: pt1.x + d01.x, y: pt1.y + d01.y };
                  var ptx12 = { x: pt1.x + d12.x, y: pt1.y + d12.y };
                  var ptx2  = { x: pt2.x + d12.x, y: pt2.y + d12.y };

                  // find the intersection of the two lines, and
                  // add it to the expanded polygon
                  expanded.push(intersect([ptx0, ptx10], [ptx12, ptx2]));
              }
              return expanded;
          }

          // drawing and animating a sample polygon on a canvas

          function drawPoly(p) {
              context.beginPath();
              context.moveTo(p[0].x, p[0].y);
              for (var i = 0; i < p.length; ++i) {
                  context.lineTo(p[i].x, p[i].y);
              }
              context.closePath();
              context.fill();
              context.stroke();
          }

          function drawPolyWithMargin(p, margin) {
              context.fillStyle = "rgb(255,255,255)";  
              context.strokeStyle = "rgb(200,150,150)";  
              drawPoly(expandPoly(p, margin));

              context.fillStyle = "rgb(150,100,100)";  
              context.strokeStyle = "rgb(200,150,150)";  
              drawPoly(p);
          }

          var p = [{ x: 100, y: 100 }, { x: 200, y: 120 }, { x: 80, y: 200 }];
          setInterval(function() {
              for (var i in p) {
                  var pt = p[i];
                  if (pt.vx === undefined) {
                      pt.vx = 5 * (Math.random() - 0.5);
                      pt.vy = 5 * (Math.random() - 0.5);
                  }

                  pt.x += pt.vx;
                  pt.y += pt.vy;

                  if (pt.x < 0 || pt.x > 400) { pt.vx = -pt.vx; }
                  if (pt.y < 0 || pt.y > 400) { pt.vy = -pt.vy; }
              }
              context.clearRect(0, 0, 800, 400);
              drawPolyWithMargin(p, 10);
          }, 50);
      }
  });
<html>
    <head>
            <script type="text/javascript" src="http://ajax.googleapis.com/ajax/libs/jquery/1.4.2/jquery.min.js"></script>
    </head>
    <body>
        <canvas id="canvas" width="400" height="400"></canvas>  
    </body>
</html>

sample code disclaimers:

  • the sample sacrifices some efficiency for the sake of clarity. In your code, you may want to compute each edge's expanded parallel just once, and not twice as in here

  • the canvas's y coordinate grows downward, which inverts the CW/CCW logic. Things keep on working though as we just need to turn the outward normals in a direction opposite to the polygon's -- and both get flipped.

like image 170
Oren Trutner Avatar answered Oct 30 '22 03:10

Oren Trutner


For each line segment of the original, find the midpoint m and (unit length) outward normal u of the segment. The corresponding segment of the expanded polygon will then lie on the line through m+n*u (where you want to expand the original by n) with normal u. To find the vertices of the expanded polygon you then need to find the intersection of pairs of successive lines.

like image 39
dmuir Avatar answered Oct 30 '22 01:10

dmuir