random walker segmentation

[Emmanuelle Gouillart]

	Hello,

	a while ago, I contributed to skimage an implementation of the
random walker segmentation algorithm (which has been improved and
extended by many others since then). This algorithm computes a multilabel
segmentation using seeds (already labeled pixels), by determining for an
unlabeled pixel the probability that a seed diffuses to the pixel (with
an anisotropic diffusion coefficient depending on the gradient between
neighboring pixels).

	In the current implementation in skimage, the computation of the
probability map is done by inverting a large sparse linear system
(involving the Laplacian of the graph of pixels). Different methods can
be chosen to solve the linear system: a brute force inversion only
works for tiny images; a conjugate gradient method works well but is
quite slow. If the package pyamg is available, a multigrid method is used
to compute a preconditioner, which speeds up the algorithm -- but it
requires pyamg. Also, the memory cost of the algorithm is important
(linear, I think, though. I haven't yet taken the time to use a memory
profiler but I should).

	Recently, a colleague brought to my attention that the linear
system was just a set of fixed-point equations, that could be solved
iteratively. Indeed, the solution verifies that the probability of a
pixel is the weighted sum (with weights on edges that are a decreasing
function of gradients) of the probabilities of its neighbors. I have
written a quick and dirty implementation (only for 2-D grayscale images
and for 2 labels) of this "local" version, available on
https://github.com/emmanuelle/scikits.image/blob/local_random_walker/skimage/segmentation/local_random_walker.py

	It turns out that this implementation is slower than the
conjugate gradient with multigrid acceleration (typically 2 to three
times slower), but it has several advantages. First, it can be as fast as
the "simple" conjugate gradient (without pyamg's muligrid acceleration),
which is the mode that most users will use (we don't expect users to
install pymag when they are just trying out algorithms). Second, its
memory cost is lower (for example, the weight of an edge is stored only
once, while it appears twice in the Laplacian matrix). Finally, because
the operations only involve neighboring pixels, it is possible that
further speed optimization can be done (using cython... or maybe a GPU
implementation, even if we're not that far yet with skimage).

	So, should we replace the linear algebra implementation with this
simpler local and iterative implementation ? I'd be interested in knowing
about your opinion.

	Cheers,
	Emmanuelle