如何通过第一列的行相加，其他列归零的方式从另一列建立稀疏矩阵？

Question

• 当我看完 随机行走算法 关于 scikit-image 我尝试从头开始实现Omega、Laplacian和A矩阵，其数学定义如下。
• 其中CI::I是指邻域，即4个相连的邻域(i-1,i+1,j-1,j+1), eij是两个像素之间的边缘，除了它们之间的强度差之外，什么都没有，α和β是任意的标量。
• 我已经在附件的代码下来实现了Laplacian和Omega，但是我不知道如何实现矩阵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

• weights 函数。

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

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))