0

我正在使用 Python 包ripser进行持久性同源性。我想利用它来帮助分割 2D 点云。

例如，我正在关注Elizabeth Munch: Python Tutorial on Topological Data Analysis。在这里，我取DoubleAnnulus并增加两者之间的分隔：

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import persim
import ripser
import teaspoon.MakeData.DynSysLib.DynSysLib as DSL
import teaspoon.MakeData.PointCloud as makePtCloud
import teaspoon.TDA.Draw as Draw
from teaspoon.parameter_selection.MsPE import MsPE_tau
from teaspoon.SP.network import ordinal_partition_graph
from teaspoon.SP.network_tools import make_network
from teaspoon.TDA.PHN import PH_network


def DoubleAnnulus(r1=1, R1=2, r2=0.8, R2=1.3, xshift=3):
    P = makePtCloud.Annulus(r=r1, R=R1)
    Q = makePtCloud.Annulus(r=r2, R=R2)

    Q[:, 0] = Q[:, 0] + xshift

    return np.concatenate((P, Q))


def drawTDAtutorial(P, diagrams, R=2):
    fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(20, 5))

    plt.sca(axes[0])
    plt.title("Point Cloud")
    plt.scatter(P[:, 0], P[:, 1])

    plt.sca(axes[1])
    plt.title("0-dim Diagram")
    Draw.drawDgm(diagrams[0])

    plt.sca(axes[2])
    plt.title("1-dim Diagram")
    Draw.drawDgm(diagrams[1])
    plt.axis([0, R, 0, R])

    plt.show()


if __name__ == "__main__":

    P = DoubleAnnulus(r1=1, R1=2, r2=0.5, R2=1.3, xshift=6)
    plt.scatter(*zip(*P))
    plt.show()

    diagrams = ripser.ripser(P)["dgms"]

    drawTDAtutorial(P, diagrams)

双环

如何获取与图上最大的 2 个图点相对应的点云顶点1-dim Diagram？

4

1 回答 1

0

K-means 聚类在这里可以很好地工作（我们的数据相当凸，算法不需要标签）。由于持久性图上的关键特征与线的 y 距离很大y = x，我们可以生成与线的1-dim Diagramy 距离的直方图，y = x并使用高于平均值 3 个标准差的点数作为我们的 -k值。

代码

import pprint

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import persim
import ripser
import teaspoon.MakeData.DynSysLib.DynSysLib as DSL
import teaspoon.MakeData.PointCloud as makePtCloud
import teaspoon.TDA.Draw as Draw
from sklearn.cluster import KMeans
from teaspoon.parameter_selection.MsPE import MsPE_tau
from teaspoon.SP.network import ordinal_partition_graph
from teaspoon.SP.network_tools import make_network
from teaspoon.TDA.PHN import PH_network


def DoubleAnnulus(r1=1, R1=2, r2=0.8, R2=1.3, xshift=3):
    P = makePtCloud.Annulus(r=r1, R=R1)
    Q = makePtCloud.Annulus(r=r2, R=R2)

    Q[:, 0] = Q[:, 0] + xshift

    return np.concatenate((P, Q))


def drawTDAtutorial(P, diagrams, R=2):
    fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(20, 5))

    plt.sca(axes[0])
    plt.title("Point Cloud")
    plt.scatter(P[:, 0], P[:, 1])

    plt.sca(axes[1])
    plt.title("0-dim Diagram")
    Draw.drawDgm(diagrams[0])

    plt.sca(axes[2])
    plt.title("1-dim Diagram")
    Draw.drawDgm(diagrams[1])
    plt.axis([0, R, 0, R])

    plt.show()


if __name__ == "__main__":

    P = DoubleAnnulus(r1=1, R1=2, r2=0.5, R2=1.3, xshift=6)
    # plt.scatter(*zip(*P))
    # plt.show()

    result = ripser.ripser(P)

    diagrams = result["dgms"]
    # drawTDAtutorial(P, diagrams, R=3)

    loop_points = diagrams[1]  # H1 class
    y_heights = loop_points[:, 1] - loop_points[:, 0]
    features = y_heights[y_heights - np.mean(y_heights) >= 3 * np.std(y_heights)]
    k = len(features)

    print(f"Number of Features: {k}")

    kmeans = KMeans(n_clusters=k).fit(P)
    labels = kmeans.labels_

    points1 = P[labels == 0]
    points2 = P[labels == 1]

    fig, ax = plt.subplots()

    ax.scatter(*zip(*points1), color="red")
    ax.scatter(*zip(*points2), color="blue")

    plt.show()

输出

Number of Features: 2

集群的

于 2021-09-02T05:06:37.640 回答