Mailman 3 August 2013 - scikit-image

random walker segmentation
by Emmanuelle Gouillart 05 Aug '13

05 Aug '13

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/skimag… 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

3 3

Re: Fail to import peak_local_max
by Guillaume Gay 05 Aug '13

05 Aug '13

That did the trick, thanks Juan! G. Le 05/08/2013 12:14, Juan Nunez-Iglesias a écrit : > This sort of stuff keeps popping up (see latest issues) because a lot > of functions have moved to use memory views. If you do a `make clean` > and then try again, it should work. > > > On Mon, Aug 5, 2013 at 7:51 PM, Guillaume Gay > <guillaume(a)mitotic-machine.org <mailto:guillaume@mitotic-machine.org>> > wrote: > > Hi every one, > > I encountered a bug when importing |peak_local_max| (with the > latest github version). The traceback goes up to feature.template. > Running |test_template.py| gives this: > > |python /usr/local/src/scikit-image/skimage/feature/tests/test_template.py > Traceback (most recent call last): > File"/usr/local/src/scikit-image/skimage/feature/tests/test_template.py", line5,in <module> > from skimage.featureimport match_template, peak_local_max > File"/usr/local/lib/python2.7/dist-packages/scikit_image-0.9dev-py2.7-linux-i686.egg/skimage/feature/__init__.py", line9,in <module> > from .templateimport match_template > File"/usr/local/lib/python2.7/dist-packages/scikit_image-0.9dev-py2.7-linux-i686.egg/skimage/feature/template.py", line4,in <module> > from .import _template > File"_template.pyx", line1,in init skimage.feature._template (skimage/feature/_template.c:4042) > TypeError: C function skimage._shared.transform.integrate has wrong signature (expected float (PyArrayObject *, Py_ssize_t, Py_ssize_t, Py_ssize_t, Py_ssize_t), got float (__Pyx_memviewslice, Py_ssize_t, Py_ssize_t, Py_ssize_t, Py_ssize_t))| > > Any hint how to fix this? > > Cheers, > > Guillaume > > -- > You received this message because you are subscribed to the Google > Groups "scikit-image" group. > To unsubscribe from this group and stop receiving emails from it, > send an email to scikit-image+unsubscribe(a)googlegroups.com > <mailto:scikit-image%2Bunsubscribe@googlegroups.com>. > For more options, visit https://groups.google.com/groups/opt_out. > > > > -- > You received this message because you are subscribed to the Google > Groups "scikit-image" group. > To unsubscribe from this group and stop receiving emails from it, send > an email to scikit-image+unsubscribe(a)googlegroups.com. > For more options, visit https://groups.google.com/groups/opt_out. > >

1 0

Fail to import peak_local_max
by Guillaume Gay 05 Aug '13

05 Aug '13

Hi every one, I encountered a bug when importing |peak_local_max| (with the latest github version). The traceback goes up to feature.template. Running |test_template.py| gives this: |python /usr/local/src/scikit-image/skimage/feature/tests/test_template.py Traceback (most recent call last): File"/usr/local/src/scikit-image/skimage/feature/tests/test_template.py", line5,in <module> from skimage.featureimport match_template, peak_local_max File"/usr/local/lib/python2.7/dist-packages/scikit_image-0.9dev-py2.7-linux-i686.egg/skimage/feature/__init__.py", line9,in <module> from .templateimport match_template File"/usr/local/lib/python2.7/dist-packages/scikit_image-0.9dev-py2.7-linux-i686.egg/skimage/feature/template.py", line4,in <module> from .import _template File"_template.pyx", line1,in init skimage.feature._template (skimage/feature/_template.c:4042) TypeError: C function skimage._shared.transform.integrate has wrong signature (expected float (PyArrayObject *, Py_ssize_t, Py_ssize_t, Py_ssize_t, Py_ssize_t), got float (__Pyx_memviewslice, Py_ssize_t, Py_ssize_t, Py_ssize_t, Py_ssize_t))| Any hint how to fix this? Cheers, Guillaume

2 1

Re: cython -a
by Johannes Schönberger 05 Aug '13

05 Aug '13

Hi Juan, > 1. about optimising Cython code in skimage: how do I compile with -a? Because of local imports and so on, calling it on the source file directly is not an option, and I haven't figured out an obvious place to make it spit out the html otherwise… cython -a does not depend on any imports. At least on my system it works for all files, regardless of imports… ? > 2. about running tests: how do I easily figure out which tests were skipped when running "make test" and why they were skipped? Since make test just invokes nosetests, you could simply use nosetests --no-skip.

2 1

cython -a
by Juan Nunez-Iglesias 05 Aug '13

05 Aug '13

Hey skimagers, 2 quick questions: 1. about optimising Cython code in skimage: how do I compile with -a? Because of local imports and so on, calling it on the source file directly is not an option, and I haven't figured out an obvious place to make it spit out the html otherwise... 2. about running tests: how do I easily figure out which tests were skipped when running "make test" and why they were skipped? I'll also take this chance to thank Stefan for pushing me to modify SLIC to use memory views. I dragged my feet but in the end the code is way better and it's making my life much easier with the spacing modifications I'm working on now...! Don't you hate the unreasonable effectiveness of peer review? ;) Thanks, Juan.

1 0

Re: Blog updated
by Chintak Sheth 04 Aug '13

04 Aug '13

Aah, great to find you on the scikit-image mailing list! I can certainly learn a lot from you. =) Coming to the aspect you pointed out, I was thinking if we could probably get rid of `valid_template` all together from this evaluation and introduce another einsum product of `ssd` with `valid_mask`? For example, this is what we are doing, `c*a**2 + c*b**2 - 2c*a*b`, which really is equivalent to `c*(a-b)**2`. However, this does add another call to einsum. And I had performed tests with about 500 pixels which amount to 1500 einsum calls. The bottleneck turns out to be einsum calls. (PR<https://github.com/scikit-image/scikit-image/pull/670>) This is why I refrained from adding another einsum call, since I'd have 2000 calls then. However, your tests indicate that calls with 2 parameters do run considerably faster, so may be I'll go ahead and make this change ? What are your thoughts on this ? Chintak

1 0

Re: Blog updated
by Jaime Fernández del Río 01 Aug '13

01 Aug '13

Hi Chintak, You have processed the StackOverflow answers impresively well: great blog post! Just a quick note on performance of np.einsum. I have found that it performs much better when handed only two parameters. So you may want to benchmark whether applying the mask to the template before the call to np.einsum makes your code run faster. I don't think there is a way out of this 3 parameter call: ssd += np.einsum('ijkl, ijkl, kl->ij', y, y, valid_mask) But there is a good chance that: ssd = np.einsum('ijkl, kl, kl->ij', y, template, valid_mask, dtype=np.float) runs noticeably faster as: ssd = np.einsum('ijkl, kl->ij', y, template*valid_mask, dtype=np.float) A quick test on my system, with a 1000x1000 image and a 9x9 template and mask, all of floats, show it's 25% faster. And this is where about half of your processing time is being spent, so that little change would give you a 10% performance boost for free in this particular case. You may want to test a wider variety of parameter sizes, to see if the improvement holds. The third call to np.einsum has a negligible impact on overall performance, but if you store the value of template*valid_mask, it also runs faster with a two parameter call, i.e. as: ssd += np.einsum('ij, ij', template, cached_template_times_valid_mask) Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.

1 0

PEP8 Update
by Johannes Schönberger 01 Aug '13

01 Aug '13

For anyone interested: http://hg.python.org/peps/rev/fb24c80e9afb Johannes Schönberger

1 0

Re: Blog updated
by Chintak Sheth 01 Aug '13

01 Aug '13

Hi Stefan On Wed, Jul 31, 2013 at 1:53 AM, Stéfan van der Walt <stefan(a)sun.ac.za>wrote: > > You can often write einsum as some operation + a reduction over one or > more axes. Can this operation be expressed in both ways? If so, > which do you find to be more readable and understandable? > Yes we can. Sum the view and the template and then reduction over 2nd and 3rd axes. However, a simple addition would result in an actual 4D array since it would no longer be able to store it as a view. So we would really loose the advantage of being memory efficient with the as_strided view. In fact, this was the initial answer I got on StackOverflow after which it was corrected as this. I'll however add a comment of 2 lines, explaining in short what einsum does here. Chintak

1 0