diff --git a/doc/examples/segmentation/plot_rag_boundary.py b/doc/examples/segmentation/plot_rag_boundary.py new file mode 100644 index 00000000..4d909360 --- /dev/null +++ b/doc/examples/segmentation/plot_rag_boundary.py @@ -0,0 +1,27 @@ +""" +========================== +Region Boundary based RAGs +========================== + +This example demonstrates construction of region boundary based RAGs with the +`rag_boundary` function. +""" +from skimage.future import graph +from skimage import data, segmentation, color, filters, io +from skimage.util.colormap import viridis + + +img = data.coffee() +gimg = color.rgb2gray(img) + +labels = segmentation.slic(img, compactness=30, n_segments=400) +edges = filters.sobel(gimg) +edges_rgb = color.gray2rgb(edges) + +g = graph.rag_boundary(labels, edges) + +out = graph.draw_rag(labels, g, edges_rgb, node_color="#999999", + colormap=viridis) + +io.imshow(out) +io.show() diff --git a/skimage/future/graph/__init__.py b/skimage/future/graph/__init__.py index 6b7e0aaa..d685a793 100644 --- a/skimage/future/graph/__init__.py +++ b/skimage/future/graph/__init__.py @@ -1,5 +1,5 @@ from .graph_cut import cut_threshold, cut_normalized -from .rag import rag_mean_color, RAG, draw_rag +from .rag import rag_mean_color, RAG, draw_rag, rag_boundary from .graph_merge import merge_hierarchical ncut = cut_normalized @@ -9,4 +9,5 @@ __all__ = ['rag_mean_color', 'ncut', 'draw_rag', 'merge_hierarchical', + 'rag_boundary', 'RAG'] diff --git a/skimage/future/graph/rag.py b/skimage/future/graph/rag.py index 6a45447b..dae5bbd1 100644 --- a/skimage/future/graph/rag.py +++ b/skimage/future/graph/rag.py @@ -2,6 +2,7 @@ import networkx as nx import numpy as np from numpy.lib.stride_tricks import as_strided from scipy import ndimage as ndi +from scipy import sparse import math from ... import draw, measure, segmentation, util, color try: @@ -11,6 +12,43 @@ except ImportError: pass +def _edge_generator_from_csr(csr_matrix): + """Yield weighted edge triples for use by NetworkX from a CSR matrix. + + This function is a straight rewrite of + `networkx.convert_matrix._csr_gen_triples`. Since that is a private + function, it is safer to include our own here. + + Parameters + ---------- + csr_matrix : scipy.sparse.csr_matrix + The input matrix. An edge (i, j, w) will be yielded if there is a + data value for coordinates (i, j) in the matrix, even if that value + is 0. + + Yields + ------ + i, j, w : (int, int, float) tuples + Each value `w` in the matrix along with its coordinates (i, j). + + Examples + -------- + + >>> dense = np.eye(2, dtype=np.float) + >>> csr = sparse.csr_matrix(dense) + >>> edges = _edge_generator_from_csr(csr) + >>> list(edges) + [(0, 0, 1.0), (1, 1, 1.0)] + """ + nrows = csr_matrix.shape[0] + values = csr_matrix.data + indptr = csr_matrix.indptr + col_indices = csr_matrix.indices + for i in range(nrows): + for j in range(indptr[i], indptr[i + 1]): + yield i, col_indices[j], values[j] + + def min_weight(graph, src, dst, n): """Callback to handle merging nodes by choosing minimum weight. @@ -311,6 +349,70 @@ def rag_mean_color(image, labels, connectivity=2, mode='distance', return graph +def rag_boundary(labels, edge_map, connectivity=2): + """ Comouter RAG based on region boundaries + + Given an image's initial segmentation and its edge map this method + constructs the corresponding Region Adjacency Graph (RAG). Each node in the + RAG represents a set of pixels within the image with the same label in + `labels`. The weight between two adjacent regions is the average value + in `edge_map` along their boundary. + + labels : ndarray + The labelled image. + edge_map : ndarray + This should have the same shape as that of `labels`. For all pixels + along the boundary between 2 adjacent regions, the average value of the + corresponding pixels in `edge_map` is the edge weight between them. + connectivity : int, optional + Pixels with a squared distance less than `connectivity` from each other + are considered adjacent. It can range from 1 to `labels.ndim`. Its + behavior is the same as `connectivity` parameter in + `scipy.ndimage.filters.generate_binary_structure`. + + Examples + -------- + >>> from skimage import data, segmentation, filters, color + >>> from skimage.future import graph + >>> img = data.chelsea() + >>> labels = segmentation.slic(img) + >>> edge_map = filters.sobel(color.rgb2gray(img)) + >>> rag = graph.rag_boundary(labels, edge_map) + + """ + + conn = ndi.generate_binary_structure(labels.ndim, connectivity) + eroded = ndi.grey_erosion(labels, footprint=conn) + dilated = ndi.grey_dilation(labels, footprint=conn) + boundaries0 = (eroded != labels) + boundaries1 = (dilated != labels) + labels_small = np.concatenate((eroded[boundaries0], labels[boundaries1])) + labels_large = np.concatenate((labels[boundaries0], dilated[boundaries1])) + n = np.max(labels_large) + 1 + + # use a dummy broadcast array as data for RAG + ones = as_strided(np.ones((1,), dtype=np.float), shape=labels_small.shape, + strides=(0,)) + count_matrix = sparse.coo_matrix((ones, (labels_small, labels_large)), + dtype=np.int_, shape=(n, n)).tocsr() + data = np.concatenate((edge_map[boundaries0], edge_map[boundaries1])) + + data_coo = sparse.coo_matrix((data, (labels_small, labels_large))) + graph_matrix = data_coo.tocsr() + graph_matrix.data /= count_matrix.data + + rag = RAG() + rag.add_weighted_edges_from(_edge_generator_from_csr(graph_matrix), + weight='weight') + rag.add_weighted_edges_from(_edge_generator_from_csr(count_matrix), + weight='count') + + for n in rag.nodes(): + rag.node[n].update({'labels': [n]}) + + return rag + + def draw_rag(labels, rag, img, border_color=None, node_color='#ffff00', edge_color='#00ff00', colormap=None, thresh=np.inf, desaturate=False, in_place=True): diff --git a/skimage/future/graph/tests/test_rag.py b/skimage/future/graph/tests/test_rag.py index e5882e20..4b2e4cff 100644 --- a/skimage/future/graph/tests/test_rag.py +++ b/skimage/future/graph/tests/test_rag.py @@ -196,3 +196,23 @@ def test_generic_rag_3d(): assert h.has_edge(0, 1) and h.has_edge(0, 3) and not h.has_edge(0, 7) k = graph.RAG(labels, connectivity=3) assert k.has_edge(0, 1) and k.has_edge(1, 2) and k.has_edge(2, 5) + + +def test_rag_boundary(): + labels = np.zeros((16, 16), dtype='uint8') + edge_map = np.zeros_like(labels, dtype=float) + + edge_map[8, :] = 0.5 + edge_map[:, 8] = 1.0 + + labels[:8, :8] = 1 + labels[:8, 8:] = 2 + labels[8:, :8] = 3 + labels[8:, 8:] = 4 + + g = graph.rag_boundary(labels, edge_map, connectivity=1) + assert set(g.nodes()) == set([1, 2, 3, 4]) + assert set(g.edges()) == set([(1, 2), (1, 3), (2, 4), (3, 4)]) + assert g[1][3]['weight'] == 0.25 + assert g[2][4]['weight'] == 0.34375 + assert g[1][3]['count'] == 16