peak detection in (x,y,z) data
in [Matlab]
|
From: highegg on 26 Sep 2009 13:16 On 25 záÅ, 22:25, Rahul <nos...(a)nospam.invalid> wrote: > I have a numerical (x,y,z) dataset with multiple peaks (local maxima) in it > that I want to detect. I've seen many routines for (x,y) peak detection but > none that I know so far do it for (x,y,z) data. > > Any routines around that I am missing? > > -- > Rahul The problem basically is that in 3D, it's no longer so simple to define a peak of the sample, unless your data is on a regular grid. If it suffices for your needs, you can first interpolate your data to a grid, using griddata or some OctaveForge package (such as OctGPR or data-smoothing), then find peaks on the grid. On a grid, you can simply find the peaks as those points that are either less or greater than all of their neighbours. hth
From: Rahul on 26 Sep 2009 16:50 highegg <highegg(a)gmail.com> wrote in news:a7f3394e-7f16-4fe1-93a9- bbca8702505a(a)o36g2000vbl.googlegroups.com: > On a grid, you can simply find the peaks as those points that are > either less or greater than all of their neighbours. > Good point. I guess first I need to pass it through some sort of noise reduction filter. Otherwise I'll have every tiny bump of noise on a scanline show up as a peak. Correct? -- Rahul
From: Rahul on 26 Sep 2009 16:55 "Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in news:h9jc0f$5vg$1(a)fred.mathworks.com: >% Data > A=rand(100,100,100); > >% Engine > Amin=minmaxfilt(A,3,'min','same'); % alternatively use imerode in > image processing [i j k]=ind2sub(size(A),find(Amin==A)); % <- index of > local minima > >% Remove index in the borders if interested in inside local minima only >% Keep this steps otherwise > idxkeep=find(i>1 & i<100 & j>1 & j<100 & k>1 & k<100); > i=i(idxkeep); > j=j(idxkeep); > k=k(idxkeep); > > I guess I still haven't understood the usage. Let's say I have make a matrix with a single sharp peaks: A=ones(10,10) A(5,5)=40 A(5,4)=30 A(5,6)=30 A(4,5)=30 A(6,5)=30 A(4,4)=20 A(6,6)=20 A(4,6)=20 A(6,4)=20 Now how do I get the function to tell me (5,5) by operating on A. 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 20 30 20 1 1 1 1 1 1 1 30 40 30 1 1 1 1 1 1 1 20 30 20 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 I tried your invocation and varients on it but I didn't really understand it. -- Rahul
From: Bruno Luong on 27 Sep 2009 03:07 Rahul <nospam(a)nospam.invalid> wrote in message <Xns9C92A20228BE66650A1FC0D7811DDBC81(a)188.40.43.213>... > > 1 1 1 1 1 1 1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 > 1 1 1 20 30 20 1 1 1 1 > 1 1 1 30 40 30 1 1 1 1 > 1 1 1 20 30 20 1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 This example does look simple, but it's really not for what one would use to understand in the context of peak detection. It contains many (degenerated) maxima: all the 1s surounding by 1s. In addition to the big peak, they all will be detected by the method I propose, I have no idea how your data look like. Usually experimental data will contains noise, and it will introduce many maxima that you need to filter out after detection based on the expected amplitude of the peak. If you want simple data, I propose to use smooth data come from PEAKS function. After you understand how it works, then tackle the noise and degenerated cases. A = peaks; Amax=minmaxfilt(A,3,'max','same'); ip = find(Amax==A); [y x]=ind2sub(size(A),ip); idxkeep=find(x>1 & x<size(A,2) & y>1 & y<size(A,1)); x=x(idxkeep); y=y(idxkeep); ip=ip(idxkeep); figure(1); clf; surf(A); hold on; plot3(x,y,A(ip),'or','Markersize',10,'Linewidth',3); Bruno
From: ImageAnalyst on 27 Sep 2009 14:05 On Sep 26, 4:55 pm, Rahul <nos...(a)nospam.invalid> wrote: > "Bruno Luong" <b.lu...(a)fogale.findmycountry> wrote innews:h9jc0f$5vg$1(a)fred.mathworks.com: > > > > > > >% Data > > A=rand(100,100,100); > > >% Engine > > Amin=minmaxfilt(A,3,'min','same'); % alternatively use imerode in > > image processing [i j k]=ind2sub(size(A),find(Amin==A)); % <- index of > > local minima > > >% Remove index in the borders if interested in inside local minima only > >% Keep this steps otherwise > > idxkeep=find(i>1 & i<100 & j>1 & j<100 & k>1 & k<100); > > i=i(idxkeep); > > j=j(idxkeep); > > k=k(idxkeep); > > I guess I still haven't understood the usage. > > Let's say I have make a matrix with a single sharp peaks: > > A=ones(10,10) > > A(5,5)=40 > > A(5,4)=30 > A(5,6)=30 > A(4,5)=30 > A(6,5)=30 > > A(4,4)=20 > A(6,6)=20 > A(4,6)=20 > A(6,4)=20 > > Now how do I get the function to tell me > > (5,5) by operating on A. > > 1 1 1 1 1 1 1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 > 1 1 1 20 30 20 1 1 1 1 > 1 1 1 30 40 30 1 1 1 1 > 1 1 1 20 30 20 1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 > 1 1 1 1 1 1 1 1 1 1 > > I tried your invocation and varients on it but I didn't really understand > it. > > -- > Rahul- -------------------------------------------------------------------------------------------------------------- Rahul: What do you consider the "peak"? Just the apex, or the complete peak out to the background level of 1? Here is some code to calculate it both ways. It will work for your simple case but probably isn't robust for real world cases with complicated peaks and noise: -ImageAnalyst A=ones(10,10) A(5,5)=40 A(5,4)=30 A(5,6)=30 A(4,5)=30 A(6,5)=30 A(4,4)=20 A(6,6)=20 A(4,6)=20 A(6,4)=20 % Get the complete peak out to the uniform background. completePeak = A>1 % Get just the apex: the pixel with the value of 40. em = imextendedmax(A, 2) % Get local maxima em2 = A .* em % Get the values from the original image. JustTheApex = (em2 == A)
First
|
Prev
|
Next
|
Last
Pages: 1 2 3 4 5 Prev: finding eigenvalues of a sparce matrix Next: Removing a data cursor programmatically |