diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/indexamajig.c | 145 | ||||
| -rw-r--r-- | src/peaks.c | 144 | ||||
| -rw-r--r-- | src/peaks.h | 1 |
3 files changed, 146 insertions, 144 deletions
diff --git a/src/indexamajig.c b/src/indexamajig.c index 41cf5b22..7c2e1b5d 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -26,6 +26,7 @@ #include "dirax.h" #include "intensities.h" #include "ewald.h" +#include "peaks.h" static void show_help(const char *s) @@ -44,150 +45,6 @@ static void show_help(const char *s) } -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; - 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; -} - - int main(int argc, char *argv[]) { int c; 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; diff --git a/src/peaks.h b/src/peaks.h index d7bbca4d..e914f9b4 100644 --- a/src/peaks.h +++ b/src/peaks.h @@ -18,6 +18,7 @@ #endif +extern int image_fom(struct image *image); extern void search_peaks(struct image *image, int dump_peaks); |