diff options
Diffstat (limited to 'src/indexamajig.c')
-rw-r--r-- | src/indexamajig.c | 145 |
1 files changed, 1 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; |