1
$\begingroup$
  • 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: weights Lattices Matrices

  • 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
$\endgroup$

1 Answer 1

0
$\begingroup$

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))
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.