speed up convex hull by filtering out points inside quadrilateral formed by 4 extreme points
961825
function fastConvexHull(points) {
    var left = points[0];
    var top = points[0];
    var right = points[0];
    var bottom = points[0];

    // find the leftmost, rightmost, topmost and bottommost points
    for (var i = 0; i < points.length; i++) {
        var p = points[i];
        if (p[0] < left[0]) left = p;
        if (p[0] > right[0]) right = p;
        if (p[1] < top[1]) top = p;
        if (p[1] > bottom[1]) bottom = p;
    }

    // filter out points that are inside the resulting quadrilateral
    var cull = [left, top, right, bottom];
    var filtered = cull.slice();
    for (i = 0; i < points.length; i++) {
        if (!pointInPolygon(points[i], cull)) filtered.push(points[i]);
    }

    // get convex hull around the filtered points
    var indices = convexHull(filtered);

    // return the hull as array of points (rather than indices)
    var hull = [];
    for (i = 0; i < indices.length; i++) hull.push(filtered[indices[i]]);
    return hull;
}
1070270
def convexhull(data, col, fill=True, point_size=4):
    
    from geoplotlib.layers import ConvexHullLayer
    _global_config.layers.append(ConvexHullLayer(data, col, fill, point_size))
1316092
def bound_perturbed_gmres(pseudo, p, epsilon, deltas):
    '''
    '''
    if not numpy.all(numpy.array(deltas) > epsilon):
        raise ArgumentError('all deltas have to be greater than epsilon')

    bound = []
    for delta in deltas:
        # get boundary paths
        paths = pseudo.contour_paths(delta)

        # get vertices on boundary
        vertices = paths.vertices()

        # evaluate polynomial
        supremum = numpy.max(numpy.abs(p(vertices)))

        # compute bound
        bound.append(epsilon/(delta-epsilon)
                     * paths.length()/(2*numpy.pi*delta)
                     * supremum)
    return bound
1321876
def convex_hull(self):
        
        points = util.vstack_empty([m.vertices for m in self.dump()])
        hull = convex.convex_hull(points)
        return hull
1290749
function _quadRetrievePoint(point, quad) {
    if (quad.level < quad.maxLevel) {
      var index = _quadIndex(point, quad.bounds);

      // If node does not exist we return an empty list
      if (quad.nodes[index] !== undefined) {
        return _quadRetrievePoint(point, quad.nodes[index]);
      }
      else {
        return [];
      }
    }
    else {
      return quad.elements;
    }
  }
636150
private static List<Point> simplifyRadialDist(List<Point> points, double sqTolerance) {
    Point prevPoint = points.get(0);
    ArrayList<Point> newPoints = new ArrayList<>();
    newPoints.add(prevPoint);
    Point point = null;

    for (int i = 1, len = points.size(); i < len; i++) {
      point = points.get(i);

      if (getSqDist(point, prevPoint) > sqTolerance) {
        newPoints.add(point);
        prevPoint = point;
      }
    }

    if (!prevPoint.equals(point)) {
      newPoints.add(point);
    }
    return newPoints;
  }
334356
def compute_convex_hull(feed: "Feed") -> Polygon:
    
    m = sg.MultiPoint(feed.stops[["stop_lon", "stop_lat"]].values)
    return m.convex_hull
226848
function trimCanvasXY(xy) {
    if( xy.length === 0 ) return;
    var last = xy[xy.length-1], i, point;

    var c = 0;
    for( i = xy.length-2; i >= 0; i-- ) {
        point = xy[i];
        if( Math.abs(last.x - point.x) === 0 && Math.abs(last.y - point.y) === 0 ) {
            xy.splice(i, 1);
            c++;
        } else {
            last = point;
        }
    }

    if( xy.length <= 1 ) {
        xy.push(last);
        c--;
    }
}
150495
def __regpoly_coords(self, x0, y0, x1, y1, sides, start, extent):
        """"""

        coords = []
        if extent == 0:
            return coords

        xm = (x0 + x1) / 2.
        ym = (y0 + y1) / 2.
        rx = xm - x0
        ry = ym - y0

        n = sides
        if n == 0: # 0 sides => circle
            n = round((rx + ry) * .5)
            if n < 2:
                n = 4

        # Extent can be negative
        dirv = 1 if extent > 0 else -1
        if abs(extent) > 360:
            extent = dirv * abs(extent) % 360

        step = dirv * 360 / n
        numsteps = 1 + extent / float(step)
        numsteps_int = int(numsteps)

        i = 0
        while i < numsteps_int:
            rad = (start - i * step) * DEG2RAD
            x = rx * math.cos(rad)
            y = ry * math.sin(rad)
            coords.append((xm+x, ym-y))
            i += 1

        # Figure out where last segment should end
        if numsteps != numsteps_int:
            # Vecter V1 is last drawn vertext (x,y) from above
            # Vector V2 is the edge of the polygon
            rad2 = (start - numsteps_int * step) * DEG2RAD
            x2 = rx * math.cos(rad2) - x
            y2 = ry * math.sin(rad2) - y

            # Vector V3 is unit vector in direction we end at
            rad3 = (start - extent) * DEG2RAD
            x3 = math.cos(rad3)
            y3 = math.sin(rad3)

            # Find where V3 crosses V1+V2 => find j s.t.  V1 + kV2 = jV3
            j = (x*y2 - x2*y) / (x3*y2 - x2*y3)

            coords.append((xm + j * x3, ym - j * y3))

        return coords
130778
function(x1, y1, x2, y2, width, height, centerX, centerY, padding) {
      var points = nodeShapes['square'].points;
      
      /*
      return renderer.boxIntersectPolygon(
        x1, y1, x2, y2,
        points, 
      */
    }
780954
def edge_mask(self, angle):
        
        self.point_arrays['point_ind'] = np.arange(self.n_points)
        featureEdges = vtk.vtkFeatureEdges()
        featureEdges.SetInputData(self)
        featureEdges.FeatureEdgesOn()
        featureEdges.BoundaryEdgesOff()
        featureEdges.NonManifoldEdgesOff()
        featureEdges.ManifoldEdgesOff()
        featureEdges.SetFeatureAngle(angle)
        featureEdges.Update()
        edges = _get_output(featureEdges)
        orig_id = vtki.point_scalar(edges, 'point_ind')

        return np.in1d(self.point_arrays['point_ind'], orig_id,
                       assume_unique=True)
548804
def local_maxima(vector,min_distance = 4, brd_mode = "wrap"):
    
    fits = gaussian_filter(numpy.asarray(vector,dtype=numpy.float32),1., mode=brd_mode)
    for ii in range(len(fits)):
        if fits[ii] == fits[ii-1]:
            fits[ii-1] = 0.0
    maxfits     = maximum_filter(fits, size=min_distance, mode=brd_mode)
    maxima_mask = fits == maxfits
    maximum     = numpy.transpose(maxima_mask.nonzero())
    return numpy.asarray(maximum)
130135
function() {
            if (!CSG.Polygon.verticesConvex(this.vertices, this.plane.normal)) {
                CSG.Polygon.verticesConvex(this.vertices, this.plane.normal);
                throw new Error("Not convex!");
            }
        }
903985
function() {
            if (!this.cachedBoundingBox) {
                var minpoint = new CSG.Vector3D(0, 0, 0);
                var maxpoint = new CSG.Vector3D(0, 0, 0);
                var polygons = this.polygons;
                var numpolygons = polygons.length;
                for (var i = 0; i < numpolygons; i++) {
                    var polygon = polygons[i];
                    var bounds = polygon.boundingBox();
                    if (i === 0) {
                        minpoint = bounds[0];
                        maxpoint = bounds[1];
                    } else {
                        minpoint = minpoint.min(bounds[0]);
                        maxpoint = maxpoint.max(bounds[1]);
                    }
                }
                this.cachedBoundingBox = [minpoint, maxpoint];
            }
            return this.cachedBoundingBox;
        }
687639
def slow_minimum_distance2(hull_a, hull_b):
    ''''''
    d2_min = np.iinfo(int).max
    for a in hull_a:
        if within_hull(a, hull_b):
            return 0
    for b in hull_b:
        if within_hull(b, hull_a):
            return 0
    for pt_a in hull_a:
        for pt_b in hull_b:
            d2_min = min(d2_min, np.sum((pt_a - pt_b)**2))
            
    for h1, h2 in ((hull_a, hull_b), (hull_b, hull_a)):
        # Find the distance from a vertex in h1 to an edge in h2
        for pt1 in h1:
            prev_pt2 = h2[-1,:]
            for pt2 in h2:
                if (np.dot(pt2-prev_pt2,pt1-prev_pt2) > 0 and
                    np.dot(prev_pt2-pt2,pt1-pt2) > 0):
                    # points form an acute triangle, so edge is closer
                    # than vertices
                    d2_min = min(d2_min, distance2_to_line(pt1, prev_pt2, pt2))
                prev_pt2 = pt2
    return d2_min
391089
function(polygon) {
            var result = {
                type: null,
                front: null,
                back: null
            };
            // cache in local vars (speedup):
            var planenormal = this.normal;
            var vertices = polygon.vertices;
            var numvertices = vertices.length;
            if (polygon.plane.equals(this)) {
                result.type = 0;
            } else {
                var EPS = CSG.Plane.EPSILON;
                var thisw = this.w;
                var hasfront = false;
                var hasback = false;
                var vertexIsBack = [];
                var MINEPS = -EPS;
                for (var i = 0; i < numvertices; i++) {
                    var t = planenormal.dot(vertices[i].pos) - thisw;
                    var isback = (t < 0);
                    vertexIsBack.push(isback);
                    if (t > EPS) hasfront = true;
                    if (t < MINEPS) hasback = true;
                }
                if ((!hasfront) && (!hasback)) {
                    // all points coplanar
                    var t = planenormal.dot(polygon.plane.normal);
                    result.type = (t >= 0) ? 0 : 1;
                } else if (!hasback) {
                    result.type = 2;
                } else if (!hasfront) {
                    result.type = 3;
                } else {
                    // spanning
                    result.type = 4;
                    var frontvertices = [],
                        backvertices = [];
                    var isback = vertexIsBack[0];
                    for (var vertexindex = 0; vertexindex < numvertices; vertexindex++) {
                        var vertex = vertices[vertexindex];
                        var nextvertexindex = vertexindex + 1;
                        if (nextvertexindex >= numvertices) nextvertexindex = 0;
                        var nextisback = vertexIsBack[nextvertexindex];
                        if (isback == nextisback) {
                            // line segment is on one side of the plane:
                            if (isback) {
                                backvertices.push(vertex);
                            } else {
                                frontvertices.push(vertex);
                            }
                        } else {
                            // line segment intersects plane:
                            var point = vertex.pos;
                            var nextpoint = vertices[nextvertexindex].pos;
                            var intersectionpoint = this.splitLineBetweenPoints(point, nextpoint);
                            var intersectionvertex = new CSG.Vertex(intersectionpoint);
                            if (isback) {
                                backvertices.push(vertex);
                                backvertices.push(intersectionvertex);
                                frontvertices.push(intersectionvertex);
                            } else {
                                frontvertices.push(vertex);
                                frontvertices.push(intersectionvertex);
                                backvertices.push(intersectionvertex);
                            }
                        }
                        isback = nextisback;
                    } // for vertexindex
                    // remove duplicate vertices:
                    var EPS_SQUARED = CSG.Plane.EPSILON * CSG.Plane.EPSILON;
                    if (backvertices.length >= 3) {
                        var prevvertex = backvertices[backvertices.length - 1];
                        for (var vertexindex = 0; vertexindex < backvertices.length; vertexindex++) {
                            var vertex = backvertices[vertexindex];
                            if (vertex.pos.distanceToSquared(prevvertex.pos) < EPS_SQUARED) {
                                backvertices.splice(vertexindex, 1);
                                vertexindex--;
                            }
                            prevvertex = vertex;
                        }
                    }
                    if (frontvertices.length >= 3) {
                        var prevvertex = frontvertices[frontvertices.length - 1];
                        for (var vertexindex = 0; vertexindex < frontvertices.length; vertexindex++) {
                            var vertex = frontvertices[vertexindex];
                            if (vertex.pos.distanceToSquared(prevvertex.pos) < EPS_SQUARED) {
                                frontvertices.splice(vertexindex, 1);
                                vertexindex--;
                            }
                            prevvertex = vertex;
                        }
                    }
                    if (frontvertices.length >= 3) {
                        result.front = new CSG.Polygon(frontvertices, polygon.shared, polygon.plane);
                    }
                    if (backvertices.length >= 3) {
                        result.back = new CSG.Polygon(backvertices, polygon.shared, polygon.plane);
                    }
                }
            }
            return result;
        }
750377
function constrainedMDSLayout(areas, params) {
        params = params || {};
        var restarts = params.restarts || 10;

        // bidirectionally map sets to a rowid  (so we can create a matrix)
        var sets = [],
            setids = {},
            i;
        for (i = 0; i < areas.length; ++i) {
            var area = areas[i];
            if (area.sets.length == 1) {
                setids[area.sets[0]] = sets.length;
                sets.push(area);
            }
        }

        var matrices = getDistanceMatrices(areas, sets, setids),
            distances = matrices.distances,
            constraints = matrices.constraints;

        // keep distances bounded, things get messed up otherwise.
        // TODO: proper preconditioner?
        var norm = norm2(distances.map(norm2)) / distances.length;
        distances = distances.map(function (row) {
            return row.map(function (value) {
                return value / norm;
            });
        });

        var obj = function obj(x, fxprime) {
            return constrainedMDSGradient(x, fxprime, distances, constraints);
        };

        var best, current;
        for (i = 0; i < restarts; ++i) {
            var initial = zeros(distances.length * 2).map(Math.random);

            current = conjugateGradient(obj, initial, params);
            if (!best || current.fx < best.fx) {
                best = current;
            }
        }
        var positions = best.x;

        // translate rows back to (x,y,radius) coordinates
        var circles = {};
        for (i = 0; i < sets.length; ++i) {
            var set = sets[i];
            circles[set.sets[0]] = {
                x: positions[2 * i] * norm,
                y: positions[2 * i + 1] * norm,
                radius: Math.sqrt(set.size / Math.PI)
            };
        }

        if (params.history) {
            for (i = 0; i < params.history.length; ++i) {
                scale(params.history[i].x, norm);
            }
        }
        return circles;
    }