diff options
author | Thomas White <taw@physics.org> | 2010-01-12 15:32:09 +0100 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2010-01-12 15:32:09 +0100 |
commit | f1fa363013da613fd913d7920d491a31c612fd72 (patch) | |
tree | 8617710e358ad77c074d582644465bf0c300e0a5 /src/peaks.c | |
parent | 37be1343c612c6abf22814c6eb4b67f7f39c61c8 (diff) |
Move hit finder to separate file
Diffstat (limited to 'src/peaks.c')
-rw-r--r-- | src/peaks.c | 144 |
1 files changed, 144 insertions, 0 deletions
diff --git a/src/peaks.c b/src/peaks.c index afc05401..984d4d3b 100644 --- a/src/peaks.c +++ b/src/peaks.c @@ -26,6 +26,150 @@ #define PEAK_WINDOW_SIZE (10) +struct peak { + int x; + int y; + int i; + int invalid; +}; + + +int image_fom(struct image *image) +{ + int x, y; + int integr, n; + float fintegr, mean, sd, th; + struct peak peaks[1024]; + int i, n_peaks, n_nearby, n_valid; + + /* Measure mean */ + integr = 0; + n = 0; + for ( x=0; x<1024; x++ ) { + for ( y=600; y<1024; y++ ) { + + int val; + if ( (x>400) && (x<600) ) continue; + val = image->data[x+image->height*y]; + if ( val < 0 ) continue; + integr += val; + n++; + + } + } + mean = (float)integr / n; /* As integer to keep maths fast */ + + /* Standard deviation */ + fintegr = 0; + for ( x=0; x<1024; x++ ) { + for ( y=600; y<1024; y++ ) { + + float val; + + if ( (x>400) && (x<600) ) continue; + val = (float)image->data[x+image->height*y]; + if ( val < 0 ) continue; + + val -= mean; + val = powf(val, 2.0); + fintegr += val; + + } + } + sd = sqrtf(fintegr / n); + + /* Threshold */ + th = mean + 5*sd; + STATUS("mean=%f ,sd=%f, th=%f\n", mean, sd, th); + + /* Find pixels above threshold */ + n_peaks = 0; + for ( x=0; x<1024; x++ ) { + for ( y=600; y<1024; y++ ) { + + int val; + + /* Chop out streaky region */ + if ( (x>400) && (x<600) ) continue; + + val = image->data[x+image->height*y]; + + if ( val > th ) { + peaks[n_peaks].x = x; + peaks[n_peaks].y = y; + peaks[n_peaks].i = val; + peaks[n_peaks].invalid = 0; + n_peaks++; + } + + } + } + + do { + + int max, max_i = -1; + int adjacent; + + n_nearby = 0; + + /* Find brightest peak */ + max = 0; + for ( i=0; i<n_peaks; i++ ) { + if ( peaks[i].i > max ) { + max = peaks[i].i; + max_i = i; + } + } + + /* Must be at least one adjacent bright pixel */ + adjacent = 0; + for ( i=0; i<n_peaks; i++ ) { + + int td; + + td = abs(peaks[i].x - peaks[max_i].x) + + abs(peaks[i].y - peaks[max_i].y); + if ( td == 1 ) { + adjacent++; + break; /* One is enough */ + } + } + if ( adjacent < 1 ) { + peaks[i].invalid = 1; + continue; + } + + /* Remove nearby (non-invalidated) peaks from list */ + n_nearby = 0; + for ( i=0; i<n_peaks; i++ ) { + + int dx, dy, ds; + + if ( peaks[i].invalid ) continue; + + dx = abs(peaks[i].x - peaks[max_i].x); + dy = abs(peaks[i].y - peaks[max_i].y); + ds = dx*dx + dy*dy; + if ( ds < 36 ) { + n_nearby++; + peaks[i].invalid = 1; + } + + } + + } while ( n_nearby ); + + n_valid = 0; + for ( i=0; i<n_peaks; i++ ) { + if ( peaks[i].invalid ) continue; + //printf("%i %i\n", peaks[i].x, peaks[i].y); + n_valid++; + } + + return n_valid; +} + + void search_peaks(struct image *image, int dump_peaks) { int x, y, width, height; |