使用稀疏矩阵表示聚类二值图像答案

【问题标题】：Clustering binary image using sparse matrix representation使用稀疏矩阵表示聚类二值图像
【发布时间】：2022-01-15 23:01:54
【问题描述】：

我在 StackOverflow 上看到了这个 question，我确实想使用 sparse matrix representation 解决上述问题中的二进制图像聚类问题，例如这个问题：

我知道有更简单有效的聚类方法（KMeans、Mean-shift 等），我知道这个问题可以通过其他解决方案来解决，例如connected components，

但我的目标是使用稀疏矩阵表示方法。

到目前为止我尝试过的是：

读取图像并定义距离函数（以欧式距离为例）：

#!/usr/bin/python3
# -*- coding: utf-8 -*-

import cv2
import numpy as np
from math import sqrt
from scipy.sparse import csr_matrix

original = cv2.imread("km.png", 0)

img = original.copy()

def distance(p1, p2):
    """
    Euclidean Distance
    """
    return sqrt((p1[0]-p2[0])**2 + (p1[1]-p2[1])**2)

根据二值图像的非零点构造稀疏矩阵，为了简单起见我考虑了两个簇，以后可以扩展为k簇：

data = np.nonzero(img)
data = tuple(zip(data[0],data[1]))
data = np.asarray(data).reshape(-1,2)
# print(data.shape)

l = data.shape[0] # 68245

# num clusters
k = 2
cov_mat = csr_matrix((l, l, k), dtype = np.int8).toarray()

第一个循环将质心作为彼此最远的点：

# inefficient loop!
for i in range(l):
    for j in range(l):
        if(i!=j):
            cov_mat[i, j, 0] = distance(data[i], data[j])
# Centroids
# TODO: check if more than two datapoints with same max distance then take first!
ci = cov_mat[...,0].argmax()

计算距每个质心的距离：

# inefficient loop!
# TODO: repeat for k clusters!
for i in range(l):
    for j in range(l):
        if(i!=j):
            cov_mat[i, j, 0] = distance(data[i], data[ci[0]])
            cov_mat[i, j, 1] = distance(data[i], data[ci[1]])

根据最小距离进行聚类：

# Labeling
cov_mat[cov_mat[...,0]>cov_mat[...,1]] = +1
cov_mat[cov_mat[...,0]<cov_mat[...,1]] = -1
# Labeling Centroids
cov_mat[ci[0], ci[0]] = +1
cov_mat[ci[1], ci[1]] = -1

获取集群的索引：

# clusters Indicies
cl1 = cov_mat[...,0].argmax()
cl2 = cov_mat[...,0].argmin()

# TODO: pass back the indices to the image to cluster it.

这种方法很费时间，请问如何提高效率？提前致谢。

【问题讨论】：

是的，不要编写自己的循环。使用同时对整个数据集进行操作的函数（主要是 numpy），而不是对单个点进行操作。每次你有一个循环遍历序列/列表/数组的每个元素时，都要认真思考如何摆脱它。 -- 第一个提示：np.linalg.norm
为什么首先要使用稀疏矩阵？它会让你的循环变得异常缓慢，并防止使用 Numpy 进行任何矢量化。由于 Scipy 它也不是很优化，因此使用 Scipy 函数计算稀疏矩阵不会好很多（至少可能比循环好得多）。在 Python 中处理稀疏矩阵在性能方面是一团糟。像 C/C++ 这样的原生语言可以更有效地做到这一点。
此外，完全构建协方差矩阵是这里的实际问题：它生成一个由2*68245**2 = 9_314_760_050 项组成的矩阵，在 RAM 中应该占用大约 70 GiB。这真的是巨大的，因此创建/填充/读取速度很慢。您的机器上是否有至少 128 GiB 的 RAM？
@JérômeRichard "你为什么要首先使用稀疏矩阵？" 因为它们可以根据 Article 在 GPU 上执行，它使用 close用稀疏矩阵做二进制分割的方法，关于RAM分配的大小，我不想把这个稀疏矩阵转换成密集矩阵，我相信这是稀疏的核心概念，否则，我们将使用numpy数组直接！
好的。对于 RAM 的使用，我只考虑使用当前代码中的 9_314_760_050 项填充 cov_mat 矩阵的循环作为参考（即“低效循环”）。该矩阵是否为稀疏矩阵这一事实并不会改变它在 RAM 中很大的事实（除非距离大多为 0，这不太可能）。稀疏矩阵只有在大多数值为 0 时才有用。在提供的代码中显然不是这种情况，因此它们的效率较低（在 RAM 空间和速度方面）。

标签： python numpy image-processing cluster-analysis sparse-matrix

【解决方案1】：

正如 cmets 中所指出的，在使用 NumPy 数组时，向量化代码应优先于标量代码（即具有显式 for 循环的代码）。此外，正如我所看到的，在这个特定问题中，密集的 NumPy 数组比 Scipy 稀疏数组更好。如果您不受限制使用稀疏数组，则此代码可以：

import numpy as np
from skimage.io import imread
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt

n_clusters = 2
palette = [[255, 0, 0], [0, 255, 0]]

img = imread('https://i.stack.imgur.com/v1NDe.png')
rows, cols = img.nonzero()
coords = np.stack([rows, cols], axis=1)

labels = KMeans(n_clusters=n_clusters, random_state=0).fit_predict(coords)

result = np.zeros(shape=img.shape+(3,), dtype=np.uint8)
for label in range(n_clusters):
    idx = (labels == label)
    result[rows[idx], cols[idx]] = palette[label]

fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(12, 8))
ax0.imshow(img, cmap='gray')
ax0.set_axis_off()
ax0.set_title('Binary image')
ax1.imshow(result)
ax1.set_axis_off()
ax1.set_title('Cluster labels');

【讨论】：