| raybutler |
09/04/2015 08:10AM (Read 2051 times)
|
|
|
Status: offline
Registered: 06/24/2008
Posts: 7
|
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
Dr. Ray Butler
School of Physics (Lecturer) | Centre for Astronomy
National University of Ireland - Galway, Ireland
|
| |
|
|
| fitz |
|
09/05/2015 10:03PM
|
|
|
Status: offline
Registered: 09/30/2005
Posts: 4040
|
Try the script below, adjust the code as needed. The comments should explain what it does but post if you have questions.
PHP Formatted Code
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
|
| |
|
|
| |
|
New Topic
Post Reply
Need to find 2nd highest peak in a double-peaked graph/histogram
Quote
