diff options
| -rw-r--r-- | src/indexamajig.c | 94 |
1 files changed, 85 insertions, 9 deletions
diff --git a/src/indexamajig.c b/src/indexamajig.c index 1b1620c1..69813118 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -40,11 +40,21 @@ static void show_help(const char *s) } -static int sum_of_peaks(struct image *image) +struct peak { + int x; + int y; + int i; + int invalid; +}; + + +static 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; @@ -86,8 +96,8 @@ static int sum_of_peaks(struct image *image) th = mean + 5*sd; STATUS("mean=%f ,sd=%f, th=%f\n", mean, sd, th); - /* Sum intensity above threshold */ - integr = 0; + /* Find pixels above threshold */ + n_peaks = 0; for ( x=0; x<1024; x++ ) { for ( y=600; y<1024; y++ ) { @@ -98,12 +108,78 @@ static int sum_of_peaks(struct image *image) val = image->data[x+image->height*y]; - if ( val > th ) integr+=val; + 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++; + } } } - return integr; + do { + + int max, max_i; + 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; + n_valid++; + } + + return n_valid; } @@ -169,7 +245,7 @@ int main(int argc, char *argv[]) char line[1024]; struct hdfile *hdfile; struct image image; - int integr; + int fom; rval = fgets(line, 1023, fh); chomp(line); @@ -185,9 +261,9 @@ int main(int argc, char *argv[]) hdf5_read(hdfile, &image); - integr = sum_of_peaks(&image); - printf("%6i %i\n", n_images, integr); - if ( integr > 200000 ) { + fom = image_fom(&image); + printf("%6i %i\n", n_images, fom); + if ( fom > 0 ) { STATUS("Hit: %s\n", line); |