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
