toolkit/javascript/d3/src/geom/hull.js
changeset 47 c0b4a8b5a012
equal deleted inserted replaced
46:efd9c589177a 47:c0b4a8b5a012
       
     1 /**
       
     2  * Computes the 2D convex hull of a set of points using Graham's scanning
       
     3  * algorithm. The algorithm has been implemented as described in Cormen,
       
     4  * Leiserson, and Rivest's Introduction to Algorithms. The running time of
       
     5  * this algorithm is O(n log n), where n is the number of input points.
       
     6  *
       
     7  * @param vertices [[x1, y1], [x2, y2], …]
       
     8  * @returns polygon [[x1, y1], [x2, y2], …]
       
     9  */
       
    10 d3.geom.hull = function(vertices) {
       
    11   if (vertices.length < 3) return [];
       
    12 
       
    13   var len = vertices.length,
       
    14       plen = len - 1,
       
    15       points = [],
       
    16       stack = [],
       
    17       i, j, h = 0, x1, y1, x2, y2, u, v, a, sp;
       
    18 
       
    19   // find the starting ref point: leftmost point with the minimum y coord
       
    20   for (i=1; i<len; ++i) {
       
    21     if (vertices[i][1] < vertices[h][1]) {
       
    22       h = i;
       
    23     } else if (vertices[i][1] == vertices[h][1]) {
       
    24       h = (vertices[i][0] < vertices[h][0] ? i : h);
       
    25     }
       
    26   }
       
    27 
       
    28   // calculate polar angles from ref point and sort
       
    29   for (i=0; i<len; ++i) {
       
    30     if (i === h) continue;
       
    31     y1 = vertices[i][1] - vertices[h][1];
       
    32     x1 = vertices[i][0] - vertices[h][0];
       
    33     points.push({angle: Math.atan2(y1, x1), index: i});
       
    34   }
       
    35   points.sort(function(a, b) { return a.angle - b.angle; });
       
    36 
       
    37   // toss out duplicate angles
       
    38   a = points[0].angle;
       
    39   v = points[0].index;
       
    40   u = 0;
       
    41   for (i=1; i<plen; ++i) {
       
    42     j = points[i].index;
       
    43     if (a == points[i].angle) {
       
    44       // keep angle for point most distant from the reference
       
    45       x1 = vertices[v][0] - vertices[h][0];
       
    46       y1 = vertices[v][1] - vertices[h][1];
       
    47       x2 = vertices[j][0] - vertices[h][0];
       
    48       y2 = vertices[j][1] - vertices[h][1];
       
    49       if ((x1*x1 + y1*y1) >= (x2*x2 + y2*y2)) {
       
    50         points[i].index = -1;
       
    51       } else {
       
    52         points[u].index = -1;
       
    53         a = points[i].angle;
       
    54         u = i;
       
    55         v = j;
       
    56       }
       
    57     } else {
       
    58       a = points[i].angle;
       
    59       u = i;
       
    60       v = j;
       
    61     }
       
    62   }
       
    63 
       
    64   // initialize the stack
       
    65   stack.push(h);
       
    66   for (i=0, j=0; i<2; ++j) {
       
    67     if (points[j].index !== -1) {
       
    68       stack.push(points[j].index);
       
    69       i++;
       
    70     }
       
    71   }
       
    72   sp = stack.length;
       
    73 
       
    74   // do graham's scan
       
    75   for (; j<plen; ++j) {
       
    76     if (points[j].index === -1) continue; // skip tossed out points
       
    77     while (!d3_geom_hullCCW(stack[sp-2], stack[sp-1], points[j].index, vertices)) {
       
    78       --sp;
       
    79     }
       
    80     stack[sp++] = points[j].index;
       
    81   }
       
    82 
       
    83   // construct the hull
       
    84   var poly = [];
       
    85   for (i=0; i<sp; ++i) {
       
    86     poly.push(vertices[stack[i]]);
       
    87   }
       
    88   return poly;
       
    89 }
       
    90 
       
    91 // are three points in counter-clockwise order?
       
    92 function d3_geom_hullCCW(i1, i2, i3, v) {
       
    93   var t, a, b, c, d, e, f;
       
    94   t = v[i1]; a = t[0]; b = t[1];
       
    95   t = v[i2]; c = t[0]; d = t[1];
       
    96   t = v[i3]; e = t[0]; f = t[1];
       
    97   return ((f-b)*(c-a) - (d-b)*(e-a)) > 0;
       
    98 }