aboutsummaryrefslogtreecommitdiff
path: root/src/geoptimiser.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/geoptimiser.c')
-rw-r--r--src/geoptimiser.c384
1 files changed, 349 insertions, 35 deletions
diff --git a/src/geoptimiser.c b/src/geoptimiser.c
index f63ba81b..9b03d1b9 100644
--- a/src/geoptimiser.c
+++ b/src/geoptimiser.c
@@ -85,20 +85,25 @@ static void show_help(const char *s)
" -c, --connected=<rg_coll> Rigid group collection for connected\n"
" ASICs.\n"
" --no-error-maps Do not generate error map PNGs.\n"
+" --stretch-map Generate stretch map PNG (panels distance).\n"
" -x, --min-num-peaks-per-pixel=<num> Minimum number of peaks per pixel.\n"
" Default: 3. \n"
+" -z, --max-num-peaks-per-pixel=<num> Maximum number of peaks per pixel.\n"
+" Default is num_patterns / 10. \n"
" -p, --min-num-pixels-per-conn-group=<num> Minimum number of useful pixels per\n"
" connected group.\n"
-" f Default: 100.\n"
+" Default: 100.\n"
" -l, --most-freq-clen Use only the most frequent camera\n"
" length.\n"
" -s, --individual-dist-offset Use a distance offset for each panel.\n"
" Default: whole-detector offset.\n"
" --no-stretch Do not optimize distance offset.\n"
" Default: distance offset is optimized\n"
+" --no-rotation Do not optimize rotation of connectedd groups.\n"
+" Default: rotation is optimized\n"
" -m --max-peak-dist=<num> Maximum distance between predicted and\n"
" detected peaks (in pixels)\n"
-" Default: half of minimal inter-Bragg distance\n"
+" Default: half of minimal inter-Bragg distance\n"
);
}
@@ -108,11 +113,14 @@ struct geoptimiser_params
char *outfile;
char *geometry_filename;
int min_num_peaks_per_pix;
+ int max_num_peaks_per_pix;
int min_num_pix_per_conn_group;
int only_best_distance;
int nostretch;
+ int norotation;
int individual_coffset;
int error_maps;
+ int stretch_map;
int enforce_cspad_layout;
int no_cspad;
double max_peak_dist;
@@ -808,11 +816,13 @@ static int compute_pixel_displacements(struct image *images, int n_images,
static int compute_avg_pix_displ(struct gpanel *gp, int idx,
- int num_peaks_per_pixel)
+ int num_peaks_per_pixel,
+ int max_peaks_per_pixel)
{
int ret;
- if ( gp->num_pix_displ[idx] >= num_peaks_per_pixel ) {
+ if ( gp->num_pix_displ[idx] >= num_peaks_per_pixel &&
+ gp->num_pix_displ[idx] < max_peaks_per_pixel ) {
ret = fill_avg_pixel_displ(gp, idx);
if ( ret != 0 ) return ret;
@@ -833,12 +843,14 @@ static int compute_avg_pix_displ(struct gpanel *gp, int idx,
static int compute_avg_displacements(struct detector *det,
struct rg_collection *connected,
struct connected_data *conn_data,
- struct gpanel *gpanels)
+ struct gpanel *gpanels,
+ struct geoptimiser_params *gparams)
{
int di, ip, ifs, iss;
- int pix_index, ret;
+ int pix_index, ret, max_peaks;
struct panel *p;
+ max_peaks = 0;
for ( di=0; di<connected->n_rigid_groups; di++ ) {
for ( ip=0; ip<connected->rigid_groups[di]->n_panels; ip++ ) {
@@ -854,9 +866,12 @@ static int compute_avg_displacements(struct detector *det,
for ( ifs=0; ifs<p->w; ifs++ ) {
pix_index = ifs+p->w*iss;
+ if ( max_peaks < gp->num_pix_displ[pix_index] )
+ max_peaks = gp->num_pix_displ[pix_index];
ret = compute_avg_pix_displ(gp, pix_index,
- conn_data[di].num_peaks_per_pixel);
+ conn_data[di].num_peaks_per_pixel,
+ gparams->max_num_peaks_per_pix);
if ( ret != 0 ) return ret;
@@ -865,6 +880,7 @@ static int compute_avg_displacements(struct detector *det,
}
}
+ STATUS("Maximum peaks per pixel is: %d\n", max_peaks);
return 0;
}
@@ -1085,20 +1101,23 @@ static void correct_rotation_and_stretch(struct rg_collection *connected,
/* Calculate corner adjustment
* NB All panels follow the first one */
- if ( ip > 0 && connected->rigid_groups[di]->n_panels == 2 && !gparams->no_cspad ) {
+ if ( ip > 0 ) {
struct panel *p0;
- double delta_x, delta_y, delta;
+ double dx_old, dy_old, dx_new, dy_new;
p0 = connected->rigid_groups[di]->panels[0];
- delta_x = p->cnx - p0->cnx / cs;
- delta_y = p->cny - p0->cny / cs;
+ dx_old = p->cnx - p0->cnx / cs;
+ dy_old = p->cny - p0->cny / cs;
- delta = sqrt(delta_x*delta_x + delta_y*delta_y);
+ dx_new = dx_old*cos(conn_data[di].cang)-
+ dy_old*sin(conn_data[di].cang);
+ dy_new = dx_old*sin(conn_data[di].cang)+
+ dy_old*cos(conn_data[di].cang);
- new_cnx = p0->cnx + delta*p0->fsx;
- new_cny = p0->cny + delta*p0->fsy;
+ new_cnx = p0->cnx + dx_new;
+ new_cny = p0->cny + dy_new;
} else {
@@ -1577,6 +1596,14 @@ static double compute_rotation_and_stretch(struct rg_collection *connected,
free(stretches);
}
+ if ( gparams->norotation ) {
+ STATUS("However, rotation will not be optimized, so the "
+ "rotation angle has been set to 0\n");
+ for ( di=0; di<connected->n_rigid_groups; di++ ) {
+ conn_data[di].cang = 0.0;
+ }
+ }
+
stretch_cf = 1.0;
printf("Computing overall stretch coefficient.\n");
@@ -1642,7 +1669,8 @@ static double compute_rotation_and_stretch(struct rg_collection *connected,
#ifdef CAN_SAVE_TO_PNG
static void draw_panel(struct image *image, cairo_t *cr,
- cairo_matrix_t *basic_m, GdkPixbuf **pixbufs, int i) {
+ cairo_matrix_t *basic_m, GdkPixbuf **pixbufs, int i)
+{
struct panel p = image->det->panels[i];
int w = gdk_pixbuf_get_width(pixbufs[i]);
int h = gdk_pixbuf_get_height(pixbufs[i]);
@@ -1816,8 +1844,119 @@ static int save_data_to_png(char *filename, struct detector *det,
return 0;
}
+
+static int save_stretch_to_png(char *filename, struct detector *det,
+ struct rg_collection *connected,
+ struct connected_data *conn_data)
+
+{
+ struct image im;
+ int i, di, ip;
+ float min_str, max_str;
+ struct rectangle rect;
+ GdkPixbuf *col_scale;
+ cairo_t *cr;
+
+ cairo_status_t r;
+ cairo_surface_t *surf;
+
+ im.det = det;
+ im.bad = NULL;
+ im.dp = malloc(det->n_panels*sizeof(float *));
+ if ( im.dp == NULL ) {
+ ERROR("Failed to allocate data\n");
+ return 1;
+ }
+
+ min_str = 10000;
+ max_str = 0;
+ for ( di=0; di<connected->n_rigid_groups; di++ ) {
+ if (conn_data[di].cstr > max_str) max_str = conn_data[di].cstr;
+ if (conn_data[di].cstr < min_str) min_str = conn_data[di].cstr;
+ }
+
+ i = 0;
+ for ( di=0; di<connected->n_rigid_groups; di++ ) {
+ for ( ip=0; ip<connected->rigid_groups[di]->n_panels; ip++ ) {
+
+ int fs, ss;
+ float val;
+ struct panel *p;
+
+ p = connected->rigid_groups[di]->panels[ip];
+
+ im.dp[i] = calloc(p->w * p->h, sizeof(float));
+ if ( im.dp[i] == NULL ) {
+ ERROR("Failed to allocate data\n");
+ return 1;
+ }
+
+ val = (conn_data[di].cstr - min_str) / (max_str - min_str);
+ val *= 10.0; /* render_panels sets this as max */
+
+ for ( ss=0; ss<p->h; ss++ ) {
+ for ( fs=0; fs<p->w; fs++ ) {
+
+ im.dp[i][fs+p->w*ss] = val;
+
+ }
+ }
+ i++;
+ }
+ }
+
+ get_pixel_extents(im.det, &rect.min_x, &rect.min_y, &rect.max_x,
+ &rect.max_y);
+
+ if (rect.min_x > 0.0) rect.min_x = 0.0;
+ if (rect.max_x < 0.0) rect.max_x = 0.0;
+ if (rect.min_y > 0.0) rect.min_y = 0.0;
+ if (rect.max_y < 0.0) rect.max_y = 0.0;
+
+ rect.width = rect.max_x - rect.min_x;
+ rect.height = rect.max_y - rect.min_y;
+
+ /* Add a thin border */
+ rect.width += 2.0;
+ rect.height += 2.0;
+ surf = cairo_image_surface_create(CAIRO_FORMAT_ARGB32, rect.width + 20,
+ rect.height);
+
+ draw_detector(surf, &im, rect);
+
+ col_scale = render_get_colour_scale(20, rect.height, SCALE_GEOPTIMISER);
+
+ cr = cairo_create(surf);
+ cairo_identity_matrix(cr);
+ cairo_translate(cr, rect.width, 0.0);
+ cairo_rectangle(cr, 0.0, 0.0, 20.0, rect.height);
+ gdk_cairo_set_source_pixbuf(cr, col_scale, 0.0, 0.0);
+ cairo_fill(cr);
+ cairo_destroy(cr);
+
+ for ( i=0; i<det->n_panels; i++ ) {
+ free(im.dp[i]);
+ }
+ free(im.dp);
+
+ r = cairo_surface_write_to_png(surf, filename);
+ if ( r != CAIRO_STATUS_SUCCESS ) return 1;
+
+ return 0;
+}
+
#else /* CAN_SAVE_TO_PNG */
+static int save_stretch_to_png(char *filename, struct detector *det,
+ struct rg_collection *connected,
+ struct connected_data *conn_data)
+{
+ ERROR("WARNING: geoptimiser was compiled without gdk-pixbuf and cairo "
+ "support. Error maps will not be saved.\n");
+ return 0; /* Return "success" so that program continues */
+}
+
+
static int save_data_to_png(char *filename, struct detector *det,
struct gpanel *gpanels)
{
@@ -1902,6 +2041,87 @@ int check_and_enforce_cspad_dist(struct geoptimiser_params *gparams,
}
}
+ return num_errors_found;
+}
+
+
+int check_and_enforce_agipd_dist(struct geoptimiser_params *gparams,
+ struct detector *det)
+{
+ int np = 0;
+ int npp = 0;
+ int num_errors_found = 0;
+
+ double dist_to_check = 66.0;
+ double tol = 0.1;
+
+ for ( np=0; np<det->n_panels; np = np+8 ) {
+
+ double dist2;
+ double cnx2, cny2;
+
+ struct panel *ep = &det->panels[np];
+
+ for ( npp=0; npp<8; npp++ ) {
+
+ struct panel *op = &det->panels[np+npp];
+
+ cnx2 = ep->cnx + npp*dist_to_check*ep->ssx;
+ cny2 = ep->cny + npp*dist_to_check*ep->ssy;
+
+ dist2 = (( cnx2 - op->cnx )*( cnx2 - op->cnx ) +
+ ( cny2 - op->cny )*( cny2 - op->cny ));
+
+ if ( dist2 > (tol*tol)) {
+
+ num_errors_found += 1;
+
+ STATUS("Warning: distance between panels %s and %s "
+ "is outside acceptable margins (Corners are "
+ "more than 0.2 pixels away: %3.2f).\n", ep->name,
+ op->name, sqrt(dist2));
+
+ if ( gparams->enforce_cspad_layout ) {
+
+ double new_op_cx, new_op_cy;
+
+ new_op_cx = ep->cnx + ep->ssx*npp*dist_to_check;
+ new_op_cy = ep->cny + ep->ssy*npp*dist_to_check;
+
+ op->cnx = new_op_cx;
+ op->cny = new_op_cy;
+
+ STATUS("Enforcing distance....\n");
+ }
+
+ }
+
+ if ( ep->fsx != op->fsx || ep->ssx != op->ssx ||
+ ep->fsy != op->fsy || ep->ssx != op->ssx ) {
+
+ num_errors_found += 1;
+
+ STATUS("Warning: relative orientation between panels "
+ "%s and %s is incorrect.\n", ep->name, op->name);
+
+ if ( gparams->enforce_cspad_layout ) {
+
+ STATUS("Enforcing relative orientation....\n");
+
+ op->fsx = ep->fsx;
+ op->ssx = ep->ssx;
+ op->fsy = ep->fsy;
+ op->ssy = ep->ssy;
+
+ op->xfs = ep->xfs;
+ op->xss = ep->xss;
+ op->yfs = ep->yfs;
+ op->yss = ep->yss;
+ }
+
+ }
+ }
+ }
return num_errors_found;
}
@@ -2076,9 +2296,9 @@ int optimize_geometry(struct geoptimiser_params *gparams,
double avg_res;
double clen_to_use;
- // for angles and stretch calculation use
- // only pixels which are distco*size_panel
- // away
+ /* for angles and stretch calculation use
+ * only pixels which are distco*size_panel
+ * away */
double dist_coeff_for_rot_str = 0.2;
double total_error;
@@ -2099,6 +2319,7 @@ int optimize_geometry(struct geoptimiser_params *gparams,
"accurate estimation of position/orientation: %i\n",
gparams->min_num_pix_per_conn_group);
+ /* CS-PAD */
if ( (det->n_panels == 64) && !gparams->no_cspad ) {
int num_errors = 0;
@@ -2142,6 +2363,51 @@ int optimize_geometry(struct geoptimiser_params *gparams,
}
}
+ /* AGIPD */
+ if ( (det->n_panels == 128) && det->panels[0].res > 4999 &&
+ det->panels[0].res < 5001 && !gparams->no_cspad ) {
+
+ int num_errors = 0;
+
+ STATUS("It looks like the detector is an AGIPD. "
+ "Checking relative distance and orientation of "
+ "connected ASICS.\n");
+ STATUS("If the detector is not an AGIPD , please rerun the "
+ "program with the --no-agipd option.\n");
+
+ STATUS("Enforcing AGIPD layout...\n");
+ num_errors = check_and_enforce_agipd_dist(gparams, det);
+
+ if ( gparams->enforce_cspad_layout ) {
+
+ int geom_wr;
+
+ STATUS("Saving geometry with enforced AGIPD layout.\n"
+ "Please restart geometry optimization using the "
+ "optimized geometry from this run as input "
+ "geometry file.\n");
+ geom_wr = write_detector_geometry_2(
+ gparams->geometry_filename,
+ gparams->outfile, det,
+ gparams->command_line, 1);
+ if ( geom_wr != 0 ) {
+ ERROR("Error in writing output geometry file.\n");
+ return 1;
+ }
+ STATUS("All done!\n");
+ return 0;
+ }
+
+ if ( !gparams->enforce_cspad_layout && num_errors > 0 ) {
+ ERROR("Relative distance and orientation of connected "
+ "ASICS do not respect the AGIPD layout.\n"
+ "Geometry optimization cannot continue.\n"
+ "Please rerun the program with the "
+ "--enforce-agipd-layout option.\n");
+ return 1;
+ }
+ }
+
images = read_patterns_from_stream(gparams->infile, det, &n_images);
if ( (n_images < 1) || (images == NULL) ) {
ERROR("Error reading stream file\n");
@@ -2164,6 +2430,14 @@ int optimize_geometry(struct geoptimiser_params *gparams,
avg_cell);
if ( clen_to_use < 0.0 ) return 1;
+ if ( gparams->max_num_peaks_per_pix == 0 && n_images > 100 ) {
+ gparams->max_num_peaks_per_pix = n_images / 10;
+ } else if ( gparams->max_num_peaks_per_pix == 0 ) {
+ gparams->max_num_peaks_per_pix = n_images;
+ }
+ STATUS("Maximum number of measurements for a pixel to be included in "
+ "the refinement: %i\n", gparams->max_num_peaks_per_pix);
+
gpanels = calloc(det->n_panels, sizeof(struct gpanel));
if ( gpanels == NULL ) {
ERROR("Failed to allocate panels.\n");
@@ -2206,7 +2480,7 @@ int optimize_geometry(struct geoptimiser_params *gparams,
adjust_min_peaks_per_conn(connected, gpanels, det, gparams, conn_data);
- if ( compute_avg_displacements(det, connected, conn_data, gpanels) ) {
+ if ( compute_avg_displacements(det, connected, conn_data, gpanels, gparams) ) {
free(conn_data);
free(images);
return 1;
@@ -2235,25 +2509,53 @@ int optimize_geometry(struct geoptimiser_params *gparams,
}
- stretch_coeff = compute_rotation_and_stretch(connected, conn_data,
- det, gpanels,
+ if ( !gparams->nostretch || !gparams->norotation ) {
+
+ stretch_coeff = compute_rotation_and_stretch(connected,
+ conn_data, det, gpanels,
dist_coeff_for_rot_str,
gparams);
- if ( stretch_coeff < 0.0 ) {
- free(conn_data);
- return 1;
- }
+ if ( stretch_coeff < 0.0 ) {
+ free(conn_data);
+ return 1;
+ }
- ret = compute_rot_stretch_for_empty_panels(quadrants, connected,
- gparams->min_num_pix_per_conn_group, conn_data);
- if ( ret ) {
- free(conn_data);
- return 1;
- }
+ ret = compute_rot_stretch_for_empty_panels(quadrants, connected,
+ gparams->min_num_pix_per_conn_group,
+ conn_data);
+ if ( ret ) {
+ free(conn_data);
+ return 1;
+ }
+
+ if ( gparams->stretch_map ) {
+
+#ifdef HAVE_SAVE_TO_PNG
+
+ STATUS("Saving stretch map - for out-of-plane rotations.\n");
- correct_rotation_and_stretch(connected, det, gpanels, conn_data,
- clen_to_use, stretch_coeff,
- gparams);
+ ret = save_stretch_to_png("stretch_co.png", det,
+ connected, conn_data);
+
+ if ( ret ) {
+ ERROR("Error while writing data to file.\n");
+ free(conn_data);
+ return 1;
+ }
+
+#else /* HAVE_SAVE_TO_PNG */
+
+ STATUS("ERROR: geoptimiser was compiled without GTK and "
+ "cairo support.\n Stretch map will not be saved.\n");
+
+#endif /* HAVE_SAVE_TO_PNG */
+
+ }
+
+ correct_rotation_and_stretch(connected, det, gpanels,
+ conn_data, clen_to_use,
+ stretch_coeff, gparams);
+ }
ret = compute_shift(connected, conn_data, det, gparams, gpanels);
if ( ret != 0 ) {
@@ -2341,13 +2643,16 @@ int main(int argc, char *argv[])
gparams->infile = NULL;
gparams->geometry_filename = NULL;
gparams->min_num_peaks_per_pix = 3;
+ gparams->max_num_peaks_per_pix = 0;
gparams->min_num_pix_per_conn_group = 100;
gparams->only_best_distance = 0;
gparams->enforce_cspad_layout = 0;
gparams->nostretch = 0;
+ gparams->norotation = 0;
gparams->individual_coffset = 0;
gparams->no_cspad = 0;
gparams->error_maps = 1;
+ gparams->stretch_map = 0;
gparams->max_peak_dist = 0.0;
const struct option longopts[] = {
@@ -2361,6 +2666,7 @@ int main(int argc, char *argv[])
{"quadrants", 1, NULL, 'q'},
{"connected", 1, NULL, 'c'},
{"min-num-peaks-per-pixel", 1, NULL, 'x'},
+ {"max-num-peaks-per-pixel", 0, NULL, 'z'},
{"min-num-pixels-per-conn-group", 1, NULL, 'p'},
{"most-few-clen", 0, NULL, 'l'},
{"max-peak-dist", 1, NULL, 'm'},
@@ -2368,9 +2674,13 @@ int main(int argc, char *argv[])
/* Long-only options with no arguments */
{"no-stretch", 0, &gparams->nostretch, 1},
+ {"no-rotation", 0, &gparams->norotation, 1},
{"no-error-maps", 0, &gparams->error_maps, 0},
+ {"stretch-map", 0, &gparams->stretch_map, 1},
{"enforce-cspad-layout", 0, &gparams->enforce_cspad_layout, 1},
+ {"enforce-agipd-layout", 0, &gparams->enforce_cspad_layout, 1},
{"no-cspad", 0, &gparams->no_cspad, 1},
+ {"no-agipd", 0, &gparams->no_cspad, 1},
/* Long-only options with arguments */
{"min-num-peaks-per-panel", 1, NULL, 11},
@@ -2380,7 +2690,7 @@ int main(int argc, char *argv[])
};
/* Short options */
- while ((c = getopt_long(argc, argv, "ho:i:g:q:c:o:x:p:lsm:",
+ while ((c = getopt_long(argc, argv, "ho:i:g:q:c:o:x:z:p:lsm:",
longopts, NULL)) != -1) {
switch (c) {
@@ -2425,6 +2735,10 @@ int main(int argc, char *argv[])
gparams->min_num_peaks_per_pix = atoi(optarg);
break;
+ case 'z' :
+ gparams->max_num_peaks_per_pix = atoi(optarg);
+ break;
+
case 'p' :
gparams->min_num_pix_per_conn_group = atoi(optarg);
break;
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-30 18:09:56 +0000