Python scipy.spatial.Voronoi() Examples

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 wigner_seitz(self):
        """
        
        Returns
        -------
        TYPE
            Using the Wigner-Seitz Method, this function finds the 1st 
            Brillouin Zone in terms of vertices and faces 
        """
        
        kpoints = []
        for i in range(-1, 2):
            for j in range(-1, 2):
                for k in range(-1, 2):
                    vec = i * self.reciprocal[0] + j * \
                        self.reciprocal[1] + k * self.reciprocal[2]
                    kpoints.append(vec)
        brill = Voronoi(np.array(kpoints))
        faces = []
        for idict in brill.ridge_dict:
            if idict[0] == 13 or idict[1] == 13:
                faces.append(brill.ridge_dict[idict])
                
        verts = brill.vertices
        return np.array(verts), np.array(faces) 

def get_wigner_seitz(recLat):
    kpoints = []
    for i in range(-1, 2):
        for j in range(-1, 2):
            for k in range(-1, 2):
                vec = i * recLat[0] + j * recLat[1] + k * recLat[2]
                kpoints.append(vec)
    brill = Voronoi(np.array(kpoints))
    faces = []
    for idict in brill.ridge_dict:
        if idict[0] == 13 or idict[1] == 13:
            faces.append(brill.ridge_dict[idict])
    verts = brill.vertices
    poly = []
    for ix in range(len(faces)):
        temp = []
        for iy in range(len(faces[ix])):
            temp.append(verts[faces[ix][iy]])
        poly.append(np.array(temp))
    return np.array(poly) 

def plot_cluster_voronoi(data, algo):
    # passing input space to set up voronoi regions.
    points = np.hstack((np.reshape(data[:,0], (len(data[:,0]), 1)), np.reshape(data[:,1], (len(data[:,1]), 1))))
    vor = Voronoi(points)
    # use helper Voronoi
    regions, vertices = voronoi_finite_polygons_2d(vor)
    fig, ax = plt.subplots()
    plot = ax.scatter([], [])
    indice = 0
    for region in regions:
        ax.plot(data[:,0][indice], data[:,1][indice], 'ko')
        polygon = vertices[region]
        # if it isn't gradient based we just color red or blue depending on whether that point uses the machine in question
        color = algo.labels_[indice]
        # we assume only two 
        if color == 0:
            color = 'r'
        else:
            color = 'b'
        ax.fill(*zip(*polygon), alpha=0.4, color=color, label="")
        indice += 1
    ax.axis('equal')
    plt.xlim(vor.min_bound[0] - 0.1, vor.max_bound[0] + 0.1)
    plt.ylim(vor.min_bound[1] - 0.1, vor.max_bound[1] + 0.1) 

def test_ridges(self, name):
        # Check that the ridges computed by Voronoi indeed separate
        # the regions of nearest neighborhood, by comparing the result
        # to KDTree.

        points = DATASETS[name]

        tree = KDTree(points)
        vor = qhull.Voronoi(points)

        for p, v in vor.ridge_dict.items():
            # consider only finite ridges
            if not np.all(np.asarray(v) >= 0):
                continue

            ridge_midpoint = vor.vertices[v].mean(axis=0)
            d = 1e-6 * (points[p[0]] - ridge_midpoint)

            dist, k = tree.query(ridge_midpoint + d, k=1)
            assert_equal(k, p[0])

            dist, k = tree.query(ridge_midpoint - d, k=1)
            assert_equal(k, p[1]) 

def voronoi_from_points_numpy(points):
    """Generate a voronoi diagram from a set of points.

    Parameters
    ----------
    points : list of list of float
        XYZ coordinates of the voronoi sites.

    Returns
    -------

    Examples
    --------
    >>>

    """
    points = asarray(points)
    voronoi = Voronoi(points)
    return voronoi


# ==============================================================================
# Main
# ============================================================================== 

def test_from_3d_voronoi():

    grid = hexa_grid3d(6, 4, 3)
    datasets = from_3d_voronoi(Voronoi(grid))
    assert datasets["vert"].shape[0] == 139
    assert datasets["edge"].shape[0] == 1272
    assert datasets["face"].shape[0] == 141
    assert datasets["cell"].shape[0] == 70
    bulk = Epithelium("bulk", datasets, config.geometry.bulk_spec())
    bulk.reset_index()
    bulk.reset_topo()
    BulkGeometry.update_all(bulk)
    bulk.sanitize()
    assert bulk.validate() 

def planar_sheet_3d(cls, identifier, nx, ny, distx, disty, noise=None):
        """Creates a planar sheet from an hexagonal grid of cells.

        Parameters
        ----------
        identifier : string
        nx, ny : int
          number of cells in the x and y axes
        distx, disty : float,
          the distances in x and y between the cells
        noise : float, default None
          position noise on the hexagonal grid

        Returns
        -------
        flat_sheet: a 2.5D :class:`Sheet` instance
        """

        from scipy.spatial import Voronoi
        from ..config.geometry import flat_sheet
        from ..generation import hexa_grid2d, from_2d_voronoi

        grid = hexa_grid2d(nx, ny, distx, disty, noise)
        datasets = from_2d_voronoi(Voronoi(grid))
        datasets["vert"]["z"] = 0
        datasets["face"]["z"] = 0

        return cls(identifier, datasets, specs=flat_sheet(), coords=["x", "y", "z"]) 

def planar_sheet_2d(cls, identifier, nx, ny, distx, disty, noise=None):
        """Creates a planar sheet from an hexagonal grid of cells.

        Parameters
        ----------
        identifier : string
        nx, ny : int
          number of cells in the x and y axes
        distx, disty : float,
          the distances in x and y between the cells
        noise : float, default None
          position noise on the hexagonal grid

        Returns
        -------
        planar_sheet: a 2D :class:`Sheet` instance
          in the (x, y) plane

        """
        from scipy.spatial import Voronoi
        from ..config.geometry import planar_spec
        from ..generation import hexa_grid2d, from_2d_voronoi

        grid = hexa_grid2d(nx, ny, distx, disty, noise)
        datasets = from_2d_voronoi(Voronoi(grid))
        return cls(identifier, datasets, specs=planar_spec(), coords=["x", "y"]) 

def test_from_2d_voronoi():

    grid = hexa_grid2d(6, 4, 1, 1)
    datasets = from_2d_voronoi(Voronoi(grid))
    assert datasets["vert"].shape[0] == 32
    assert datasets["edge"].shape[0] == 82
    assert datasets["face"].shape[0] == 24 

def test_multisheet():

    base_specs = tyssue.config.geometry.flat_sheet()
    specs = base_specs.copy()
    specs["face"]["layer"] = 0
    specs["vert"]["layer"] = 0
    specs["vert"]["depth"] = 0.0
    specs["edge"]["layer"] = 0
    specs["settings"]["geometry"] = "flat"
    specs["settings"]["interpolate"] = {"function": "multiquadric", "smooth": 0}

    layer_args = [
        (24, 24, 1, 1, 0.4),
        (16, 16, 2, 2, 1),
        (24, 24, 1, 1, 0.4),
        (24, 24, 1, 1, 0.4),
    ]
    dz = 1.0

    layer_datasets = []
    for i, args in enumerate(layer_args):
        centers = hexa_grid2d(*args)
        data = from_2d_voronoi(Voronoi(centers))
        data["vert"]["z"] = i * dz
        layer_datasets.append(data)

    msheet = MultiSheet("more", layer_datasets, specs)
    bbox = [[0, 25], [0, 25]]
    for sheet in msheet:
        edge_out = sheet.cut_out(bbox, coords=["x", "y"])
        sheet.remove(edge_out)

    assert len(msheet) == 4
    datasets = msheet.concat_datasets()
    assert np.all(np.isfinite(datasets["vert"][["x", "y", "z"]]))

    msheet.update_interpolants()
    for interp in msheet.interpolants:
        assert np.isfinite(interp(10, 10)) 

def voronoi_plot(clusters: list, columns: list):
	check_types([
		("clusters", clusters, [list], False),
		("columns", columns, [list], False)])
	from scipy.spatial import voronoi_plot_2d, Voronoi
	v = Voronoi(clusters)
	voronoi_plot_2d(v, show_vertices = 0)
	plt.xlabel(columns[0])
	plt.ylabel(columns[1])
	plt.show() 

def Voronoi_tessalate(atoms):

	base_coords = [a.position for a in atoms]
	repeat_unit_cell = atoms.get_cell().T
	inv_ruc = np.linalg.inv(repeat_unit_cell)
	basis = [np.array([1,0,0]), np.array([0,1,0]), np.array([0,0,1])]
	mesh = []

	for coord in base_coords:

		mesh.append(coord)
		
		fcoord = np.dot(inv_ruc, coord)
		zero_threshold_indices = fcoord < 1e-6
		fcoord[zero_threshold_indices] = 0.0
		one_threshold_indices = abs(fcoord - 1.0) < 1e-6
		fcoord[one_threshold_indices] = 0.0

		if np.all(fcoord):
			trans_vecs = [-1 * b for b in basis] + basis
		else:
			trans_vecs = [basis[dim] for dim in (0,1,2) if fcoord[dim] == 0.0]
			combs = list(itertools.combinations(trans_vecs, 2)) + list(itertools.combinations(trans_vecs, 3))
			for comb in combs:
				compound = np.array([0.0,0.0,0.0])
				for vec in comb:
					compound += vec
				trans_vecs.append(compound)

		for vec in trans_vecs:
			trans_coord = [np.round(i, 6) for i in np.dot(repeat_unit_cell, fcoord + vec)]
			mesh.append(trans_coord)

	mesh = np.asarray(mesh)
	vor = Voronoi(mesh)

	return vor, mesh 

def test_hexagrid3d_noise():
    np.random.seed(1)
    grid = hexa_grid3d(6, 4, 3, noise=0.1)
    datasets = from_3d_voronoi(Voronoi(grid))
    assert datasets["vert"].shape[0] == 318
    assert datasets["edge"].shape[0] == 3300
    assert datasets["face"].shape[0] == 335
    assert datasets["cell"].shape[0] == 72 

def _get_Voronoi_edges(vor):
    r"""
    Given a Voronoi object as produced by the scipy.spatial.Voronoi class,
    this function calculates the start and end points of eeach edge in the
    Voronoi diagram, in terms of the vertex indices used by the received
    Voronoi object.

    Parameters
    ----------
    vor : scipy.spatial.Voronoi object

    Returns
    -------
    A 2-by-N array of vertex indices, indicating the start and end points of
    each vertex in the Voronoi diagram.  These vertex indices can be used to
    index straight into the ``vor.vertices`` array to get spatial positions.
    """
    edges = [[], []]
    for facet in vor.ridge_vertices:
        # Create a closed cycle of vertices that define the facet
        edges[0].extend(facet[:-1]+[facet[-1]])
        edges[1].extend(facet[1:]+[facet[0]])
    edges = np.vstack(edges).T  # Convert to scipy-friendly format
    mask = np.any(edges == -1, axis=1)  # Identify edges at infinity
    edges = edges[~mask]  # Remove edges at infinity
    edges = np.sort(edges, axis=1)  # Move all points to upper triangle
    # Remove duplicate pairs
    edges = edges[:, 0] + 1j*edges[:, 1]  # Convert to imaginary
    edges = np.unique(edges)  # Remove duplicates
    edges = np.vstack((np.real(edges), np.imag(edges))).T  # Back to real
    edges = np.array(edges, dtype=int)
    return edges 

def get_points_to_poly_assignments(poly_to_pt_assignments):
    """
    Reverse of poly to points assignments: Returns a list of size N, which is the number of unique points in
    `poly_to_pt_assignments`. Each list element is an index into the list of Voronoi regions.
    """
    pt_poly = [(i_pt, i_vor)
               for i_vor, pt_indices in enumerate(poly_to_pt_assignments)
               for i_pt in pt_indices]

    return [i_vor for _, i_vor in sorted(pt_poly, key=lambda x: x[0])] 

def polygon_shapes_from_voronoi_lines(poly_lines, geo_shape=None, shapes_from_diff_with_min_area=None):
    """
    Form shapely Polygons objects from a list of shapely LineString objects in `poly_lines` by using
    [`polygonize`](http://toblerity.org/shapely/manual.html#shapely.ops.polygonize). If `geo_shape` is not None, then
    the intersection between any generated polygon and `geo_shape` is taken in case they overlap (i.e. the Voronoi
    regions at the border are "cut" to the `geo_shape` polygon that represents the geographic area holding the
    Voronoi regions). Setting `shapes_from_diff_with_min_area` fixes rare errors where the Voronoi shapes do not fully
    cover `geo_shape`. Set this to a small number that indicates the minimum valid area of "fill up" Voronoi region
    shapes.
    Returns a list of shapely Polygons objects.
    """

    # generate shapely Polygon objects from the LineStrings of the Voronoi shapes in `poly_lines`
    poly_shapes = []
    for p in polygonize(poly_lines):
        if geo_shape is not None and not geo_shape.contains(p):    # if `geo_shape` contains polygon `p`,
            p = p.intersection(geo_shape)                          # intersect it with `geo_shape` (i.e. "cut" it)

        if not p.is_empty:
            poly_shapes.append(p)

    if geo_shape is not None and shapes_from_diff_with_min_area is not None:
        # fix rare cases where the generated polygons of the Voronoi regions don't fully cover `geo_shape`
        vor_polys_union = cascaded_union(poly_shapes)   # union of Voronoi regions
        diff = np.array(geo_shape.difference(vor_polys_union), dtype=object)    # "gaps"
        diff_areas = np.array([p.area for p in diff])    # areas of "gaps"
        # use only those "gaps" bigger than `shapes_from_diff_with_min_area` because very tiny areas are generated
        # at the borders due to floating point errors
        poly_shapes.extend(diff[diff_areas >= shapes_from_diff_with_min_area])

    return poly_shapes 

def test_to_nd():
    grid = hexa_grid2d(6, 4, 3, 3)
    datasets = from_2d_voronoi(Voronoi(grid))
    sheet = Sheet("test", datasets)
    result = utils._to_3d(sheet.face_df["x"])
    assert (sheet.face_df[['x', 'y', 'z']] * result).shape[1] == 3 

def test_single_cell():
    grid = hexa_grid3d(6, 4, 3)
    datasets = from_3d_voronoi(Voronoi(grid))
    sheet = Sheet("test", datasets)
    eptm = utils.single_cell(sheet, 1)
    assert len(eptm.edge_df) == len(sheet.edge_df[sheet.edge_df["cell"] == 1]) 

def test_multisheet():

    base_specs = tyssue.config.geometry.flat_sheet()
    specs = base_specs.copy()
    specs["face"]["layer"] = 0
    specs["vert"]["layer"] = 0
    specs["vert"]["depth"] = 0.0
    specs["edge"]["layer"] = 0
    specs["settings"]["geometry"] = "flat"
    specs["settings"]["interpolate"] = {"function": "multiquadric", "smooth": 0}
    layer_args = [
        (24, 24, 1, 1, 0.4),
        (16, 16, 2, 2, 1),
        (24, 24, 1, 1, 0.4),
        (24, 24, 1, 1, 0.4),
    ]
    dz = 1.0

    layer_datasets = []
    for i, args in enumerate(layer_args):
        centers = hexa_grid2d(*args)
        data = from_2d_voronoi(Voronoi(centers))
        data["vert"]["z"] = i * dz
        layer_datasets.append(data)

    msheet = MultiSheet("more", layer_datasets, specs)
    bbox = [[0, 25], [0, 25]]
    for sheet in msheet:
        edge_out = sheet.cut_out(bbox, coords=["x", "y"])
        sheet.remove(edge_out)

    MultiSheetGeometry.update_all(msheet)
    assert np.all(np.isfinite(msheet[0].face_df["area"])) 

def test_plot24(self):
        np.random.seed(seed)
        X = [[np.random.random(), np.random.random()] for i in range(25)]
        vor = Voronoi(X)
        voronoi_plot_2d(vor)
        self.newsave()
        plt.close() 

def plot_2d_clusters(X, labels, centers):
    """
    Given an observation array, a label vector, and the location of the centers
    plot the clusters
    """

    clabels = set(labels)
    K = len(clabels)

    if len(centers) != K:
        raise ValueError("Expecting the number of unique labels and centres to"
                         " be the same!")

    # Plot the true clusters
    figure(figsize=(10, 10))
    ax = gca()

    vor = Voronoi(centers)

    voronoi_plot_2d(vor, ax)

    colors = cm.hsv(np.arange(K)/float(K))
    for k, col in enumerate(colors):
        my_members = labels == k
        scatter(X[my_members, 0], X[my_members, 1], c=col, marker='o', s=20)

    for k, col in enumerate(colors):
        cluster_center = centers[k]
        scatter(cluster_center[0], cluster_center[1], c=col, marker='o', s=200)

    axis('tight')
    axis('equal')
    title('Clusters') 

def relaxLloyd(pts,strength):
    for i in range(strength):
        vor = Voronoi(pts)
        newpts = []
        for idx in range(len(vor.points)):
            pt = vor.points[idx,:]
            region = vor.regions[vor.point_region[idx]]
            if -1 in region:
                newpts.append(pt)
            else:
                vxs = np.asarray([vor.vertices[i,:] for i in region])
                newpt = centroidnp(vxs)
                newpts.append(newpt)
        pts = np.array(newpts) 

def _voronoi(structure, limit=5, tol=1e-2):
    lps = structure.lat_params
    limits = np.array([ int(np.ceil(limit/lps[i])) for i in range(3)])

    # Create a new cell large enough to find neighbors at least `limit` away.
    new_struct = structure.transform(limits+1, in_place=False)
    nns = defaultdict(list)

    # look for neighbors for each atom in the structure
    for i, atom in enumerate(structure.atoms):
        atom.neighbors = []

        # center on the atom so that its voronoi cell is fully described
        new_struct.recenter(atom, middle=True)
        tess = Voronoi(new_struct.cartesian_coords)
        for j, r in enumerate(tess.ridge_points):
            if not i in r: # only ridges involving the specified atom matter
                continue
            inds = tess.ridge_vertices[j]
            verts = tess.vertices[inds]

            # check that the area of the facet is large enough
            if _get_facet_area(verts) < tol:
                continue

            # Get the indices of all points which this atom shares a ridge with
            ind = [k for k in tess.ridge_points[j] if k != i ][0]
            # map the atom index back into the original cell. 
            atom.neighbors.append(atom)
        ## nns[atom] = atom.neighbors
        ## `qmpy.Atom` objects that are not saved are not hashable, and hence
        ## cannot be used as dictionary keys. Unit tests will fail.
        nns[i] = atom.neighbors
    return nns 

def get_voronoi_polygons(self, image: np.ndarray):
        points = self.get_keypoints(image, include_corners=False)
        voronoi = spatial.Voronoi(points)
        regions, vertices = self.finitize_voronoi(voronoi)
        self.constrain_points(image.shape[:2], vertices)
        max_polygon_points = len(max(regions, key=lambda x: len(x)))
        polygons = np.full((len(regions), max_polygon_points, 2), -1)
        for i, region in enumerate(regions):
            polygon_points = vertices[region]
            polygons[i][:polygon_points.shape[0]] = polygon_points
        return polygons 

def _compare_qvoronoi(self, points, output, **kw):
        """Compare to output from 'qvoronoi o Fv < data' to Voronoi()"""

        # Parse output
        output = [list(map(float, x.split())) for x in output.strip().splitlines()]
        nvertex = int(output[1][0])
        vertices = list(map(tuple, output[3:2+nvertex]))  # exclude inf
        nregion = int(output[1][1])
        regions = [[int(y)-1 for y in x[1:]]
                   for x in output[2+nvertex:2+nvertex+nregion]]
        nridge = int(output[2+nvertex+nregion][0])
        ridge_points = [[int(y) for y in x[1:3]]
                        for x in output[3+nvertex+nregion:]]
        ridge_vertices = [[int(y)-1 for y in x[3:]]
                          for x in output[3+nvertex+nregion:]]

        # Compare results
        vor = qhull.Voronoi(points, **kw)

        def sorttuple(x):
            return tuple(sorted(x))

        assert_allclose(vor.vertices, vertices)
        assert_equal(set(map(tuple, vor.regions)),
                     set(map(tuple, regions)))

        p1 = list(zip(list(map(sorttuple, ridge_points)), list(map(sorttuple, ridge_vertices))))
        p2 = list(zip(list(map(sorttuple, vor.ridge_points.tolist())),
                 list(map(sorttuple, vor.ridge_vertices))))
        p1.sort()
        p2.sort()

        assert_equal(p1, p2) 

def test_simple(self):
        # Simple case with known Voronoi diagram
        points = [(0, 0), (0, 1), (0, 2),
                  (1, 0), (1, 1), (1, 2),
                  (2, 0), (2, 1), (2, 2)]

        # qhull v o Fv Qbb Qc Qz < dat
        output = """
        2
        5 10 1
        -10.101 -10.101
           0.5    0.5
           1.5    0.5
           0.5    1.5
           1.5    1.5
        2 0 1
        3 3 0 1
        2 0 3
        3 2 0 1
        4 4 3 1 2
        3 4 0 3
        2 0 2
        3 4 0 2
        2 0 4
        0
        12
        4 0 3 0 1
        4 0 1 0 1
        4 1 4 1 3
        4 1 2 0 3
        4 2 5 0 3
        4 3 4 1 2
        4 3 6 0 2
        4 4 5 3 4
        4 4 7 2 4
        4 5 8 0 4
        4 6 7 0 2
        4 7 8 0 4
        """
        self._compare_qvoronoi(points, output) 

def test_masked_array_fails(self):
        masked_array = np.ma.masked_all(1)
        assert_raises(ValueError, qhull.Voronoi, masked_array) 

def test_issue_8051(self):
        points = np.array([[0, 0], [0, 1], [0, 2], [1, 0], [1, 1], [1, 2],[2, 0], [2, 1], [2, 2]])
        Voronoi(points) 

def _get_voronoi_vertices_and_ridges(self):
        borders = self._get_densified_borders()

        voronoi_diagram = Voronoi(borders)
        vertices = voronoi_diagram.vertices
        ridges = voronoi_diagram.ridge_vertices

        return vertices, ridges