python-3.x - 如何通过对第一列中的行求和并将其他列归零来从另一个构建稀疏矩阵？

Question

• 在我阅读了scikit-image库网站上的random walkerAlgorithm的代码后，我尝试从头开始实现 Omega、Laplacian 和 A 矩阵，其数学定义如下：
• 其中 CI::I 指的是邻域，即 4 个相连的邻域（i-1,i+1,j-1,j+1），eij 是两个像素之间的边缘，除了它们的强度之间的差异，什么都不是， alpha 和 beta 是任意标量。
• 我已经在附加的代码中实现了拉普拉斯算子和欧米茄，但我不知道如何实现矩阵 A，因为我不知道如何为稀疏矩阵的切片赋值。

import numpy as np
import time
from scipy import sparse


def make_graph_edges(image):
    if(len(image.shape)==2):
        n_x, n_y = image.shape
        vertices = np.arange(n_x * n_y ).reshape((n_x, n_y))
        edges_horizontal = np.vstack(( vertices[:, :-1].ravel(), vertices[:, 1:].ravel()))   # X *(Y-1)
        edges_vertical   = np.vstack(( vertices[   :-1].ravel(), vertices[1:   ].ravel()))   #(X-1)* Y  
        edges = np.hstack((edges_horizontal, edges_vertical))
    return edges

• 权重功能：

def compute_weights(image,mask,alpha, beta, eps=1.e-6):
    intra_gradients = np.concatenate([np.diff(image, axis=ax).ravel()
     for ax in [1, 0] ], axis=0) ** 2            # gradient ^2
    # 5-Connected
    inter_gradients = np.concatenate([np.diff(mask, axis=ax).ravel()
    for ax in [1, 0] ], axis=0)**2 
    # inter_gradients = np.concatenate((inter_gradients,(mask-image).ravel()),axis=0)**2  # gradient ^2
    # print('inter_gradients shape',inter_gradients.shape)
    #----------------------------------------
    # 1-Connected
    # inter_gradients = (image - mask)**2
    #----------------------------------------
    # Normalize gradients
    intra_gradients = (intra_gradients - np.amin(intra_gradients))/(np.amax(intra_gradients)- np.amin(intra_gradients))
    inter_gradients = (inter_gradients - np.amin(inter_gradients))/(np.amax(inter_gradients)- np.amin(inter_gradients))
    #------------------------------------------------------
    intra_scale_factor  = -beta  / (10 * image.std())
    intra_weights = np.exp(intra_scale_factor * intra_gradients)
    intra_weights += eps
    #------------------------------------------------------
    inter_scale_factor  = -alpha / (10 * image.std())
    inter_weights = np.exp(inter_scale_factor * inter_gradients)
    inter_weights += eps
    #------------------------------------------------------
    return -intra_weights, inter_weights

• 构建矩阵：

def build_matrices(image, mask, alpha=90, beta=130): 
   edges_2D = make_graph_edges(image)
   intra_weights,inter_weights=compute_weights(image=image,mask=mask,alpha=alpha ,beta=beta, eps=1.e-6 )
    #================
    # Matrix Laplace
    #================    
    # Build the sparse linear system
    pixel_nb  = edges_2D.shape[1]  # N = n_x * (n_y - 1) * +  (n_x - 1) * n_y
    print('Edges Shape: ',edges_2D.shape,'intra-Weights shape: ',intra_weights.shape)
    i_indices = edges_2D.ravel()   # Src - Dest
    print('i',i_indices.shape)
    j_indices = edges_2D[::-1].ravel() # Same list in reverse order ( Dest - Src)
    print('j',j_indices.shape)    
    stacked_intra = np.hstack((intra_weights, intra_weights)) # weights (S-->D, D-->S) are same because graph is undirected
    lap = sparse.coo_matrix((2*stacked_intra, (i_indices, j_indices)), shape=(pixel_nb, pixel_nb))
    lap.setdiag(-2*np.ravel(lap.sum(axis=0)))
    print('Lap',lap.shape)
    Laplace = lap.tocsr()
    #================
    # Matrix Omega
    #================
    # Build the sparse linear system   
    stacked_inter = np.hstack((inter_weights, inter_weights)) # weights (S-->D, D-->S) are same because graph is undirected
    Omeg = sparse.coo_matrix((2*stacked_inter, (i_indices, j_indices)), shape=(pixel_nb, pixel_nb))
    Omeg.setdiag(2*np.ravel((image-mask)**2))
    print('Omeg',Omeg.shape)
    Omega = Omeg.tocsr()
    #================
    # Matrix A
    #================     
    # Build the sparse linear system  
    Mat_A = 0
    return Laplace, Omega, Mat_A

标签： python-3.ximage-processingscipysparse-matrix

Answer 1

答案是：

weights = Omega.copy() 
firstColumn  = weights.sum(axis=1)/2
otherColumns = sparse.csr_matrix((weights.shape[0],weights.shape[1]-1))
Mat_A = sparse.hstack((firstColumn, otherColumns))