aboutsummaryrefslogtreecommitdiff
path: root/src/indexamajig.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/indexamajig.c
parent37be1343c612c6abf22814c6eb4b67f7f39c61c8 (diff)
Move hit finder to separate file
Diffstat (limited to 'src/indexamajig.c')
-rw-r--r--src/indexamajig.c145
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;