ENH: addressed Tony's comments (renaming of functions, doc)

This commit is contained in:
emmanuelle
2012-01-21 11:43:53 +01:00
parent b9d6090c5e
commit 528187c2b0
3 changed files with 35 additions and 20 deletions
@@ -3,20 +3,21 @@
Random walker segmentation
==========================
The random walker algorithm (*Random walks for image segmentation*, Leo Grady,
IEEE Trans Pattern Anal Mach Intell. 2006 Nov;28(11):1768-83) determines the
segmentation of an image from a set of markers labeling several phases (2 or
more). An anisotropic diffusion equation is solved with tracers initiated
at the markers' position. The local diffusivity coefficient is greater if
neighboring pixels have similar values, so that diffusion is difficult across
high gradients. The label of each unknown pixel is attributed to the label
of the known marker that has the highest probability to be reached first
during this diffusion process.
The random walker algorithm [1]_ determines the segmentation of an image from
a set of markers labeling several phases (2 or more). An anisotropic diffusion
equation is solved with tracers initiated at the markers' position. The local
diffusivity coefficient is greater if neighboring pixels have similar values,
so that diffusion is difficult across high gradients. The label of each unknown
pixel is attributed to the label of the known marker that has the highest
probability to be reached first during this diffusion process.
In this example, two phases are clearly visible, but the data are too
noisy to perform the segmentation from the histogram only. We determine
markers of the two phases from the extreme tails of the histogram of gray
values, and use the random walker for the segmentation.
.. [1] *Random walks for image segmentation*, Leo Grady, IEEE Trans. Pattern
Anal. Mach. Intell. 2006 Nov; 28(11):1768-83
"""
print __doc__
@@ -31,18 +31,24 @@ from scipy.sparse.linalg import cg
#-----------Laplacian--------------------
def _make_edges_3d(n_x, n_y, n_z):
def _make_graph_edges_3d(n_x, n_y, n_z):
"""
Returns a list of edges for a 3D image.
Parameters
===========
----------
n_x: integer
The size of the grid in the x direction.
n_y: integer
The size of the grid in the y direction
n_z: integer
The size of the grid in the z direction
Returns
-------
edges : ndarray of shape (2, n_x * n_y * (nz - 1) +
n_x * (n_y - 1) * nz + (n_x - 1) * n_y * nz)
Graph edges with each column describing a node-id pair.
"""
vertices = np.arange(n_x * n_y * n_z).reshape((n_x, n_y, n_z))
edges_deep = np.vstack((vertices[:, :, :-1].ravel(),
@@ -55,7 +61,6 @@ def _make_edges_3d(n_x, n_y, n_z):
def _compute_weights_3d(edges, data, beta=130, eps=1.e-6):
l_x, l_y, l_z = data.shape
gradients = _compute_gradients_3d(data) ** 2
beta /= 10 * data.std()
gradients *= beta
@@ -65,7 +70,6 @@ def _compute_weights_3d(edges, data, beta=130, eps=1.e-6):
def _compute_gradients_3d(data):
l_x, l_y, l_z = data.shape
gr_deep = np.abs(data[:, :, :-1] - data[:, :, 1:]).ravel()
gr_right = np.abs(data[:, :-1] - data[:, 1:]).ravel()
gr_down = np.abs(data[:-1] - data[1:]).ravel()
@@ -101,7 +105,6 @@ def _buildAB(lap_sparse, labels):
Build the matrix A and rhs B of the linear system to solve.
A and B are two block of the laplacian of the image graph.
"""
l_x, l_y, l_z = labels.shape
labels = labels[labels >= 0]
indices = np.arange(labels.size)
unlabeled_indices = indices[labels == 0]
@@ -119,25 +122,31 @@ def _buildAB(lap_sparse, labels):
return lap_sparse, rhs
def _trim_edges_weights(edges, weights, mask):
def _mask_edges_weights(edges, weights, mask):
"""
Remove edges of the graph connected to masked nodes, as well as
corresponding weights of the edges.
"""
mask0 = np.hstack((mask[:, :, :-1].ravel(), mask[:, :-1].ravel(),
mask[:-1].ravel()))
mask1 = np.hstack((mask[:, :, 1:].ravel(), mask[:, 1:].ravel(),
mask[1:].ravel()))
ind_mask = np.logical_and(mask0, mask1)
edges, weights = edges[:, ind_mask], weights[ind_mask]
maxval = edges.max()
order = np.searchsorted(np.unique(edges.ravel()), np.arange(maxval + 1))
max_node_index = edges.max()
# Reassign edges labels to 0, 1, ... edges_number - 1
order = np.searchsorted(np.unique(edges.ravel()),
np.arange(max_node_index + 1))
edges = order[edges]
return edges, weights
def _build_laplacian(data, mask=None, beta=50):
l_x, l_y, l_z = data.shape
edges = _make_edges_3d(l_x, l_y, l_z)
edges = _make_graph_edges_3d(l_x, l_y, l_z)
weights = _compute_weights_3d(edges, data, beta=beta, eps=1.e-10)
if mask is not None:
edges, weights = _trim_edges_weights(edges, weights, mask)
edges, weights = _mask_edges_weights(edges, weights, mask)
lap = _make_laplacian_sparse(edges, weights)
del edges, weights
return lap
@@ -188,7 +197,8 @@ def random_walker(data, labels, beta=130, mode='bf', tol=1.e-3, copy=True):
(http://code.google.com/p/pyamg/) is installed. For images of
size > 512x512, this is the recommended (fastest) mode.
tol : tolerance to achieve when solving the linear system, in
tol : float
tolerance to achieve when solving the linear system, in
cg' and 'cg_mg' modes.
copy : bool
@@ -105,3 +105,7 @@ def test_3d_inactive():
labels = random_walker(data, labels, mode='cg')
assert (labels.reshape(data.shape)[13:17, 13:17, 13:17] == 2).all()
return data, labels, old_labels, after_labels
if __name__ == '__main__':
from numpy import testing
testing.run_module_suite()