Re: [Numpy-discussion] Functions for finding the relative extrema of numeric data
What is your application?
The most common case is looking at Fourier transforms and identifying spectral peaks. I've also analyzed images looking at 1D slices (usually very regular data) and looked for peaks there.
That stackoverflow page had a nice link to a comparison of different algorithms here: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2631518/. That paper is focused on mass-spectrometry data, but the approach would generalize to any 1D data set. Unless somebody feels otherwise, I'll close this pull request and start working on an implementation of peak finding via continuous wavelet transform (the best and most computationally intensive approach of those described above).
-Jacob
------------------------------
Message: 4 Date: Tue, 13 Sep 2011 22:34:01 +0200 From: Ralf Gommers ralf.gommers@googlemail.com Subject: Re: [Numpy-discussion] Functions for finding the relative extrema of numeric data To: Discussion of Numerical Python numpy-discussion@scipy.org Message-ID: <CABL7CQhxCX0LKFENMW6-4ZSbdieGxz04zbsrnY4bXYVxVL78Dw@mail.gmail.com
Content-Type: text/plain; charset="iso-8859-1"
Hi Jacob,
On Fri, Sep 9, 2011 at 11:57 PM, Jacob Silterra jsilter@gmail.com wrote:
Hello all,
I'd like to see functions for calculating the relative extrema in a set
of
data included in numpy. I use that functionality frequently, and always
seem
to be writing my own version. It seems like this functionality would be useful to the community at large, as it's a fairly common operation.
What is your application?
For numeric data (which is presumably noisy), the definition of a
relative
extrema isn't completely obvious. The implementation I am proposing finds
a
point in an ndarray along an axis which is larger (or smaller) than it's `order` nearest neighbors (`order` being an optional parameter, default
1).
This is likely to find more points than may be desired, which I believe
is
preferable to the alternative. The output is formatted the same as numpy.where.
Code available here: https://github.com/numpy/numpy/pull/154
I'm not sure whether this belongs in numpy or scipy, that question is somewhat debatable. More sophisticated peak-finding functions (in N dimensions, as opposed to 1) may also be useful to the community, and
those
would definitely belong in scipy.
I have the feeling this belongs in scipy. Although if it's just these two functions I'm not sure where exactly to put them. Looking at the functionality, this is quite a simple approach. For any data of the type I'm usually working with it will not be able to find the right local extrema. The same is true for your alternative definition below.
A more powerful peak detection function would be a very good addition to scipy imho (perhaps in scipy.interpolate?). See also
http://stackoverflow.com/questions/1713335/peak-finding-algorithm-for-python...
Cheers, Ralf
An alternative implementation would be to require that function be continuously descending (or ascending) for `order` points, which would enforce a minimum width on a peak.
-Jacob Silterra
NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
I've written some peak picking functions that work in N dimensions for a module for looking at NMR data in python, http://code.google.com/p/nmrglue/. I'd be glad to polish up the code if people think it would be a useful addition to scipy.ndimage or scipy.interpolate? The methods are not based on any formal algorithms I know of, just some fast and relatively simple methods that I found seem to work decently.
The methods are contained in the peakpick.py and segmentation.py files in the analysis directory (specifically see the find_all_connected, find_all_ downward and find_all_thres): http://code.google.com/p/nmrglue/source/browse/trunk/nmrglue/analysis/peakpi... http://code.google.com/p/nmrglue/source/browse/trunk/nmrglue/analysis/segmen...
Let me know if there is an interest in including these in scipy or numpy.
-Jonathan Helmus
Jacob Silterra wrote:
What is your application?
The most common case is looking at Fourier transforms and identifying spectral peaks. I've also analyzed images looking at 1D slices (usually very regular data) and looked for peaks there.
That stackoverflow page had a nice link to a comparison of different algorithms here: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2631518/. That paper is focused on mass-spectrometry data, but the approach would generalize to any 1D data set. Unless somebody feels otherwise, I'll close this pull request and start working on an implementation of peak finding via continuous wavelet transform (the best and most computationally intensive approach of those described above).
-Jacob
------------------------------ Message: 4 Date: Tue, 13 Sep 2011 22:34:01 +0200 From: Ralf Gommers <ralf.gommers@googlemail.com <mailto:ralf.gommers@googlemail.com>> Subject: Re: [Numpy-discussion] Functions for finding the relative extrema of numeric data To: Discussion of Numerical Python <numpy-discussion@scipy.org <mailto:numpy-discussion@scipy.org>> Message-ID: <CABL7CQhxCX0LKFENMW6-4ZSbdieGxz04zbsrnY4bXYVxVL78Dw@mail.gmail.com <mailto:CABL7CQhxCX0LKFENMW6-4ZSbdieGxz04zbsrnY4bXYVxVL78Dw@mail.gmail.com>> Content-Type: text/plain; charset="iso-8859-1" Hi Jacob, On Fri, Sep 9, 2011 at 11:57 PM, Jacob Silterra <jsilter@gmail.com <mailto:jsilter@gmail.com>> wrote: > Hello all, > > I'd like to see functions for calculating the relative extrema in a set of > data included in numpy. I use that functionality frequently, and always seem > to be writing my own version. It seems like this functionality would be > useful to the community at large, as it's a fairly common operation. > What is your application? > > For numeric data (which is presumably noisy), the definition of a relative > extrema isn't completely obvious. The implementation I am proposing finds a > point in an ndarray along an axis which is larger (or smaller) than it's > `order` nearest neighbors (`order` being an optional parameter, default 1). > This is likely to find more points than may be desired, which I believe is > preferable to the alternative. The output is formatted the same as > numpy.where. > > Code available here: https://github.com/numpy/numpy/pull/154 > > I'm not sure whether this belongs in numpy or scipy, that question is > somewhat debatable. More sophisticated peak-finding functions (in N > dimensions, as opposed to 1) may also be useful to the community, and those > would definitely belong in scipy. > I have the feeling this belongs in scipy. Although if it's just these two functions I'm not sure where exactly to put them. Looking at the functionality, this is quite a simple approach. For any data of the type I'm usually working with it will not be able to find the right local extrema. The same is true for your alternative definition below. A more powerful peak detection function would be a very good addition to scipy imho (perhaps in scipy.interpolate?). See also http://stackoverflow.com/questions/1713335/peak-finding-algorithm-for-python-scipy Cheers, Ralf > An alternative implementation would be to require that function be > continuously descending (or ascending) for `order` points, which would > enforce a minimum width on a peak. > > -Jacob Silterra > > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion@scipy.org <mailto:NumPy-Discussion@scipy.org> > http://mail.scipy.org/mailman/listinfo/numpy-discussion > >
On 09/15/2011 03:32 PM, Jacob Silterra wrote:
I'll close this pull request and start working on an implementation of peak finding via continuous wavelet transform (the best and most computationally intensive approach of those described above).
Just for information, which tools are you going to use for the CWT? I may be interested in providing some help.
Cheers,
Davide Lasagna
Hi all,
I am not sure if this is of help for anyone. I wrote some code to find the relative maxima in a 1D array for my own purpose. Maybe someone is interested or even finds a bug *g*. I post the code here and appreciate any feedback. Even "stop spamming your buggy code" :-)
from numpy import diff, sign, convolve, array, where, around, int32, alen
def localmaxima_at(x): '''Returns the indices of local maxima in the 1D array x.
If several elements in x have the same value, then the index of the element in the middle is returned. If there are two adjacent elements with the same value, one of them is returned. x[0] and x[-1] are never returned as an index for the local maximum. @Author: Samuel John @copyright: http://creativecommons.org/licenses/by-nc-sa/3.0/ @todo: unittests ''' assert len(x) > 2, "Length of x should be greater than two in order to define a meaningful local maximum." assert x.ndim == 1, "Expected 1D array." #print 'x=\n',x filled=sign(diff(x)).astype(int32) # fill zeros: has_zeros = (filled == 0).any() last = 0 if has_zeros: for i in xrange(alen(filled)): if filled[i] == 0: filled[i] = last else: last = filled[i] #print 'filled\n',filled left = where( convolve( filled, array([-1,1]), mode='full' ) -2 == 0 )[0] if has_zeros: filled=sign(diff(x)).astype(int32) last = 0 for i in reversed(xrange(len(filled))): if filled[i] == 0: filled[i] = last else: last = filled[i] right = where( convolve( filled, array([-1,1]), mode='full' ) -2 == 0 )[0] #print 'left\n',left #print 'right\n',right return around( (left + right) / 2.0).astype(int32)
bests, Samuel
participants (4)
-
Davide -
Jacob Silterra -
Jonathan Helmus -
Samuel John