aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2011-02-22 01:18:05 -0800
committerThomas White <taw@physics.org>2012-02-22 15:27:15 +0100
commitb57080f25fc077aa8289f20c6856dddb32101f8a (patch)
tree978ff08a3c8fc2a111fb40c936969990be8d650c
parent06ffec86ab14cedb4880158a2161fd0ad616d302 (diff)
Detect peaks in each panel individually
This avoids straying off a panel when doing the search.
-rw-r--r--src/peaks.c129
1 files changed, 68 insertions, 61 deletions
diff --git a/src/peaks.c b/src/peaks.c
index 99aa4539..676a4387 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -275,14 +275,15 @@ int integrate_peak(struct image *image, int xp, int yp,
}
-void search_peaks(struct image *image, float threshold, float min_gradient)
+static void search_peaks_in_panel(struct image *image, float threshold,
+ float min_gradient, struct panel *p)
{
- int x, y, width, height;
+ int fs, ss, stride;
float *data;
double d;
int idx;
- float fx = 0.0;
- float fy = 0.0;
+ float f_fs = 0.0;
+ float f_ss = 0.0;
float intensity = 0.0;
int nrej_dis = 0;
int nrej_hot = 0;
@@ -293,42 +294,28 @@ void search_peaks(struct image *image, float threshold, float min_gradient)
int ncull;
data = image->data;
- width = image->width;
- height = image->height;
-
- if ( image->features != NULL ) {
- image_feature_list_free(image->features);
- }
- image->features = image_feature_list_new();
+ stride = image->width;
- for ( x=1; x<image->width-1; x++ ) {
- for ( y=1; y<image->height-1; y++ ) {
+ for ( fs = p->min_fs+1; fs <= p->max_fs-1; fs++ ) {
+ for ( ss = p->min_ss+1; ss <= p->max_ss-1; ss++ ) {
double dx1, dx2, dy1, dy2;
double dxs, dys;
double grad;
- int mask_x, mask_y;
- int sx, sy;
+ int mask_fs, mask_ss;
+ int s_fs, s_ss;
double max;
unsigned int did_something;
int r;
- struct panel *p;
/* Overall threshold */
- if ( data[x+width*y] < threshold ) continue;
-
- p = find_panel(image->det, x, y);
- if ( !p ) continue;
- if ( p->no_index ) continue;
-
- /* Ignore streak */
- if ( in_streak(x, y) ) continue;
+ if ( data[fs+stride*ss] < threshold ) continue;
/* Get gradients */
- dx1 = data[x+width*y] - data[(x+1)+width*y];
- dx2 = data[(x-1)+width*y] - data[x+width*y];
- dy1 = data[x+width*y] - data[(x+1)+width*(y+1)];
- dy2 = data[x+width*(y-1)] - data[x+width*y];
+ dx1 = data[fs+stride*ss] - data[(fs+1)+stride*ss];
+ dx2 = data[(fs-1)+stride*ss] - data[fs+stride*ss];
+ dy1 = data[fs+stride*ss] - data[(fs+1)+stride*(ss+1)];
+ dy2 = data[fs+stride*(ss-1)] - data[fs+stride*ss];
/* Average gradient measurements from both sides */
dxs = ((dx1*dx1) + (dx2*dx2)) / 2;
@@ -339,25 +326,29 @@ void search_peaks(struct image *image, float threshold, float min_gradient)
if ( grad < min_gradient ) continue;
- mask_x = x;
- mask_y = y;
+ mask_fs = fs;
+ mask_ss = ss;
do {
- max = data[mask_x+width*mask_y];
+ max = data[mask_fs+stride*mask_ss];
did_something = 0;
- for ( sy=biggest(mask_y-PEAK_WINDOW_SIZE/2, 0);
- sy<smallest(mask_y+PEAK_WINDOW_SIZE/2, height-1);
- sy++ ) {
- for ( sx=biggest(mask_x-PEAK_WINDOW_SIZE/2, 0);
- sx<smallest(mask_x+PEAK_WINDOW_SIZE/2, width-1);
- sx++ ) {
-
- if ( data[sx+width*sy] > max ) {
- max = data[sx+width*sy];
- mask_x = sx;
- mask_y = sy;
+ for ( s_ss=biggest(mask_ss-PEAK_WINDOW_SIZE/2,
+ p->min_ss);
+ s_ss<=smallest(mask_ss+PEAK_WINDOW_SIZE/2,
+ p->max_ss);
+ s_ss++ ) {
+ for ( s_fs=biggest(mask_fs-PEAK_WINDOW_SIZE/2,
+ p->min_fs);
+ s_fs<=smallest(mask_fs+PEAK_WINDOW_SIZE/2,
+ p->max_fs);
+ s_fs++ ) {
+
+ if ( data[s_fs+stride*s_ss] > max ) {
+ max = data[s_fs+stride*s_ss];
+ mask_fs = s_fs;
+ mask_ss = s_ss;
did_something = 1;
}
@@ -365,35 +356,30 @@ void search_peaks(struct image *image, float threshold, float min_gradient)
}
/* Abort if drifted too far from the foot point */
- if ( distance(mask_x, mask_y, x, y) > p->peak_sep ) {
+ if ( distance(mask_fs, mask_ss, fs, ss)
+ > p->peak_sep ) {
break;
}
} while ( did_something );
/* Too far from foot point? */
- if ( distance(mask_x, mask_y, x, y) > p->peak_sep ) {
+ if ( distance(mask_fs, mask_ss, fs, ss) > p->peak_sep ) {
nrej_dis++;
continue;
}
/* Should be enforced by bounds used above. Muppet check. */
- assert(mask_x < image->width);
- assert(mask_y < image->height);
- assert(mask_x >= 0);
- assert(mask_y >= 0);
-
- /* Isolated hot pixel? */
- if ( is_hot_pixel(image, mask_x, mask_y) ) {
- nrej_hot++;
- continue;
- }
+ assert(mask_fs <= p->max_fs);
+ assert(mask_ss <= p->max_ss);
+ assert(mask_fs >= p->min_fs);
+ assert(mask_ss >= p->min_ss);
/* Centroid peak and get better coordinates.
* Don't bother doing polarisation/SA correction, because the
* intensity of this peak is only an estimate at this stage. */
- r = integrate_peak(image, mask_x, mask_y,
- &fx, &fy, &intensity, NULL, NULL, 0, 1);
+ r = integrate_peak(image, mask_fs, mask_ss,
+ &f_fs, &f_ss, &intensity, NULL, NULL, 0, 1);
if ( r ) {
/* Bad region - don't detect peak */
nrej_bad++;
@@ -401,21 +387,21 @@ void search_peaks(struct image *image, float threshold, float min_gradient)
}
/* It is possible for the centroid to fall outside the image */
- if ( (fx < 0.0) || (fx > image->width)
- || (fy < 0.0) || (fy > image->height) ) {
+ if ( (f_fs < p->min_fs) || (f_fs > p->max_fs)
+ || (f_ss < p->min_ss) || (f_ss > p->max_ss) ) {
nrej_fra++;
continue;
}
/* Check for a nearby feature */
- image_feature_closest(image->features, fx, fy, &d, &idx);
+ image_feature_closest(image->features, f_fs, f_ss, &d, &idx);
if ( d < p->peak_sep ) {
nrej_pro++;
continue;
}
/* Add using "better" coordinates */
- image_add_feature(image->features, fx, fy, image, intensity,
+ image_add_feature(image->features, f_fs, f_ss, image, intensity,
NULL);
nacc++;
@@ -431,12 +417,33 @@ void search_peaks(struct image *image, float threshold, float min_gradient)
ncull = 0;
}
- STATUS("%i accepted, %i box, %i hot, %i proximity, %i outside frame, "
+ STATUS("%i accepted, %i box, %i hot, %i proximity, %i outside panel, "
"%i in bad regions, %i badrow culled.\n",
nacc, nrej_dis, nrej_hot, nrej_pro, nrej_fra, nrej_bad, ncull);
}
+void search_peaks(struct image *image, float threshold, float min_gradient)
+{
+ int i;
+
+ if ( image->features != NULL ) {
+ image_feature_list_free(image->features);
+ }
+ image->features = image_feature_list_new();
+
+
+ for ( i=0; i<image->det->n_panels; i++ ) {
+
+ struct panel *p = &image->det->panels[i];
+
+ if ( p->no_index ) continue;
+ search_peaks_in_panel(image, threshold, min_gradient, p);
+
+ }
+}
+
+
void dump_peaks(struct image *image, FILE *ofh, pthread_mutex_t *mutex)
{
int i;