Python scipy.spatial.Delaunay() Examples

def generate_cls_boundary(cls_input,cntr_id_field,boundary_output,cpu_core):
    arcpy.env.parallelProcessingFactor=cpu_core
    arcpy.SetProgressorLabel('Generating Delaunay Triangle...')
    arrays=arcpy.da.FeatureClassToNumPyArray(cls_input,['SHAPE@XY',cntr_id_field])
   
    cid_field_type=[f.type for f in arcpy.Describe(cls_input).fields if f.name==cntr_id_field][0]
    delaunay=Delaunay(arrays['SHAPE@XY']).simplices.copy()
    arcpy.CreateFeatureclass_management('in_memory','boundary_temp','POLYGON',spatial_reference=arcpy.Describe(cls_input).spatialReference)
    fc=r'in_memory\boundary_temp'
    arcpy.AddField_management(fc,cntr_id_field,cid_field_type)
    cursor = arcpy.da.InsertCursor(fc, [cntr_id_field,"SHAPE@"])
    arcpy.SetProgressor("step", "Copying Delaunay Triangle to Temp Layer...",0, delaunay.shape[0], 1)
    for tri in delaunay:
        arcpy.SetProgressorPosition()
        cid=arrays[cntr_id_field][tri[0]]
        if cid == arrays[cntr_id_field][tri[1]] and cid == arrays[cntr_id_field][tri[2]]:
            cursor.insertRow([cid,arcpy.Polygon(arcpy.Array([arcpy.Point(*arrays['SHAPE@XY'][i]) for i in tri]))])
    arcpy.SetProgressor('default','Merging Delaunay Triangle...')
    if '64 bit' in sys.version:
        arcpy.PairwiseDissolve_analysis(fc,boundary_output,cntr_id_field)
    else:
        arcpy.Dissolve_management(fc,boundary_output,cntr_id_field)
    arcpy.Delete_management(fc)
    
    return 

def build_globe(self):
        """Generates the globe as ``vtkPolyData``"""
        # NOTE: https://gitlab.kitware.com/paraview/paraview/issues/19417
        from scipy.spatial import Delaunay
        pos, tex = self.create_sphere()
        pts = self.spherical_to_cartesian(pos[:,0], pos[:,1])
        points = interface.points_to_poly_data(pts).GetPoints()
        texcoords = interface.convert_array(tex, name='Texture Coordinates')
        # Now generate triangles
        cell_connectivity = Delaunay(pos).simplices.astype(int)
        cells = vtk.vtkCellArray()
        cells.SetNumberOfCells(cell_connectivity.shape[0])
        cells.SetCells(cell_connectivity.shape[0], interface.convert_cell_conn(cell_connectivity))
        # Generate output
        output = vtk.vtkPolyData()
        output.SetPoints(points)
        output.GetPointData().SetTCoords(texcoords)
        output.SetPolys(cells)
        return output 

def test_voronoi_name_mapping(xy_of_hex):
    """Test scipy Voronoi names are mapped to landlab-style names."""
    voronoi = Voronoi(xy_of_hex)
    delaunay = Delaunay(xy_of_hex)
    graph = VoronoiDelaunay(xy_of_hex)

    assert np.all(graph.x_of_node == approx(voronoi.points[:, 0]))
    assert np.all(graph.y_of_node == approx(voronoi.points[:, 1]))

    assert np.all(graph.x_of_corner == approx(voronoi.vertices[:, 0]))
    assert np.all(graph.y_of_corner == approx(voronoi.vertices[:, 1]))

    assert np.all(graph.nodes_at_link == voronoi.ridge_points)

    assert tuple(graph.n_corners_at_cell) == tuple(
        len(region) for region in voronoi.regions
    )
    for cell, corners in enumerate(graph.corners_at_cell):
        assert np.all(corners[: graph.n_corners_at_cell[cell]] == voronoi.regions[cell])
        assert np.all(corners[graph.n_corners_at_cell[cell] :] == -1)
    assert np.all(graph.corners_at_face == voronoi.ridge_vertices)
    assert np.all(graph.nodes_at_face == voronoi.ridge_points)
    assert np.all(graph.cell_at_node == voronoi.point_region)

    assert np.all(graph.nodes_at_patch == delaunay.simplices) 

def get_circle(batch_size, masks_size, num_points, device):
    half_dim = masks_size / 2
    half_width = half_dim
    half_height = half_dim

    vert = np.array([[
        half_width + math.floor(math.cos(2 * math.pi / num_points * x) * 10),
        half_height + math.floor(math.sin(2 * math.pi / num_points * x) * 10)]
        for x in range(0, num_points)])
    vert = (vert - half_dim) / half_dim

    tri = Delaunay(vert).simplices.copy()

    vert = torch.Tensor(vert)[None, None, ...].to(device).repeat(batch_size, 1, 1, 1)
    face = torch.Tensor(tri)[None, None, ...].to(device).repeat(batch_size, 1, 1, 1).type(torch.int32)

    vert[:, :, :, 1] = -vert[:, :, :, 1]

    return vert, face 

def get_alpha_shape(pointcloud, alpha):
    pointcloud = np.asarray(pointcloud)

    print(pointcloud,'\n',pointcloud.shape[1])

    assert pointcloud.ndim == 2
    assert pointcloud.shape[1] == 3, "for now, only 3-dimensional analysis is implemented"

    triangulation = spat.Delaunay(pointcloud)

    tetrahedrons = pointcloud[triangulation.simplices]  # remove this copy step, could be fatal
    radii2 = r2_circumsphere_tetrahedron(tetrahedrons[:, 0, :], tetrahedrons[:, 1, :], tetrahedrons[:, 2, :],
                                         tetrahedrons[:, 3, :])
    reduced_triangulation = triangulation.simplices[radii2 < alpha ** 2]
    del radii2, triangulation, tetrahedrons

    outer_triangulation = get_single_faces(reduced_triangulation)

    return outer_triangulation 

def plot_in_hull(p, cloud):
    """
    plot relative to `in_hull` for 2d data
    """
    hull = Delaunay(cloud)

    # plot triangulation
    poly = PolyCollection(hull.points[hull.vertices], facecolors='grey',
                          edgecolors='grey', alpha=0.1)
    plt.clf()
    plt.title('in hull: green, out of hull: red')
    plt.gca().add_collection(poly)
    plt.plot(hull.points[:, 0], hull.points[:, 1], 'o', color='grey',
             alpha=0.2)

    # plot tested points `p` - green are inside hull, red outside
    inside = hull.find_simplex(p) >= 0
    plt.plot(p[inside, 0], p[inside, 1], 'og')
    plt.plot(p[~inside, 0], p[~inside, 1], 'or')
    plt.show() 

def create_random_circle(n, radius, seed=0):
    k = numpy.arange(n)
    boundary_pts = radius * numpy.column_stack(
        [numpy.cos(2 * numpy.pi * k / n), numpy.sin(2 * numpy.pi * k / n)]
    )

    # Compute the number of interior nodes such that all triangles can be somewhat
    # equilateral.
    edge_length = 2 * numpy.pi * radius / n
    domain_area = numpy.pi - n * (
        radius ** 2 / 2 * (edge_length - numpy.sin(edge_length))
    )
    cell_area = numpy.sqrt(3) / 4 * edge_length ** 2
    target_num_cells = domain_area / cell_area
    # Euler:
    # 2 * num_points - num_boundary_edges - 2 = num_cells
    # <=>
    # num_interior_points ~= 0.5 * (num_cells + num_boundary_edges) + 1 - num_boundary_points
    m = int(0.5 * (target_num_cells + n) + 1 - n)

    # Generate random points in circle;
    # <http://mathworld.wolfram.com/DiskPointPicking.html>.
    # Choose the seed such that the fully smoothened mesh has no random boundary points.
    if seed is not None:
        numpy.random.seed(seed)
    r = numpy.random.rand(m)
    alpha = 2 * numpy.pi * numpy.random.rand(m)

    interior_pts = numpy.column_stack(
        [numpy.sqrt(r) * numpy.cos(alpha), numpy.sqrt(r) * numpy.sin(alpha)]
    )

    pts = numpy.concatenate([boundary_pts, interior_pts])

    tri = Delaunay(pts)
    # pts = numpy.column_stack([pts[:, 0], pts[:, 1], numpy.zeros(pts.shape[0])])
    return pts, tri.simplices


# This test iterates over a few meshes that produce weird sitations that did have the
# methods choke. Mostly bugs in GhostedMesh. 

def generate_mask(vertices, num_samples=128):
    """Create a filled convex polygon mask based on the given vertices.

    Parameters
    ----------
    vertices : `iterable`
        ensemble of vertice (x,y) coordinates, in array units
    num_samples : `int`
        number of points in the output array along each dimension

    Returns
    -------
    `numpy.ndarray`
        polygon mask

    """
    vertices = e.asarray(vertices)
    unit = e.arange(num_samples)
    xxyy = e.stack(e.meshgrid(unit, unit), axis=2)

    # use delaunay to fill from the vertices and produce a mask
    triangles = Delaunay(vertices, qhull_options='QJ Qf')
    mask = ~(triangles.find_simplex(xxyy) < 0)
    return mask 

def __create_weights_delaunay_triangulation(self, stimulus):
        """!
        @brief Create weight Denlauny triangulation structure between neurons in line with stimulus.
        
        @param[in] stimulus (list): External stimulus for the chaotic neural network.
        
        """
        
        points = numpy.array(stimulus)
        triangulation = Delaunay(points)
        
        for triangle in triangulation.simplices:
            for index_tri_point1 in range(len(triangle)):
                for index_tri_point2 in range(index_tri_point1 + 1, len(triangle)):
                    index_point1 = triangle[index_tri_point1]
                    index_point2 = triangle[index_tri_point2]
                    
                    weight = self.__calculate_weight(stimulus[index_point1], stimulus[index_point2])
                    
                    self.__weights[index_point1][index_point2] = weight
                    self.__weights[index_point2][index_point1] = weight
                    
                    self.__weights_summary[index_point1] += weight
                    self.__weights_summary[index_point2] += weight 

def __init__(self, size, maxdim=1, complex="scipy"):
        super(LevelSetLayer, self).__init__()
        self.size = size
        self.maxdim = maxdim
        self.fnobj = levelsetdgm()

        # extract width and height
        width, height = size
        if complex == "scipy":
            # initialize complex to use for persistence calculations
            axis_x = np.arange(0, width)
            axis_y = np.arange(0, height)
            grid_axes = np.array(np.meshgrid(axis_x, axis_y))
            grid_axes = np.transpose(grid_axes, (1, 2, 0))

            # creation of a complex for calculations
            tri = Delaunay(grid_axes.reshape([-1, 2]))
            faces = tri.simplices.copy()
            self.complex = self.fnobj.init_filtration(faces)
        elif complex == "freudenthal":
            self.complex = init_freudenthal_2d(width, height)
        else:
            AssertionError("bad complex type") 

def _append_tmp_sources(self):

    from scipy.spatial import cKDTree as kdt
    from scipy.spatial import Delaunay as triag

    sources = row_stack([self.sources]+self.tmp_sources)
    tree = kdt(sources)
    self.sources = sources
    self.tree = tree
    self.tmp_sources = []
    self.tri = triag(
      self.sources,
      incremental=False,
      qhull_options='QJ Qc'
    )
    self.num_sources = len(self.sources)

    return len(sources) 

def from_vertices(cls, data):
        """
        Uses Delauney triangulation to compute triangle simplices for
        each point.
        """
        try:
            from scipy.spatial import Delaunay
        except:
            raise ImportError("Generating triangles from points requires, "
                              "SciPy to be installed.")
        if not isinstance(data, Points):
            data = Points(data)
        if not len(data):
            return cls(([], []))
        tris = Delaunay(data.array([0, 1]))
        return cls((tris.simplices, data)) 

def triangulate_depthmap_points(points2d, depth, pix_thresh=2, depth_thresh=0.2):
    tri = Delaunay(points2d)
    faces_ = tri.simplices.copy()
    faces = []

    for j in range(faces_.shape[0]):
        if np.linalg.norm(points2d[faces_[j, 0], :] - points2d[faces_[j, 1], :]) <= pix_thresh and \
                        np.linalg.norm(points2d[faces_[j, 2], :] - points2d[faces_[j, 1], :]) <= pix_thresh and \
                        np.linalg.norm(points2d[faces_[j, 0], :] - points2d[faces_[j, 2], :]) <= pix_thresh and \
                        np.linalg.norm(depth[faces_[j, 0]] - depth[faces_[j, 1]]) <= depth_thresh and \
                        np.linalg.norm(depth[faces_[j, 2]] - depth[faces_[j, 1]]) <= depth_thresh and \
                        np.linalg.norm(depth[faces_[j, 0]] - depth[faces_[j, 2]]) <= depth_thresh:
            # faces.append(faces_[j, :])
            faces.append(faces_[j, (2, 1, 0)])
    # faces = np.array(faces)
    return faces 

def alpha_shape(points, alpha):
    """
    Compute the alpha shape (concave hull) of a set
    of points.
    @param points: Iterable container of points.
    @param alpha: alpha value to influence the
        gooeyness of the border. Smaller numbers
        don't fall inward as much as larger numbers.
        Too large, and you lose everything!
    """
    if len(points) < 4:
        # When you have a triangle, there is no sense
        # in computing an alpha shape.
        return geometry.MultiPoint(list(points)).convex_hull

    coords = np.array([point.coords[0] for point in points])
    tri = Delaunay(coords)
    triangles = coords[tri.vertices]
    a = ((triangles[:,0,0] - triangles[:,1,0]) ** 2 + (triangles[:,0,1] - triangles[:,1,1]) ** 2) ** 0.5
    b = ((triangles[:,1,0] - triangles[:,2,0]) ** 2 + (triangles[:,1,1] - triangles[:,2,1]) ** 2) ** 0.5
    c = ((triangles[:,2,0] - triangles[:,0,0]) ** 2 + (triangles[:,2,1] - triangles[:,0,1]) ** 2) ** 0.5
    s = ( a + b + c ) / 2.0
    areas = (s*(s-a)*(s-b)*(s-c)) ** 0.5
    circums = a * b * c / (4.0 * areas)
    filtered = triangles[circums < (1.0 / alpha)]
    edge1 = filtered[:,(0,1)]
    edge2 = filtered[:,(1,2)]
    edge3 = filtered[:,(2,0)]
    edge_points = np.unique(np.concatenate((edge1,edge2,edge3)), axis = 0).tolist()
    m = geometry.MultiLineString(edge_points)
    triangles = list(polygonize(m))
    return cascaded_union(triangles), edge_points 

def construct_simplices(self, points, labels, epsilon, distfcn):
        distdict = calculate_distmatrix(points, labels, distfcn)

        delaunayobj = Delaunay(points)

        simplices = []
        for simplexnda in delaunayobj.simplices:
            simplex = tuple(simplexnda)

            detached = [contain_detachededges(face, distdict, epsilon) for face in facesiter(simplex)]

            if True in detached and len(simplex) > 2:
                simplices += [face for face, notkeep in zip(facesiter(simplex), detached)
                              if not notkeep]
            else:
                simplices.append(simplex)
        simplices = map(lambda simplex: tuple(sorted(simplex)), simplices)
        simplices = list(set(simplices))

        allpts = get_allpoints(simplices)
        simplices += [(point,) for point in (set(labels)-allpts)]

        return simplices 

def in_hull(p, hull):
    from scipy.spatial import Delaunay
    if not isinstance(hull, Delaunay):
        hull = Delaunay(hull)
    return hull.find_simplex(p) >= 0 

def morph_by_points (image, sp, dp):
    if sp.shape != dp.shape:
        raise ValueError ('morph_by_points() sp.shape != dp.shape')
    (h,w,c) = image.shape

    result_image = np.zeros(image.shape, dtype = image.dtype)

    for tri in Delaunay(dp).simplices:
        morphTriangle(result_image, image, sp[tri], dp[tri])

    return result_image 

def in_hull(p, hull):
    from scipy.spatial import Delaunay
    if not isinstance(hull, Delaunay):
        hull = Delaunay(hull)
    return hull.find_simplex(p) >= 0 

def init_tri_complex(width, height):
    """
    initialize 2d complex in dumbest possible way
    """
    # initialize complex to use for persistence calculations
    axis_x = np.arange(0, width)
    axis_y = np.arange(0, height)
    grid_axes = np.array(np.meshgrid(axis_x, axis_y))
    grid_axes = np.transpose(grid_axes, (1, 2, 0))

    # creation of a complex for calculations
    tri = Delaunay(grid_axes.reshape([-1, 2]))
    return unique_simplices(tri.simplices, 2) 

def alpha_filtration(x, maxdim=2):
    """
    compute Delaunay triangulation
    fill in simplices as appropriate

    if x is 1-dimensional, defaults to 1D alpha
    inputs:
        x - pointcloud
        maxdim - maximal simplex dimension (default = 2)
    """
    if x.shape[1] == 1:
        x = x.flatten()
        return alpha_filtration_1d(x)
    tri = Delaunay(x)
    simplices = []
    # fill in 0 - cells
    for i in range(len(x)):
        simplices.append(([i], 0.))
    # fill in higher-d - cells
    for s in tri.simplices:
        # edges
        for i, j in itertools.combinations(s, 2):
            simplices.append(([i, j], np.linalg.norm(x[i] - x[j])))
        # higher order simplices
        for dim in range(2, maxdim+1):
            # loop over faces
            for face in itertools.combinations(s, dim+1):
                # get filtration value for face
                vface = 0.
                for i, j in itertools.combinations(face, 2):
                    vface = max(vface, np.linalg.norm(x[i] - x[j]))
                simplices.append((list(face), vface))
    return d.Filtration(simplices) 

def in_hull(p, hull):
    """
    :param p: (N, K) test points
    :param hull: (M, K) M corners of a box
    :return (N) bool
    """
    try:
        if not isinstance(hull, Delaunay):
            hull = Delaunay(hull)
        flag = hull.find_simplex(p) >= 0
    except scipy.spatial.qhull.QhullError:
        print('Warning: not a hull %s' % str(hull))
        flag = np.zeros(p.shape[0], dtype=np.bool)

    return flag 

def in_hull(points, hull):
    """
    Test if a list of coordinates are inside a given convex hull

    Parameters
    ----------
    points : array_like (N x ndims)
        The spatial coordinates of the points to check

    hull : scipy.spatial.ConvexHull object **OR** array_like
        Can be either a convex hull object as returned by
        ``scipy.spatial.ConvexHull`` or simply the coordinates of the points
        that define the convex hull.

    Returns
    -------
    result : 1D-array
        A 1D-array Boolean array of length *N* indicating whether or not the
        given points in ``points`` lies within the provided ``hull``.

    """
    from scipy.spatial import Delaunay, ConvexHull
    if isinstance(hull, ConvexHull):
        hull = hull.points
    hull = Delaunay(hull)
    return hull.find_simplex(points) >= 0 

def cut_finite_fault_surfaces(geo_model, ver:dict, sim:dict):
    """Cut vertices and simplices for finite fault surfaces to finite fault ellipse

    Args:
        geo_model (gempy.core.model.Project): gempy geo_model object
        ver (dict): Dictionary with surfaces as keys and vertices ndarray as values.
        sim (dict): Dictionary with surfaces as keys and simplices ndarray as values.

    Returns:
        ver, sim (dict, dict): Updated vertices and simplices with finite fault
            surfaces cut to ellipses.
    """
    from scipy.spatial import Delaunay
    from copy import copy

    finite_ver = copy(ver)
    finite_sim = copy(sim)

    finite_fault_series = list(geo_model._faults.df[geo_model._faults.df["isFinite"] == True].index)
    finite_fault_surfaces = list(
        geo_model._surfaces.df[geo_model._surfaces.df._stack == finite_fault_series].surface.unique())

    for fault in finite_fault_surfaces:
        U, fpoints_rot, fctr_rot, a, b = get_fault_rotation_objects(geo_model, "Fault 1")
        rpoints = np.dot(ver[fault], U[0])
        # rpoints = np.dot(rpoints, U[-1])
        r = (rpoints[:, 0] - fctr_rot[0]) ** 2 / a ** 2 + (rpoints[:, 1] - fctr_rot[1]) ** 2 / b ** 2

        finite_ver[fault] = finite_ver[fault][r < 1]
        delaunay = Delaunay(finite_ver[fault])
        finite_sim[fault] = delaunay.simplices
        # finite_sim[fault] = finite_sim[fault][np.isin(sim[fault], np.argwhere(r<0.33))]

    return finite_ver, finite_sim 

def in_hull(p, hull):
    from scipy.spatial import Delaunay
    if not isinstance(hull,Delaunay):
        hull = Delaunay(hull)
    return hull.find_simplex(p)>=0 

def get_data(self, element, ranges, style):
        try:
            from scipy.spatial import Delaunay
        except:
            SkipRendering("SciPy not available, cannot plot TriSurface")
        x, y, z = (element.dimension_values(i) for i in range(3))
        points2D = np.vstack([x, y]).T
        tri = Delaunay(points2D)
        simplices = tri.simplices
        return [dict(x=x, y=y, z=z, simplices=simplices)] 

def delaunay_from_points_numpy(points):
    """Computes the delaunay triangulation for a list of points using Numpy.

    Parameters
    ----------
    points : sequence of tuple
        XYZ coordinates of the original points.
    boundary : sequence of tuples
        list of ordered points describing the outer boundary (optional)
    holes : list of sequences of tuples
        list of polygons (ordered points describing internal holes (optional)

    Returns
    -------
    list
        The faces of the triangulation.
        Each face is a triplet of indices referring to the list of point coordinates.

    Example
    -------
    >>>

    """
    xyz = asarray(points)
    d = Delaunay(xyz[:, 0:2])
    return d.simplices 

def triangulatePoints(points):
    points = toNumpy(points)
    triangulation = sp.Delaunay(points)
    triangles = []
    for i in range(len(triangulation.simplices)):
        verts = map(lambda p: shapes.Point(p[0], p[1]),
                    points[triangulation.simplices[i, :]])
        triangle = shapes.Triangle(
            verts[0], verts[1], verts[2])
        triangles.append(triangle)
    return triangles 

def hull_volume(xyz):
    """ Calculate the volume of the convex hull of 3D (X,Y,Z) LMA data.
        xyz is a (N_points, 3) array of point locations in space. """
    assert xyz.shape[1] == 3
        
    tri = Delaunay(xyz[:,0:3])
    vertices = tri.points[tri.vertices]
    
    # This is the volume formula in 
    # https://github.com/scipy/scipy/blob/master/scipy/spatial/tests/test_qhull.py#L106
    # Except the formula needs to be divided by ndim! to get the volume, cf., 
    # http://en.wikipedia.org/wiki/Simplex#Geometric_properties
    # Credit Pauli Virtanen, Oct 14, 2012, scipy-user list
    
    q = vertices[:,:-1,:] - vertices[:,-1,None,:]
    simplex_volumes = (1.0 / factorial(q.shape[-1])) * np.fromiter(
           (np.linalg.det(q[k,:,:]) for k in range(tri.nsimplex)) , dtype=float)
    # print vertices.shape # number of simplices, points per simplex, coords
    # print q.shape
    
    # The simplex volumes have negative values since they are oriented 
    # (think surface normal direction for a triangle
    volume=np.sum(np.abs(simplex_volumes))
    return volume, vertices, simplex_volumes

##############ADDED 01/05/2017 ############### 

def update(self, points, values, reduce_to_active=True):
        # Find all active points
        if reduce_to_active:
            active_simplices = self.active_simplices()
            active_point_idx = np.unique(active_simplices.flatten())
            self.points = self.points[active_point_idx]
            self.values = self.values[active_point_idx]

        self.points = np.concatenate([self.points, points], axis=0)
        self.values = np.concatenate([self.values, values], axis=0)
        self.delaunay = Delaunay(self.points) 

def update(self, points, values, reduce_to_active=True):
        # Find all active points
        if reduce_to_active:
            active_simplices = self.active_simplices()
            active_point_idx = np.unique(active_simplices.flatten())
            self.points = self.points[active_point_idx]
            self.values = self.values[active_point_idx]

        self.points = np.concatenate([self.points, points], axis=0)
        self.values = np.concatenate([self.values, values], axis=0)
        self.delaunay = Delaunay(self.points)