Thanks for the detailed response!<br><br><div class="gmail_quote">On Mon, Nov 12, 2012 at 6:57 PM, Josh Warner <span dir="ltr">&lt;<a href="mailto:silvertrumpet999@gmail.com" target="_blank">silvertrumpet999@gmail.com</a>&gt;</span> wrote:<br>

<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">I&#39;ve looked at these two algorithms, and the biggest difference seems to be in the output.  `is_local_maximum` returns a Boolean array, while `peak_local_max` returns the indices of points corresponding to maxima (a la `np.sort` vs. `np.argsort`, though you cannot just pass the output of `peak_local_max` as indices). <div>

<br></div><div>The other differences are more subtle, but significant.</div><div><br></div><div>The API for `peak_local_max` could use some cleanup - the first threshold kwarg is set to &#39;deprecated&#39;(!) and IMHO should be removed - but this algorithm allows finer control over thresholding and peak searching.</div>

</blockquote><div><br></div><div><br></div><div>I think it&#39;s marked &#39;deprecated&#39; so that it can be removed gracefully (i.e. it gives people time to change their code). That said, it was marked &#39;deprecated&#39; a few releases back, so it&#39;s probably ripe for removal.</div>

<div><br></div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div>This would be good to avoid finding phantom peaks in noise if large dark regions were present. One significant drawback of this algorithm is the min_distance kwarg is set by default to 10, rather arbitrary, and ANY input (even the minimum value of 1) excludes both neighboring pixels AND the border.  See example below.</div>

<div><br></div><div>In contrast, `is_local_maximum` has a much simpler API.  It doesn&#39;t have the finer thresholding / peak searching controls, but has a unique ability to search for peaks ONLY within arbitrary, connected, labeled regions.  This has some interesting potentials for masking etc, though I believe within each label only one peak will be found.  This algorithm also has the ability to search arbitrary local regions for peaks using something akin to a morphological structuring element, through the `footprint=` kwarg.  The documentation for this could probably be clarified.  </div>

<div><br></div><div>The way `peak_local_max` excludes borders concerns me for general use, as does its default `min_distance=10`, and personally I would prefer to wok around the limitations in `is_local_maximum`.  </div>
<div>
<br></div><div>A best-of-both-worlds combination could probably be created without overly much effort...<br></div></blockquote><div><br></div><div><br></div><div>This is a great summary of the API differences. I agree that excluding the border region is a bit of a wart. (My guess is that this could be fixed by changing the boundary condition on the maximum filter used in `peak_local_max`.) Although I agree that defaulting to `min_distance=10` is arbitrary, I&#39;m not sure there&#39;s an obvious choice. I would normally assume that peaks separated by 1 pixel are just noise.</div>

<div><br></div><div>The footprint and mask parameter would definitely be a nice addition to `peak_local_max`.</div><div><br></div><div>If you have time, a pull request to address some or all of these issues would be great. If not, maybe you could this as an issue on github?</div>

<div><br></div><div>Thanks,</div><div>-Tony</div><div><br></div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><br><br>Snippet showing border-excluding behavior of `peak_local_max`, which will only get worse with higher values of `min_distance`:</div>

<div><br></div><div><div class="im"><div>import numpy as np</div><div>import matplotlib.pyplot as plt</div></div><div>from skimage.feature import peak_local_max</div><div>from skimage.morphology import is_local_maximum</div>

<div><br></div><div># Generate standardized random data</div><div>np.random.seed(seed=1234)</div><div>testim = np.random.randint(0, 255, size=(20, 20))</div><div><br></div><div># Find peaks using both methods</div><div>ismax   = is_local_maximum(testim)                  # Boolean image returned</div>

<div>peakmax = peak_local_max(testim, min_distance=1)    # (M, 2) indices returned</div><div><br></div><div># `peakmax` not plottable - placing values in 2d array</div><div>Ipeakmax = np.zeros(testim.shape)</div><div>Ipeakmax[peakmax[:, 0], peakmax[:, 1]] = 1</div>

<div><br></div><div># Show the results</div><div>fig, ax = plt.subplots(ncols=2, nrows=1)</div><div>ax[0].imshow(ismax, cmap=&#39;gray&#39;)</div><div>ax[0].set_title(&#39;Peaks found by `is_local_maximum`&#39;)</div><div>

ax[1].imshow(Ipeakmax, cmap=&#39;gray&#39;)</div><div>ax[1].set_title(&#39;Peaks found by `peak_local_max`&#39;)</div><div><br></div><div>plt.show()</div><br><br><br>On Monday, November 12, 2012 3:57:28 PM UTC-6, Tony S Yu wrote:<div>

<div class="h5"><blockquote class="gmail_quote" style="margin:0;margin-left:0.8ex;border-left:1px #ccc solid;padding-left:1ex"><br><br><div class="gmail_quote">On Mon, Nov 12, 2012 at 7:43 AM, Frank <span dir="ltr">&lt;<a>pennek...@googlemail.com</a>&gt;</span> wrote:<br>

<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">

Dear group,<br><br>I have some issues with the watershed algorithm implemented in scikits image. I use a global threshold to segment cells from background, but some cells touch and I want them to be split. Watershed seems the appropriate way to deal with my problem, however my particles are split in too many pieces. Is there a way to adjust the sensitivity of the watershed method? <br>



<br>Many thanks for any suggestion!<br><br>The code that I use looks like below. An example image that I want to process can be downloaded here: <a href="https://dl.dropbox.com/u/10373933/test.jpg" target="_blank">https://dl.dropbox.com/u/<u></u>10373933/test.jpg</a><br>



<br># packages needed to perform image processing and analysis<br>import numpy as np<br>import scipy as scp<br>import matplotlib.pyplot as plt<br>import matplotlib.image as mpimg<br>import scipy.ndimage as nd<br>import skimage<br>



from skimage import io<br>from skimage.morphology import watershed, is_local_maximum<br>from skimage.segmentation import find_boundaries, visualize_boundaries<br>from skimage.color import gray2rgb<br><br>#read files jpeg file<br>



image = mpimg.imread(&#39;c:\\test.jpg&#39;)<br>image_thresh = image &gt; 140<br>labels = nd.label(image_thresh)[0]<br>distance = nd.distance_transform_edt(<u></u>image_thresh)<br>local_maxi = is_local_maximum(distance, labels=labels, footprint=np.ones((9, 9)))<br>



markers = nd.label(local_maxi)[0]<br>labelled_image = watershed(-distance, markers, mask=image_thresh)<br><br>#find outline of objects for plotting<br>boundaries = find_boundaries(labelled_<u></u>image)<br>img_rgb = gray2rgb(image)<br>



overlay = np.flipud(visualize_<u></u>boundaries(img_rgb,boundaries)<u></u>)<br>imshow(overlay)</blockquote><div><br></div><div>Hi Frank,</div><div><br></div><div>Actually, I don&#39;t think the issue is in the watershed segmentation. Instead, I think the problem is in the marker specification: Using local maxima creates too many marker points when a blob deviates greatly from a circle. (BTW, does anyone know if there are any differences between `is_local_maximum` and `peak_local_max`? Maybe the former should be deprecated.)</div>



<div><br></div><div>Using the centroids of blobs gives cleaner results. See slightly-modified example below.</div><div><br></div><div>Best,</div><div>-Tony</div><div><br></div><div><div># packages needed to perform image processing and analysis</div>



<div>import numpy as np</div><div>import matplotlib.pyplot as plt</div><div>import scipy.ndimage as nd</div><div><br></div><div>from skimage import io</div><div>from skimage import measure</div><div>from skimage.morphology import watershed</div>



<div>from skimage.segmentation import find_boundaries, visualize_boundaries</div><div>from skimage.color import gray2rgb</div><div><br></div><div>#read files jpeg file</div><div>image = io.imread(&#39;test.jpg&#39;)</div>



<div><br></div><div>image_thresh = image &gt; 140</div><div>labels = nd.label(image_thresh)[0]</div><div>distance = nd.distance_transform_edt(<u></u>image_thresh)</div><div><br></div><div>props = measure.regionprops(labels, [&#39;Centroid&#39;])</div>



<div>coords = np.array([np.round(p[&#39;<u></u>Centroid&#39;]) for p in props], dtype=int)</div><div># Create marker image where blob centroids are marked True</div><div>markers = np.zeros(image.shape, dtype=bool)</div><div>

markers[tuple(np.transpose(<u></u>coords))] = True</div>

<div><br></div><div>labelled_image = watershed(-distance, markers, mask=image_thresh)</div><div><br></div><div>#find outline of objects for plotting</div><div>boundaries = find_boundaries(labelled_<u></u>image)</div><div>

img_rgb = gray2rgb(image)</div>

<div>overlay = visualize_boundaries(img_rgb, boundaries, color=(1, 0, 0))</div><div><br></div><div>plt.imshow(overlay)</div><div>plt.show()</div></div></div>
</blockquote></div></div></div>

<p></p>

-- <br>
 <br>
 <br>
</blockquote></div><br>