Need to find 2nd highest peak in a double-peaked graph/histogram

Ray Butler wrote on Sep 04, 2015

Hi,

Normally I can script my way out of trouble, but I've been struggling with this for a while.

I need to be able to automatically find the 2nd highest peak in a double-peaked graph/histogram - like the peak at z ~ 700 in the following screengrab:

I tried ngaussfit with appropriate gausspars, but while it always finds the dominant peak at z ~ 0, it usually fails to find the second peak, which is the one I need. The problem is that this peak moves around a lot between frames, also changing in width and amplitude; so the initial guesses for gausspars.ampl2, gausspars.cent2 and gausspars.fwhm2 either have to be ridiculously oversized (the fit returned is WAY off), or reasonable but locating only the few frames where the peak falls in the "right" place.

It's one of those situations where to the human eye, the whereabouts of the 2nd peak is always really obvious, but making a code "see" it is trickier.

Any help appreciated.

Thanks,

Ray

---

Mike Fitzpatrick wrote on Sep 05, 2015

Try the script below, adjust the code as needed. The comments should explain what it does but post if you have questions.

procedure peak2 (imname)

string  imname                          { prompt = "image name"         }
int     nbin    = 128                   { prompt = "number of bins"     }

begin
        string  limname, hname, pname, ans
        real    bins[1024], peaks[1024]
        int     lnbin

        limname = imname                            # get params to local vars
        lnbin = nbin

        hname = mktemp ("tmp$hist")                 # create temp filenames
        pname = mktemp ("tmp$peaks")
        ans = mktemp ("tmp$ans")

        imhistogram (imname, listout=yes, > hname)  # list out histogram

        # Read the histogram into a local array.
        list = hname
        for (i=1; fscan (list, x, y) != EOF; i=i+1) {
            bins[i] = x ; peaks[i] = y
        }

        # Find the peaks, where a 'peak' is higher than both it's neighbors
        # points a little farther away to find broader peaks.
        for (i=2; i < (lnbin-5); i=i+1) {
            if (i > 5) {
              if (peaks[i] > peaks[i-1] && peaks[i] > peaks[i+1] &&
                  peaks[i] > peaks[i-4] && peaks[i] > peaks[i+4])
                    printf ("%g %g\n", bins[i], peaks[i], >> pname)
            } else {
              if (peaks[i] > peaks[i-1] && peaks[i] > peaks[i+1])
                    printf ("%g %g\n", bins[i], peaks[i], >> pname)
            }
        }

        # Sort the peaks list and pull out the second result.
        sort (pname,column=2,numeric+,reverse+) | head ("STDIN", nl=2, > ans)
        tail (ans, nl=1)

        delete (hname, ver-) ; delete (pname, ver-) ; delete (ans, ver-)
end
 

---

Last post on Sep 05, 2015