Finding Edges

Share icon

Tue, 09/30/2008 - 03:05 pm

I am having a devil of a time coming up with a way to properly identify peak edges in my data. I'd like to use FindPeak, but I don't understand how it works very well. My data consists of a long wave with amplitude data- there are a series of "peaks" (actually 1-D MRI images), but they don't have any conventional shapes like a Gaussian or Lorentzian. I can get FindPeak to get the peak value right, but I have little to no interest in that. It's the edges that I'm after, and up till now, I've been identifying them by sight. This is no longer feasible. FindPeak has the nasty habit of telling me peak edges that are no where near correct. I think it is overly sensitive to what defines a peak. My data is a tad noisy, but what FindPeak is giving me is absolutely ridiculous. I am looking at a "peak" right now that is about 500 points wide, and Igor is telling me it's about 5 points wide. That makes no sense. It appears to be identifying amplitude swings from noise as actual peak edges. I don't understand the box size functionality, and just trying random numbers for the /b=value does not seem to help.

Does anyone have any suggestions on how to find peak edges? I am including an example of my data to better show what I am talking about (as both an image and an Igor binary). I'd visually identify points 2205 (roughly) and 2747 as the edges of this peak. This is a snippet of an actual dataset, which consists of a series of these guys, one right after the other. I'm sorry to ask such a basic question, but I've already invested a day on this, and am getting nowhere. I appreciate your help and feedback very much! Peak1.ibw
This should be straightforward. The main thing is to filter the high-frequency noise before you look for peaks. There is such a difference in scale between the noise and the real peaks that it should be easy. I tried FindPeak and got the results you did. Maybe something is off with that operation? In any case, the code below works. Basically, one smooths, differentiates, and looks for the positions of the two peaks. The only thing wrong with the code here is that Wavestats rounds max and min locations to integer point numbers. One could use FindPeak or FindPeaks and get the peak position more accurately that way. Probably a better but more complicated way is to use wavestats to get the approximate peaks of the derivative. Then fit a low (even) order polynomial to the points surrounding the max or min location (left or right edges). Then get the more exact position from the fit coefficients.

Hope this helps! (Be sure to display w1 and w2 along with your original data to see what is going on.)

Function FindWidth(w, nsmooth)
    wave w; variable nsmooth
    Duplicate /o w w1
    Smooth nsmooth, w1
    Differentiate w1 /D=w2
    wavestats /q w2
    print abs(V_minloc-V_maxloc)
//  Killwaves w1,w2
End

John Bechhoefer
Department of Physics
Simon Fraser University
Burnaby, BC, Canada

September 30, 2008 at 06:56 pm - Permalink

You might also find the PulseStats operation useful. For example, I loaded your wave and did the following:
wavestats peak1
  V_npnts= 4096; V_numNaNs= 0; V_numINFs= 0; V_avg= 677.141;
  V_Sum= 2.77357e+06; V_sdev= 1922.9; V_rms= 2038.42;
  V_adev= 1218.66; V_skew= 2.5337; V_kurt= 4.69783;
  V_minloc= 2139; V_maxloc= 2486; V_min= -118.567; V_max= 7024.83;
  V_startRow= 0; V_endRow= 4095;

pulsestats/b=55 /f=0.99 /L=(V_max, peak1[V_endrow]) /M=(pnt2x(peak1, V_maxLoc) - pnt2x(peak1, V_minLoc)) Peak1
  V_Flag= 1; V_PulseLoc1= 2218.27; V_PulseLoc2= 2750.07;
  V_PulseLoc3= NaN; V_PulsePolarity= 1; V_PulseLvl0= 7024.83;
  V_PulseLvl123= 681.587; V_PulseLvl4= -23.2181;
  V_PulseAmp4_0= -7048.05; V_PulseWidth2_1= 531.797;
  V_PulseWidth3_2= NaN; V_PulseWidth3_1= NaN;
  V_Flag= 1; V_PulseLoc1= 2184.2; V_PulseLoc2= 2768.5;
  V_PulseLoc3= NaN; V_PulsePolarity= 1; V_PulseLvl0= 7024.83;
  V_PulseLvl123= 47.2623; V_PulseLvl4= -23.2181;
  V_PulseAmp4_0= -7048.05; V_PulseWidth2_1= 584.295;
  V_PulseWidth3_2= NaN; V_PulseWidth3_1= NaN;


Depending on what values you use for the /F flag, the V_PulseWidth2_1 value gets pretty close to 500, and the V_PulseLoc2 and V_PulseLoc1 values are pretty close to what one gets by hand. This might not be the best way to do this, but it might be helpful.

September 30, 2008 at 08:08 pm - Permalink

Thank you both for your input. I have tried both methods, and for my data and my needs, I believe the PulseStats command works best. It returns edges that are about what I pick by eye when using /f=.95. It was hard for me to see how pulsestats could be used on data like mine (just 2 edges), but after seeing the code, and the values in returned, and reading through the manual (specifically mentions two edges), I think I have a handle on it.

October 1, 2008 at 09:01 am - Permalink

This is the function that I came up with. It is still a work in progress.
<br />
Function GradCal(w, Size, Points1D, Points2D, Threshhold, Box, grad)<br />
    wave w; variable Size, Points1D, Points2D, Threshhold, Box; string grad; <br />
    variable i<br />
    string GradLeft=grad+"Left", GradRight=grad+"Right", GradCalib=grad+"Calib"<br />
    string GradBW=grad+"BW", GradCent=grad+"Cent"<br />
    Make /O/N=(Points2D)/D $GradLeft, $GradRight, $GradCalib, $GradBW, $GradCent<br />
    Wave Left=$(GradLeft), Right=$(GradRight), Calib=$(GradCalib), BW=$(GradBW), Cent=$(GradCent);<br />
    <br />
    for(i=0;i<Points2D;i+=1)<br />
        wavestats /R=[Points1D*i,Points1D*(i+1)-1] /Q w<br />
        pulsestats/b=(Box) /f=(Threshhold) /L=(V_max, mean(w,Points1D*i,Points1D/2+Points1D*i)) /R=[Points1D*i,Points1D*(i+1)-1] /M=20 /Q w<br />
//      printf "%g, %g\r" V_PulseLoc1, V_PulseLoc2<br />
        if(i<Points2D/2-1)<br />
            Left[i]=V_PulseLoc1; Right[i]=V_PulseLoc2<br />
        else<br />
            Left[i]=V_PulseLoc2; Right[i]=V_PulseLoc1<br />
        endif<br />
    endfor<br />
    <br />
    BW=(Left-Right)*31250/Points1D<br />
    Calib=BW/(42.5775*Size)<br />
    <br />
End<br />

October 1, 2008 at 02:11 pm - Permalink

No one seems to have mentioned the EdgeStats operation built into Igor. I have found it useful in dealing with square-ish pulses; perhaps it may work on shorter peaks too, if you properly specify the search regions.

Stephen R. Chinn

October 2, 2008 at 03:39 am - Permalink

Is there a significant difference between EdgeStats and PulseStats? Is EdgeStats just PulseStats for one edge, without returning values of NaN?

October 2, 2008 at 07:05 am - Permalink

patarroyo wrote:
Is there a significant difference between EdgeStats and PulseStats? Is EdgeStats just PulseStats for one edge, without returning values of NaN?


Yes, Edgestats is pretty much PulseStats for one edge.

Regarding FindPeak, it's rather primitive and needs an upgrade. From the help:


FindPeak computes the sliding average of the input wave using the BoxSmooth algorithm with the box parameter. The peak center is found where the derivative of this smoothed result crosses zero. The peak edges are found where the second derivative of the smoothed result crosses zero. Linear interpolation of the derivatives is used to more precisely locate the center and edges. The peak value is simply the larger of the two unsmoothed values surrounding the peak center (if /N, the smaller value).

As you can see, FindPeak is not a high-accuracy measurement routine; it is intended as a simple peak-finder. Use the PulseStats operation for more precise statistics.

A simple automatic peak-finder is implemented in the procedure file:
#include <Peak AutoFind>



What it lacks are estimates of the expected width, like PulseStats' /M=dx, and amplitude, like FindLevels.


Software Engineer, WaveMetrics, Inc.

October 2, 2008 at 11:27 am - Permalink