aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2021-01-14 17:08:25 +0100
committerThomas White <taw@physics.org>2021-01-14 17:09:35 +0100
commit7fb94888dc081749e004aeef9b2742909e1cbc8e (patch)
treee6be5cd7e784859b5bef86c1c0db3304a32fa7f3
parent0adf5fd72e96e5772288088cac71451a5f898c2e (diff)
Avoid very slow loop over pixels to create bad pixel map
Bad regions specified in terms of x/y still require an iteration over all pixels of the detector, but I don't see an easy way around that. Avoiding x/y bad regions will give best performance.
-rw-r--r--libcrystfel/src/datatemplate.c105
-rw-r--r--libcrystfel/src/datatemplate.h3
-rw-r--r--libcrystfel/src/datatemplate_priv.h5
-rw-r--r--libcrystfel/src/image.c152
-rw-r--r--src/make_pixelmap.c2
5 files changed, 169 insertions, 98 deletions
diff --git a/libcrystfel/src/datatemplate.c b/libcrystfel/src/datatemplate.c
index 4f767c45..d2257f9e 100644
--- a/libcrystfel/src/datatemplate.c
+++ b/libcrystfel/src/datatemplate.c
@@ -103,7 +103,8 @@ static struct dt_badregion *new_bad_region(DataTemplate *det, const char *name)
new->min_ss = 0;
new->max_ss = 0;
new->is_fsss = 99; /* Slightly nasty: means "unassigned" */
- new->panel = NULL;
+ new->panel_name = NULL;
+ new->panel_number = 0; /* Needs to be set after loading */
strcpy(new->name, name);
return new;
@@ -654,7 +655,7 @@ static int parse_field_bad(struct dt_badregion *badr, const char *key,
badr->max_ss = atof(val);
reject = check_badr_fsss(badr, 1);
} else if ( strcmp(key, "panel") == 0 ) {
- badr->panel = strdup(val);
+ badr->panel_name = strdup(val);
} else {
ERROR("Unrecognised field '%s'\n", key);
}
@@ -915,6 +916,30 @@ signed int find_dim(signed int *dims, int which)
}
+static int lookup_panel(const char *panel_name,
+ const DataTemplate *dt,
+ int *res)
+{
+ int i;
+
+ /* If there is exactly one panel, you can get away without
+ * specifying the panel name */
+ if ( (panel_name == NULL) && (dt->n_panels == 1) ) {
+ *res = 0;
+ return 0;
+ }
+
+ for ( i=0; i<dt->n_panels; i++ ) {
+ if ( strcmp(dt->panels[i].name, panel_name) == 0 ) {
+ *res = i;
+ return 0;
+ }
+ }
+
+ return 1;
+}
+
+
DataTemplate *data_template_new_from_string(const char *string_in)
{
DataTemplate *dt;
@@ -1241,17 +1266,31 @@ DataTemplate *data_template_new_from_string(const char *string_in)
}
for ( i=0; i<dt->n_bad; i++ ) {
+
if ( dt->bad[i].is_fsss == 99 ) {
ERROR("Please specify the coordinate ranges for"
" bad region %s\n", dt->bad[i].name);
reject = 1;
}
+
+ if ( dt->bad[i].is_fsss ) {
+ if ( lookup_panel(dt->bad[i].panel_name, dt,
+ &dt->bad[i].panel_number) )
+ {
+ ERROR("No such panel '%s' for bad "
+ "region %s\n",
+ dt->bad[i].panel_name,
+ dt->bad[i].name);
+ reject = 1;
+ }
+ }
}
free(defaults.cnz_from);
free(defaults.data);
free(defaults.mask);
+
for ( rgi=0; rgi<n_rg_definitions; rgi++) {
int pi, n1;
@@ -1495,65 +1534,3 @@ int data_template_get_slab_extents(const DataTemplate *dt,
*ph = h + 1;
return 0;
}
-
-
-/* Return non-zero if pixel fs,ss on panel p is in a bad region
- * as specified in the geometry file (regions only, not including
- * masks, NaN/inf, no_index etc */
-int data_template_in_bad_region(const DataTemplate *dtempl,
- int pn, double fs, double ss)
-{
- double rx, ry;
- double xs, ys;
- int i;
- struct panel_template *p;
-
- if ( pn >= dtempl->n_panels ) {
- ERROR("Panel index out of range\n");
- return 0;
- }
- p = &dtempl->panels[pn];
-
- /* Convert xs and ys, which are in fast scan/slow scan coordinates,
- * to x and y */
- xs = fs*p->fsx + ss*p->ssx;
- ys = fs*p->fsy + ss*p->ssy;
-
- rx = xs + p->cnx;
- ry = ys + p->cny;
-
- for ( i=0; i<dtempl->n_bad; i++ ) {
-
- struct dt_badregion *b = &dtempl->bad[i];
-
- if ( (b->panel != NULL)
- && (strcmp(b->panel, p->name) != 0) ) continue;
-
- if ( b->is_fsss ) {
-
- int nfs, nss;
-
- /* fs/ss bad regions are specified according
- * to the original coordinates */
- nfs = fs + p->orig_min_fs;
- nss = ss + p->orig_min_ss;
-
- if ( nfs < b->min_fs ) continue;
- if ( nfs > b->max_fs ) continue;
- if ( nss < b->min_ss ) continue;
- if ( nss > b->max_ss ) continue;
-
- } else {
-
- if ( rx < b->min_x ) continue;
- if ( rx > b->max_x ) continue;
- if ( ry < b->min_y ) continue;
- if ( ry > b->max_y ) continue;
-
- }
-
- return 1;
- }
-
- return 0;
-}
diff --git a/libcrystfel/src/datatemplate.h b/libcrystfel/src/datatemplate.h
index 66545786..5975f90c 100644
--- a/libcrystfel/src/datatemplate.h
+++ b/libcrystfel/src/datatemplate.h
@@ -87,9 +87,6 @@ extern void data_template_add_copy_header(DataTemplate *dt,
extern int data_template_get_slab_extents(const DataTemplate *dt, int *pw, int *ph);
-extern int data_template_in_bad_region(const DataTemplate *dtempl,
- int pn, double fs, double ss);
-
extern struct rg_collection *data_template_get_rigid_groups(const DataTemplate *dtempl,
const char *collection_name);
diff --git a/libcrystfel/src/datatemplate_priv.h b/libcrystfel/src/datatemplate_priv.h
index c980b53a..f88a6ad0 100644
--- a/libcrystfel/src/datatemplate_priv.h
+++ b/libcrystfel/src/datatemplate_priv.h
@@ -171,14 +171,15 @@ struct dt_badregion
{
char name[1024];
int is_fsss;
- char *panel;
double min_x;
double max_x;
double min_y;
double max_y;
- /* Specified INCLUSIVELY */
+ /* Coordinates are specified INCLUSIVELY */
+ int panel_number;
+ char *panel_name;
int min_fs;
int max_fs;
int min_ss;
diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c
index e8504f23..6da4b435 100644
--- a/libcrystfel/src/image.c
+++ b/libcrystfel/src/image.c
@@ -630,9 +630,67 @@ static void set_image_parameters(struct image *image,
}
-static int flag_value(float pixel, struct panel_template *p)
+static void mark_flagged_pixels_lessthan(float *dp, int *bad,
+ long int n, float val)
{
+ long int i;
+ for ( i=0; i<n; i++ ) {
+ if ( dp[i] < val ) bad[i] = 1;
+ }
+}
+
+
+static void mark_flagged_pixels_morethan(float *dp, int *bad,
+ long int n, float val)
+{
+ long int i;
+ for ( i=0; i<n; i++ ) {
+ if ( dp[i] > val ) bad[i] = 1;
+ }
+}
+
+
+static void mark_flagged_pixels_equal(float *dp, int *bad,
+ long int n, float val)
+{
+ long int i;
+ fenv_t envp;
+
+ fegetenv(&envp);
+ fesetround(1); /* Round to nearest (for flag_value) */
+
+ for ( i=0; i<n; i++ ) {
+ if ( rint(dp[i]) == val ) bad[i] = 1;
+ }
+
+ fesetenv(&envp);
+}
+
+
+static void mark_flagged_pixels_naninf(float *dp, int *bad,
+ long int n)
+{
+ long int i;
+ for ( i=0; i<n; i++ ) {
+ float val = dp[i];
+ if ( isnan(val) || isinf(val) ) bad[i] = 1;
+ }
+}
+
+
+static void mark_flagged_pixels(struct panel_template *p,
+ float *dp, int *bad)
+{
+ int p_w, p_h;
+ long int n;
int i;
+
+ p_w = p->orig_max_fs - p->orig_min_fs + 1;
+ p_h = p->orig_max_ss - p->orig_min_ss + 1;
+ n = p_w * p_h;
+
+ mark_flagged_pixels_naninf(dp, bad, n);
+
for ( i=0; i<MAX_FLAG_VALUES; i++ ) {
float fv = p->flag_values[i];
@@ -643,53 +701,88 @@ static int flag_value(float pixel, struct panel_template *p)
break;
case FLAG_LESSTHAN:
- if ( pixel < fv ) return 1;
+ mark_flagged_pixels_lessthan(dp, bad, n, fv);
break;
case FLAG_MORETHAN:
- if ( pixel > fv ) return 1;
+ mark_flagged_pixels_morethan(dp, bad, n, fv);
break;
case FLAG_EQUAL:
- if ( rint(pixel) == fv) return 1;
+ mark_flagged_pixels_equal(dp, bad, n, fv);
break;
}
}
- return 0;
}
-static void mark_bad_regions(struct panel_template *p,
- float *dp, int *bad,
- const DataTemplate *dtempl,
- int i)
+static void draw_bad_region_fsss(struct dt_badregion *region,
+ int **bad,
+ struct detgeom *detgeom)
{
- int p_w, p_h;
+ struct detgeom_panel *panel;
int fs, ss;
- fenv_t envp;
- fegetenv(&envp);
- fesetround(1); /* Round to nearest (for flag_value) */
+ panel = &detgeom->panels[region->panel_number];
- p_w = p->orig_max_fs - p->orig_min_fs + 1;
- p_h = p->orig_max_ss - p->orig_min_ss + 1;
+ for ( ss=region->min_ss; ss<=region->max_ss; ss++ ) {
+ for ( fs=region->min_fs; fs<=region->max_fs; fs++ ) {
+ bad[region->panel_number][fs+ss*panel->w] = 1;
+ }
+ }
+}
+
+
+static void draw_bad_region_xy(struct dt_badregion *region,
+ int **bad,
+ struct detgeom *detgeom)
+{
+ int i;
+
+ for ( i=0; i<detgeom->n_panels; i++ ) {
+
+ int fs, ss;
+
+ struct detgeom_panel *p = &detgeom->panels[i];
+ for ( ss=0; ss<p->h; ss++ ) {
+ for ( fs=0; fs<p->w; fs++ ) {
+
+ double x, y;
+
+ x = fs*p->fsx + ss*p->ssx + p->cnx;
+ y = fs*p->fsy + ss*p->ssy + p->cny;
+
+ if ( (x > region->min_x )
+ && (x < region->max_x)
+ && (y > region->min_y)
+ && (y < region->max_y) )
+ {
+ bad[i][fs+ss*p->w] = 1;
+ }
- for ( ss=0; ss<p_h; ss++ ) {
- for ( fs=0; fs<p_w; fs++ ) {
- float val = dp[fs+ss*p_w];
- if ( data_template_in_bad_region(dtempl, i,
- fs, ss)
- || isnan(val)
- || isinf(val)
- || flag_value(val, p) )
- {
- bad[fs+ss*p_w] = 1;
}
}
}
+}
- fesetenv(&envp);
+
+static void mark_bad_regions(struct image *image,
+ const DataTemplate *dtempl)
+{
+ int i;
+
+ for ( i=0; i<dtempl->n_bad; i++ ) {
+ if ( dtempl->bad[i].is_fsss ) {
+ draw_bad_region_fsss(&dtempl->bad[i],
+ image->bad,
+ image->detgeom);
+ } else {
+ draw_bad_region_xy(&dtempl->bad[i],
+ image->bad,
+ image->detgeom);
+ }
+ }
}
@@ -756,8 +849,8 @@ static int create_badmap(struct image *image,
/* Add bad regions (skip if panel is bad anyway) */
if ( !p->bad ) {
- mark_bad_regions(p, image->dp[i], image->bad[i],
- dtempl, i);
+ mark_flagged_pixels(p, image->dp[i],
+ image->bad[i]);
}
/* Load mask (skip if panel is bad anyway) */
@@ -775,6 +868,9 @@ static int create_badmap(struct image *image,
dtempl->mask_good, dtempl->mask_bad);
}
}
+
+ mark_bad_regions(image, dtempl);
+
return 0;
}
diff --git a/src/make_pixelmap.c b/src/make_pixelmap.c
index 099ba79f..577c1a1a 100644
--- a/src/make_pixelmap.c
+++ b/src/make_pixelmap.c
@@ -330,7 +330,7 @@ int main(int argc, char *argv[])
y[fs + w*ss] = ry;
z[fs + w*ss] = 0.0; /* 2D part only */
- if ( data_template_in_bad_region(dtempl, pn, cfs, css) ) {
+ if ( image->bad[pn][fs + w*ss] ) {
b[fs + w*ss] = bad_pixel_val;
} else {
b[fs + w*ss] = good_pixel_val;