How can I implement Guided Random algorithm walker matrices?

Question

• My Question is related to the implementation of the paper Segmentation by retrieval with guided random walks: Application to left ventricle segmentation in MRI

• After careful reading of the code of the random walkerAlgorithm on scikit-image library website, I tried to implement the Omega, Laplacian, and A matrices from scratch which is mathematically defined as follows:

• Where CI::I mean the neighbourhood, which is 4 connected neighbours (i-1,i+1,j-1,j+1), eij is the edge between two pixels, which is nothing except the difference between their intensities, alpha and beta are arbitrary scalars.
• I have implemented the Laplacian and Omega in the attached code down, but I didn't know how to implement Matrix A, since I don't know how to assign values to slices of the sparse matrices.

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 Function:

def compute_weights(image,mask,alpha, beta, eps=1.e-6):
    # Weight calculation is main difference in multispectral version
    # Original gradient**2 replaced with sum of gradients ** 2
    intra_gradients = np.concatenate([np.diff(image, axis=ax).ravel()
     for ax in [1, 0] ], axis=0) ** 2            # gradient ^2
    # print('intra_gradients shape',intra_gradients.shape)
    # 5-Connected
    inter_gradients = np.concatenate([np.diff(mask, axis=ax).ravel()
    for ax in [1, 0] ], axis=0)**2 
    #----------------------------------------
    # 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

• Building Matrices:

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 )

    # vox = np.concatenate((image[...,np.newaxis], mask[...,np.newaxis]), axis=2)
    # edges_3D = make_graph_edges(vox)
    #================
    # 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

The Answer is:

    #================
    # Matrix A
    #================     
    # Build the sparse linear system  
    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))