aboutsummaryrefslogtreecommitdiff
path: root/src/peaks.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-01-12 15:32:09 +0100
committerThomas White <taw@physics.org>2010-01-12 15:32:09 +0100
commitf1fa363013da613fd913d7920d491a31c612fd72 (patch)
tree8617710e358ad77c074d582644465bf0c300e0a5 /src/peaks.c
parent37be1343c612c6abf22814c6eb4b67f7f39c61c8 (diff)
Move hit finder to separate file
Diffstat (limited to 'src/peaks.c')
-rw-r--r--src/peaks.c144
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;