diff options
Diffstat (limited to 'src/geoptimiser.c')
-rw-r--r-- | src/geoptimiser.c | 384 |
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; |