Hi Josh, thanks for your answer! I'll probably work on the new implementation during the Euroscipy sprint, we'll see whether Cython can decrese the small performance hit :-). Cheers, Emmanuelle On Wed, Jul 31, 2013 at 03:31:19PM -0700, Josh Warner wrote:

Hi Emmanuelle,

I think this would be a good step forward, and would be happy to help! My guess is that this approach could be sped up further in Cython with code we control, making it easier to maintain and expose similar performance to all users.

The memory issue is a real one, particularly for larger (multichannel) datasets. That alone probably justifies the addition, even at the expense of a small performance hit. But we'll see how small we can make the hit ;)

Also, the iterative solution framework might be useful for other algorithms.

Josh

On Tuesday, July 30, 2013 3:50:16 PM UTC-5, Emmanuelle Gouillart wrote:

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