From 36fe26d5c45c8205f44617960cec3705bf40a9f2 Mon Sep 17 00:00:00 2001 From: Valerio Mariani Date: Fri, 14 Aug 2015 17:36:35 +0200 Subject: Major refactoring of geoptimiser --- src/geoptimiser.c | 2750 ++++++++++++++++++++++++++++------------------------- 1 file changed, 1462 insertions(+), 1288 deletions(-) (limited to 'src/geoptimiser.c') diff --git a/src/geoptimiser.c b/src/geoptimiser.c index ba38f63c..e1d59f05 100644 --- a/src/geoptimiser.c +++ b/src/geoptimiser.c @@ -68,34 +68,54 @@ static void show_help(const char *s) printf( "Refines detector geometry.\n" "\n" -" -h, --help Display this help message.\n" +" -h, --help Display this help message.\n" "\n" -" --version Print CrystFEL version number and\n" -" exit.\n" -" -i, --input= Specify stream file to be used for \n" -" geometry optimization.\n" -" -g. --geometry= Input detector geometry.\n" -" -o, --output= Output detector geometry.\n" -" -q, --quadrants= Rigid group collection for quadrants.\n" -" -c, --connected= Rigid group collection for connected\n" -" ASICs.\n" -" --no-error-maps Do not generate error map PNGs.\n" -" -x, --min-num-peaks-per-pixel= Minimum number of peaks per pixel.\n" -" Default: 3. \n" -" -p, --min-num-peaks-per-panel= Minimum number of peaks per pixel.\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" -" -m --max-peak-dist= Maximum distance between predicted and\n" -" detected peaks (in pixels)\n" -" Default: 4.0 pixels.\n" +" --version Print CrystFEL version number and\n" +" exit.\n" +" -i, --input= Specify stream file to be used for \n" +" geometry optimization.\n" +" -g. --geometry= Get detector geometry from file.\n" +" -o, --output= Output stream.\n" +" -q, --quadrants= Rigid group collection for quadrants.\n" +" -c, --connected= Rigid group collection for connected\n" +" ASICs.\n" +" --no-error-maps Do not generate error map PNGs.\n" +" -x, --min-num-peaks-per-pixel= Minimum number of peaks per pixel.\n" +" Default: 3. \n" +" --min-num-peaks-per-panel= DEPRECATED. This option has been\n" +" renamed to --min-num-pixels-per-conn-group.\n" +" -p, --min-num-pixels-per-conn-group= Minimum number of useful pixels pern" +" connected group.\n" +" f 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" +" -m --max-peak-dist= Maximum distance between predicted and\n" +" detected peaks (in pixels)\n" +" Default: 4.0 pixels.\n" ); } +struct geoptimiser_params +{ + char *infile; + char *outfile; + char *geometry_filename; + int min_num_peaks_per_pix; + int min_num_pix_per_conn_group; + int only_best_distance; + int nostretch; + int individual_coffset; + int error_maps; + int enforce_cspad_layout; + int no_cspad; + double max_peak_dist; + const char *command_line; +}; + struct connected_data { @@ -127,19 +147,76 @@ struct pattern_list { }; -struct single_pix_displ +struct single_pixel_displ { double dfs; double dss; - struct single_pix_displ* ne; + struct single_pixel_displ* ne; +}; + + +struct pixel_displ_list +{ + struct single_pixel_displ *pix_displ_list; + struct single_pixel_displ **curr_pix_displ; + int *num_pix_displ; }; struct connected_stretch_and_angles { double *stretch_coeff; - unsigned int *num_angles; - int num_coeff; + long *num_angles; + long num_coeff; +}; + +struct cell_params +{ + double a; + double b; + double c; + double alpha; + double beta; + double gamma; +}; + + +struct avg_displacements +{ + double *displ_x; + double *displ_y; + double *displ_abs; +}; + + +struct avg_rot_and_stretch +{ + double *aver_ang; + double *aver_str; + int *aver_num_ang; +}; + + +struct avg_shift +{ + double *aver_x; + double *aver_y; + int *aver_num_sh; +}; + +struct pixel_maps +{ + double *pix_to_x; + double *pix_to_y; +}; + + +struct enhanced_det +{ + struct detector *det; + int width; + int height; + int num_pix; }; @@ -176,9 +253,8 @@ static Reflection *find_closest_reflection(RefList *rlist, for ( refl = first_refl(rlist, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { + refl != NULL; + refl = next_refl(refl, iter) ) { double ds; double rfs, rss; double rx, ry; @@ -248,7 +324,7 @@ static double extract_f_from_stuff(const char *field_name, if ( strncmp(stuff->fields[i], field_name_plus_equal, strlen(field_name_plus_equal)) == 0 ) { - return atoi(stuff->fields[i]+ + return atof(stuff->fields[i]+ strlen(field_name_plus_equal)); } } @@ -261,7 +337,6 @@ static double extract_f_from_stuff(const char *field_name, static struct pattern_list *read_patterns_from_steam_file(const char *infile, struct detector *det) { - Stream *st; struct pattern_list *pattern_list; @@ -390,7 +465,7 @@ static struct pattern_list *read_patterns_from_steam_file(const char *infile, patn->unit_cells = new_unit_cells; crystal_reflist = crystal_get_reflections( - cur.crystals[i]); + cur.crystals[i]); for ( refl = first_refl(crystal_reflist, &iter); refl != NULL; @@ -402,7 +477,6 @@ static struct pattern_list *read_patterns_from_steam_file(const char *infile, get_indices(refl, &h, &k, &l); n = add_refl(patn->ref_list, h, k, l); copy_data(n, refl); - } } @@ -445,13 +519,14 @@ static struct rvec get_q_from_xyz(double rx, double ry, double dist, double l) } -static void compute_avg_cell_parameters(struct pattern_list *pattern_list, - double *avcp) +static struct cell_params *compute_avg_cell_parameters(struct pattern_list + *pattern_list) { int numavc; int j, i; double minc[6]; double maxc[6]; + double avg_cpar[6]; for (j=0; j<6; j++) { minc[j] = 1e10; @@ -484,7 +559,7 @@ static void compute_avg_cell_parameters(struct pattern_list *pattern_list, cpar[5] = rad2deg(cpar[5]); for ( j=0; j<6; j++ ) { - avcp[j] += cpar[j]; + avg_cpar[j] += cpar[j]; if (cpar[j]maxc[j]) maxc[j] = cpar[j]; } @@ -495,33 +570,41 @@ static void compute_avg_cell_parameters(struct pattern_list *pattern_list, } if ( numavc>0 ) { - for ( j=0; j<6; j++ ) avcp[j] /= numavc; + for ( j=0; j<6; j++ ) avg_cpar[j] /= numavc; } + struct cell_params *cparams = malloc(sizeof(struct cell_params)); + + cparams->a = avg_cpar[0]; + cparams->b = avg_cpar[1]; + cparams->c = avg_cpar[2]; + cparams->alpha = avg_cpar[3]; + cparams->beta = avg_cpar[4]; + cparams->gamma = avg_cpar[5]; + STATUS("Average cell coordinates:\n"); STATUS("Average a, b, c (in nm): %6.3f, %6.3f, %6.3f\n", - avcp[0],avcp[1],avcp[2]); + cparams->a, cparams->b, cparams->c); STATUS("Minimum -Maximum a, b, c:\n" "\t%6.3f - %6.3f,\n" "\t%6.3f - %6.3f,\n" "\t%6.3f - %6.3f\n", minc[0], maxc[0], minc[1], maxc[1], minc[2], maxc[2]); STATUS("Average alpha,beta,gamma in degrees: %6.3f, %6.3f, %6.3f\n", - avcp[3], avcp[4], avcp[5]); + cparams->alpha, cparams->beta, cparams->gamma); STATUS("Minimum - Maximum alpha,beta,gamma in degrees:\n" "\t%5.2f - %5.2f,\n" "\t%5.2f - %5.2f,\n" "\t%5.2f - %5.2f\n", minc[3], maxc[3], minc[4], maxc[4], minc[5], maxc[5]); + return cparams; } -static double compute_clen_to_use(struct pattern_list *pattern_list, - double istep, double *avcp, - double max_peak_distance, - int only_best_distance, - double* min_braggp_dist) +static double pick_clen_to_use(struct geoptimiser_params *gparams, + struct pattern_list *pattern_list, + double avg_res, struct cell_params *cparams) { int cp, i, u; int num_clens; @@ -530,9 +613,8 @@ static double compute_clen_to_use(struct pattern_list *pattern_list, int *clens_population; double *clens; double *lambdas; - double irecistep; + double min_braggp_dist; double clen_to_use; - double *braggp_dist; struct rvec cqu; max_clens = 1024; @@ -586,16 +668,16 @@ static double compute_clen_to_use(struct pattern_list *pattern_list, clens_new = realloc(clens_population, (max_clens+1024)*sizeof(double)); lambdas_new = realloc(clens_population, - (max_clens+1024)*sizeof(double)); + (max_clens+1024)*sizeof(double)); if ( clens_new == NULL || clens_population_new == NULL || lambdas_new == NULL) { ERROR("Failed to allocate memory for " "camera length list\n"); - free(clens); - free(clens_population); - free(lambdas); - return -1.0; + free(clens); + free(clens_population); + free(lambdas); + return -1.0; } max_clens += 1024; @@ -605,7 +687,6 @@ static double compute_clen_to_use(struct pattern_list *pattern_list, } clens[num_clens] = pattern_list->patterns[cp]->clen; - clens_population[num_clens] = 1; lambdas[num_clens] = pattern_list->patterns[cp]->lambda; num_clens++; @@ -618,37 +699,34 @@ static double compute_clen_to_use(struct pattern_list *pattern_list, } if ( num_clens == 1 ) { - STATUS("All patterns in the stream file have the same camera " - "length: %f m.\n", + STATUS("All patterns have the same camera length: %f m.\n", clens[0]); } else { STATUS("%i different camera lengths were found for the input " - "patterns in the stream file:\n", num_clens); + "patterns:\n", num_clens); } - braggp_dist = calloc(num_clens, sizeof(double)); - *min_braggp_dist = 999999.0; best_clen = 0; clen_to_use = clens[0]; for ( i=0; i0) { - cqu = get_q_from_xyz(1/istep, 0, clens[i], lambdas[i]); - - irecistep = 1/cqu.u; - - braggp_dist[i] = fmin(fmin(irecistep/avcp[0], - irecistep/avcp[1]), - irecistep/avcp[2]); - if (braggp_dist[i] < *min_braggp_dist) { - *min_braggp_dist = braggp_dist[i]; - } - STATUS("Camera length %0.4f m was found %i times in the " - "stream file.\n" - "Corresponding minimum inter-bragg peak distance " - "is (based on average cell parameters): " - " %0.1f pixels.\n", + cqu = get_q_from_xyz(1.0/avg_res, 0, clens[i], lambdas[i]); + + min_braggp_dist = fmin(fmin((1.0/cqu.u)/cparams->a, + (1.0/cqu.u)/cparams->b), + (1.0/cqu.u)/cparams->c); + STATUS("Camera length %0.4f m was found %i times.\n" + "Minimum inter-bragg peak distance (based on " + "average cell parameters): %0.1f pixels.\n", clens[i], clens_population[i], - braggp_dist[i]); + min_braggp_dist); + if ( min_braggp_dist<1.2*gparams->max_peak_dist ) { + STATUS("WARNING: The distance between Bragg " + "peaks is too small: " + "%0.1f < 1.2*%0.1f pixels.\n", + min_braggp_dist, + gparams->max_peak_dist); + } if ( clens_population[i] > clens_population[best_clen] ) { best_clen = i; clen_to_use = clens[best_clen]; @@ -656,24 +734,23 @@ static double compute_clen_to_use(struct pattern_list *pattern_list, } } - if ( only_best_distance ) { - STATUS("Only the %i patterns with clen=%0.4f m will be used.\n", + if ( gparams->only_best_distance ) { + STATUS("Only %i patterns with CLEN=%0.4f m will be used.\n", clens_population[best_clen], clen_to_use); - *min_braggp_dist = braggp_dist[best_clen]; } free(clens); free(lambdas); free(clens_population); - free(braggp_dist); + return clen_to_use; } -static double comp_median(double *arr, unsigned int n) +static double comp_median(double *arr, long n) { - int low, high, median, middle, ll, hh; + long low, high, median, middle, ll, hh; double A; if (n<1) return 0.0; @@ -739,7 +816,7 @@ static double comp_median(double *arr, unsigned int n) /* Re-set active partition */ if ( hh<=median ) low = ll; - if ( hh>=median ) high = hh-1; + if ( hh>=median ) high = hh-1; } return 0.0; @@ -767,199 +844,211 @@ static int find_quad_for_connected(struct rigid_group *rg, } -static void free_all_curr_pix_displ(struct single_pix_displ *all_pix_displ, - struct single_pix_displ **curr_pix_displ, - int num_pix_in_slab) +static int fill_avg_pixel_displ(struct pixel_displ_list *pix_displ_list, + int pix_index, + struct avg_displacements *avg_displ) { - int i; - struct single_pix_displ *curr = NULL; - struct single_pix_displ *next = NULL; + double *list_fs_displ; + double *list_ss_displ; + int count = 0; + int ei; + + list_fs_displ = calloc(pix_displ_list->num_pix_displ[pix_index], + sizeof(double)); + if ( list_fs_displ == NULL ) { + ERROR("Failed to allocate memory for pixel statistics.\n"); + return 1; + } + list_ss_displ = calloc(pix_displ_list->num_pix_displ[pix_index], + sizeof(double)); + if ( list_ss_displ == NULL ) { + ERROR("Failed to allocate memory for pixel statistics.\n"); + free(list_fs_displ); + return 1; + } + + pix_displ_list->curr_pix_displ[pix_index] = + &pix_displ_list->pix_displ_list[pix_index]; - for ( i=0; inum_pix_displ[pix_index]; ei++ ) { - curr = &all_pix_displ[i]; + struct single_pixel_displ * pixel_displ_entry; + pixel_displ_entry = pix_displ_list->curr_pix_displ[pix_index]; + if ( pixel_displ_entry->dfs == -10000.0 ) break; + list_fs_displ[count] = pixel_displ_entry->dfs; + list_ss_displ[count] = pixel_displ_entry->dss; + count++; + if ( pixel_displ_entry->ne == NULL ) { + break; + } else { + pix_displ_list->curr_pix_displ[pix_index] = + pix_displ_list->curr_pix_displ[pix_index]->ne; + } + } - if ( curr->ne != NULL ) { - curr = curr->ne; - while ( curr != NULL ) { - next = curr->ne; - free(curr); - curr = next; - } - } - } + if ( count < 1 ) { + free(list_fs_displ); + free(list_ss_displ); + return 0; + } + + avg_displ->displ_x[pix_index] = comp_median(list_fs_displ, count); + avg_displ->displ_y[pix_index] = comp_median(list_ss_displ, count); + avg_displ->displ_abs[pix_index] = modulus2d(avg_displ->displ_x[pix_index], + avg_displ->displ_y[pix_index]); + free(list_fs_displ); + free(list_ss_displ); - free(curr_pix_displ); - free(all_pix_displ); + return 0; } -static int fill_pixel_statistics(int *num_pix_displ, - struct single_pix_displ** curr_pix_displ, - struct single_pix_displ* all_pix_displ, - struct connected_data *conn_data, - int ifs, int iss, int di, int aw, int dfv, - double *displ_x, - double *displ_y, double *displ_abs) +static int allocate_next_element(struct single_pixel_displ** curr_pix_displ, + int pix_index) { - double *cPxAfs; - double *cPxAss; - int cnu = 0; - - cPxAfs = calloc(num_pix_displ[ifs+aw*iss], sizeof(double)); - if ( cPxAfs == NULL ) { + curr_pix_displ[pix_index]->ne = malloc(sizeof(struct single_pixel_displ)); + if ( curr_pix_displ[pix_index]->ne == NULL ) { ERROR("Failed to allocate memory for pixel statistics.\n"); return 1; } - cPxAss = calloc(num_pix_displ[ifs+aw*iss], sizeof(double)); - if ( cPxAss == NULL ) { - ERROR("Failed to allocate memory for pixel statistics.\n"); - free(cPxAfs); - return 1; - } - curr_pix_displ[ifs+aw*iss] = &all_pix_displ[ifs+aw*iss]; + curr_pix_displ[pix_index] = curr_pix_displ[pix_index]->ne; + + return 0; +} - while (1) { - if (curr_pix_displ[ifs+aw*iss]->dfs == dfv) break; - cPxAfs[cnu] = curr_pix_displ[ifs+aw*iss]->dfs; - cPxAss[cnu] = curr_pix_displ[ifs+aw*iss]->dss; - cnu++; - if ( curr_pix_displ[ifs+aw*iss]->ne == NULL ) break; - curr_pix_displ[ifs+aw*iss] = curr_pix_displ[ifs+aw*iss]->ne; - } - if ( cnu < 1 ) return 0; +static int add_distance_to_list(struct enhanced_det *edet, + struct imagefeature *imfe, + struct pixel_displ_list *pix_displ_list, + Reflection *refl, double fx, double fy) +{ + int pix_index; + double rfs, rss; + double crx, cry; + + pix_index = ((int)rint(imfe->fs) + edet->width*(int)rint(imfe->ss)); + + if ( pix_displ_list->num_pix_displ[pix_index]>0 ) { + + int ret; - displ_x[ifs+aw*iss] = comp_median(cPxAfs, cnu); - displ_y[ifs+aw*iss] = comp_median(cPxAss, cnu); - displ_abs[ifs+aw*iss] = modulus2d(displ_x[ifs+aw*iss], - displ_y[ifs+aw*iss]); - conn_data[di].n_peaks_in_conn++; + ret = allocate_next_element(pix_displ_list->curr_pix_displ, + pix_index); - free(cPxAfs); - free(cPxAss); + if ( ret != 0) return ret; + + } + + get_detector_pos(refl, &rfs, &rss); + compute_x_y(edet->det, rfs, rss, &crx, &cry); + pix_displ_list->curr_pix_displ[pix_index]->dfs = (fx-crx); + pix_displ_list->curr_pix_displ[pix_index]->dss = (fy-cry); + pix_displ_list->curr_pix_displ[pix_index]->ne = NULL; + pix_displ_list->num_pix_displ[pix_index]++; return 0; } - -static int compute_panel_statistics(struct rg_collection *connected, - int *num_pix_displ, - struct single_pix_displ** curr_pix_displ, - struct single_pix_displ* all_pix_displ, - struct connected_data *conn_data, - int di, int ip, int np, - double dfv, - int aw, double *displ_x, - double *displ_y, double *displ_abs) +static int count_pixels_with_min_peaks( struct panel *p, int min_num_peaks, + struct pixel_displ_list *pix_displ_list, + struct enhanced_det *edet) { - struct panel *p; + int pix_index; + int pixel_count = 0; int ifs, iss; - p = connected->rigid_groups[di]->panels[ip]; - for ( ifs=p->min_fs; ifsmax_fs+1; ifs++ ) { for ( iss=p->min_ss; issmax_ss+1; iss++ ) { - if ( num_pix_displ[ifs+aw*iss]>=np ) { - - int ret; - - ret = fill_pixel_statistics(num_pix_displ, - curr_pix_displ, - all_pix_displ, - conn_data, - ifs, iss, di, aw, - dfv, displ_x, - displ_y, displ_abs); - - if ( ret == -2 ) { - break; - } else if ( ret != 0 ) { - return ret; - } - } else { - displ_x[ifs+aw*iss] = dfv; - displ_y[ifs+aw*iss] = dfv; - displ_abs[ifs+aw*iss] = dfv; + pix_index = ifs+edet->width*iss; + + if ( pix_displ_list->num_pix_displ[pix_index] >= + min_num_peaks ) { + pixel_count += 1; } } } - return 0; + return pixel_count; } - -static int allocate_next_element(struct single_pix_displ** curr_pix_displ, - int ipx) +static void adjust_min_peaks_per_conn(struct enhanced_det *edet, + struct rg_collection *connected, + struct geoptimiser_params *gparams, + struct connected_data *conn_data, + struct pixel_displ_list *pix_displ_list) { - curr_pix_displ[ipx]->ne = malloc(sizeof(struct single_pix_displ)); - if ( curr_pix_displ[ipx]->ne == NULL ) { - ERROR("Failed to allocate memory for pixel statistics.\n"); - return 1; - } - curr_pix_displ[ipx] = curr_pix_displ[ipx]->ne; + int min_num_peaks, di, ip; + struct panel *p; - return 0; -} + STATUS("Adjusting the minimum number of measurements per pixel in " + "order to have enough measurements for each connected group.\n"); + for ( di=0; din_rigid_groups; di++ ) { + + for ( min_num_peaks=gparams->min_num_peaks_per_pix; + min_num_peaks>0; min_num_peaks-- ) { -static int compute_pixel_statistics(struct pattern_list *pattern_list, - struct detector *det, - struct rg_collection *connected, - struct rg_collection *quadrants, - int num_pix_in_slab, - int max_peak_distance, int aw, - double dfv, - double min_num_peaks_per_pixel, - double min_num_peaks_per_panel, - int only_best_distance, - double clen_to_use, - double *slab_to_x, double *slab_to_y, - struct connected_data *conn_data, - double *displ_x, - double *displ_y, double *displ_abs, - struct single_pix_displ* all_pix_displ, - struct single_pix_displ** curr_pix_displ, - int *num_pix_displ) -{ - int ipx, cp, ich, di, ip, np; + int di_count = 0; - for (di=0; din_rigid_groups; di++) { + for (ip=0; iprigid_groups[di]->n_panels; + ip++) { + + int pix_count; + + p = connected->rigid_groups[di]->panels[ip]; - conn_data[di].num_quad = find_quad_for_connected( - connected->rigid_groups[di], - quadrants); - conn_data[di].cang = 0.0; - conn_data[di].cstr = 1.0; - conn_data[di].sh_x = dfv; - conn_data[di].sh_y = dfv; - conn_data[di].num_peaks_per_pixel = 1; - conn_data[di].name = connected->rigid_groups[di]->name; - conn_data[di].n_peaks_in_conn = 0; - } + pix_count = count_pixels_with_min_peaks(p, + min_num_peaks, + pix_displ_list, + edet); + + di_count += pix_count; + } + conn_data[di].n_peaks_in_conn = di_count; + if ( di_count >= + gparams->min_num_pix_per_conn_group ) { + conn_data[di].num_peaks_per_pixel = + min_num_peaks; + break; + } + } - for ( ipx=0; ipxn_patterns; cp++ ) { ImageFeatureList *flist = pattern_list->patterns[cp]->im_list; - if ( only_best_distance ) { + if ( gparams->only_best_distance ) { if ( fabs(pattern_list->patterns[cp]->clen-clen_to_use) > - 0.0001 ) { + 0.0001 ) { continue; } } @@ -970,87 +1059,99 @@ static int compute_pixel_statistics(struct pattern_list *pattern_list, double min_dist; double fx, fy; - double rfs, rss; - double crx, cry; Reflection *refl; + int ret; RefList *rlist = pattern_list->patterns[cp]->ref_list; struct imagefeature *imfe = image_get_feature(flist, ich); - compute_x_y(det, imfe->fs, imfe->ss, &fx, &fy); + compute_x_y(edet->det, imfe->fs, imfe->ss, &fx, &fy); - refl = find_closest_reflection(rlist, fx, fy, det, + refl = find_closest_reflection(rlist, fx, fy, edet->det, &min_dist); if ( refl == NULL ) continue; - if ( min_dist < max_peak_distance ) { + if ( min_dist < gparams->max_peak_dist ) { - int ipx = ((int)rint(imfe->fs) + aw* - (int)rint(imfe->ss)); - if ( num_pix_displ[ipx]>0 ) { + ret = add_distance_to_list(edet, imfe, + pix_displ_list, + refl, fx, fy); - int ret; + if ( ret != 0 ) return ret; - ret = allocate_next_element(curr_pix_displ, - ipx); + } + } + } - if ( ret != 0) return ret; + return 0; +} - } - get_detector_pos(refl, &rfs, &rss); - compute_x_y(det, rfs, rss, &crx, &cry); - get_detector_pos(refl, &rfs, &rss); - curr_pix_displ[ipx]->dfs = (fx-crx); - curr_pix_displ[ipx]->dss = (fy-cry); - curr_pix_displ[ipx]->ne = NULL; - num_pix_displ[ipx]++; - } - } +static int compute_avg_pix_displ(struct pixel_displ_list *pix_displ_list, + int num_peaks_per_pixel,int pix_index, + struct avg_displacements *avg_displ) +{ + int ret; + + if ( pix_displ_list->num_pix_displ[pix_index] >= + num_peaks_per_pixel ) { + + ret = fill_avg_pixel_displ(pix_displ_list, + pix_index, + avg_displ); + if ( ret !=0 ) return ret; + } else { + avg_displ->displ_x[pix_index] = -10000.0; + avg_displ->displ_y[pix_index] = -10000.0; + avg_displ->displ_abs[pix_index] = -10000.0; } - for ( np=min_num_peaks_per_pixel; np>0; np-- ) { - for ( di=0; din_rigid_groups; di++ ) { - if ( conn_data[di].num_peaks_per_pixel>np ) { - continue; - } + return 0; - for (ip=0; iprigid_groups[di]->n_panels; - ip++) { +} - int ret; - - ret = compute_panel_statistics(connected, - num_pix_displ, - curr_pix_displ, - all_pix_displ, - conn_data, di, ip, - np, - dfv, - aw, - displ_x, - displ_y, - displ_abs); - - if ( ret !=0 ) return ret; - } +static int compute_avg_displacements(struct enhanced_det *edet, + struct rg_collection *connected, + struct pixel_displ_list *pix_displ_list, + struct connected_data *conn_data, + struct avg_displacements *avg_displ) +{ + int di, ip, ifs, iss; + int pix_index, ret; + struct panel *p; + + for ( di=0; din_rigid_groups; di++ ) { + + for (ip=0; iprigid_groups[di]->n_panels; + ip++) { + + p = connected->rigid_groups[di]->panels[ip]; + + for ( ifs=p->min_fs; ifsmax_fs+1; ifs++ ) { + for ( iss=p->min_ss; issmax_ss+1; iss++ ) { + + pix_index = ifs+edet->width*iss; + + ret = compute_avg_pix_displ(pix_displ_list, + conn_data[di].num_peaks_per_pixel, + pix_index, avg_displ); - if ( conn_data[di].n_peaks_in_conn >= - min_num_peaks_per_panel ) { - conn_data[di].num_peaks_per_pixel = np; + if ( ret != 0 ) return ret; + + } } } - } - return 0; + } + return 0; } static double compute_error(struct rg_collection *connected, - int aw, + struct enhanced_det* edet, struct connected_data *conn_data, int *num_pix_displ, double *displ_abs) @@ -1072,13 +1173,17 @@ static double compute_error(struct rg_collection *connected, for (ifs=p->min_fs; ifsmax_fs+1; ifs++) { for (iss=p->min_ss; issmax_ss+1; iss++) { - if ( num_pix_displ[ifs+aw*iss]>= + + int pix_index; + pix_index = ifs+edet->width*iss; + + if ( num_pix_displ[pix_index]>= conn_data[di].num_peaks_per_pixel ) { double cer; - cer = displ_abs[ifs+aw*iss]* - displ_abs[ifs+aw*iss]; + cer = displ_abs[pix_index]* + displ_abs[pix_index]; connected_error += cer; num_connected_error++; total_error += cer; @@ -1093,11 +1198,10 @@ static double compute_error(struct rg_collection *connected, connected_error /= (double)num_connected_error; connected_error = sqrt(connected_error); - STATUS("Connected group %s: \n" - " Pixels with more than %d measurements: %i " - " Error: = %0.4f pixels.\n", - conn_data[di].name, conn_data[di].num_peaks_per_pixel, - num_connected_error, + STATUS("Error for connected group %s: %d pixels with " + "more than %d peaks: = %0.4f pixels.\n", + conn_data[di].name, num_connected_error, + conn_data[di].num_peaks_per_pixel, connected_error); } } @@ -1113,120 +1217,184 @@ static double compute_error(struct rg_collection *connected, } -static void fill_coordinate_matrices(struct detector *det, int aw, - double *slab_to_x, double *slab_to_y) +static struct pixel_maps *initialize_pixel_maps(struct enhanced_det *edet) { + int pi; + struct pixel_maps *pixel_maps; - for ( pi=0; pin_panels; pi++ ) { + pixel_maps = malloc(sizeof(struct pixel_maps)); + pixel_maps->pix_to_x = malloc(edet->num_pix*sizeof(double)); + if ( pixel_maps->pix_to_x == NULL ) { + free(pixel_maps); + return NULL; + } + + pixel_maps->pix_to_y = malloc(edet->num_pix*sizeof(double)); + if ( pixel_maps->pix_to_x == NULL ) { + ERROR("Failed to allocate memory for pixel maps.\n"); + free(pixel_maps->pix_to_x); + free(pixel_maps); + return NULL; + } + + for ( pi=0; pidet->n_panels; pi++ ) { struct panel *p; int iss, ifs; - p = &det->panels[pi]; + p = &edet->det->panels[pi]; - for (iss=p->min_ss; iss < p->max_ss+1; iss++) { - for (ifs=p->min_fs; ifs < p->max_fs+1; ifs++) { + for (iss=p->min_ss; iss < p->max_ss+1; iss++) { + for (ifs=p->min_fs; ifs < p->max_fs+1; ifs++) { double xs, ys; + int pix_index; - compute_x_y(det, ifs, iss, &xs, &ys); - - slab_to_x[iss*aw+ifs] = xs; - slab_to_y[iss*aw+ifs] = ys; + pix_index = iss*edet->width+ifs; + compute_x_y(edet->det, ifs, iss, &xs, &ys); + pixel_maps->pix_to_x[pix_index] = xs; + pixel_maps->pix_to_y[pix_index] = ys; } } } + + return pixel_maps; } -static int correct_empty_panels(struct rg_collection *quadrants, - struct rg_collection *connected, - int min_num_peaks_per_panel, - struct connected_data *conn_data) +void free_pixel_maps(struct pixel_maps* pixel_maps) { - double *aver_ang; - double *aver_str; - int *aver_num_ang; + free(pixel_maps->pix_to_x); + free(pixel_maps->pix_to_y); + free(pixel_maps); +} - int di,i; - aver_ang = malloc(quadrants->n_rigid_groups*sizeof(double)); - if ( aver_ang == NULL ) { + +struct avg_rot_and_stretch* initialize_avg_rot_stretch( + int num_rigid_groups) +{ + int i; + + struct avg_rot_and_stretch *avg_rot_str; + + avg_rot_str = malloc(sizeof(struct avg_rot_and_stretch)); + if ( avg_rot_str == NULL ) { ERROR("Failed to allocate memory to correct empty panels.|n"); - return 1; + return NULL; + } + + avg_rot_str->aver_ang = malloc(num_rigid_groups*sizeof(double)); + if ( avg_rot_str->aver_ang == NULL ) { + ERROR("Failed to allocate memory to correct empty panels.|n"); + free(avg_rot_str); + return NULL; } - aver_str = malloc(quadrants->n_rigid_groups*sizeof(double)); - if ( aver_str == NULL ) { + + avg_rot_str->aver_str = malloc(num_rigid_groups*sizeof(double)); + if ( avg_rot_str->aver_str == NULL ) { ERROR("Failed to allocate memory to correct empty panels.|n"); - free(aver_ang); - return 1; + free(avg_rot_str->aver_ang); + free(avg_rot_str); + return NULL; } - aver_num_ang = malloc(quadrants->n_rigid_groups*sizeof(int)); - if ( aver_num_ang == NULL ) { + + avg_rot_str->aver_num_ang = malloc(num_rigid_groups*sizeof(int)); + if ( avg_rot_str->aver_num_ang == NULL ) { ERROR("Failed to allocate memory to correct empty panels.|n"); - free(aver_ang); - free(aver_str); - return 1; + free(avg_rot_str->aver_ang); + free(avg_rot_str->aver_str); + free(avg_rot_str); + return NULL; + } + + for (i=0; iaver_ang[i] = 0.0; + avg_rot_str->aver_str[i] = 0.0; + avg_rot_str->aver_num_ang[i] = 0; } - for (i=0; in_rigid_groups; i++) { - aver_ang[i] = 0; - aver_str[i] = 0; - aver_num_ang[i] = 0; + return avg_rot_str; +} + + +void free_avg_angle_and_stretches(struct avg_rot_and_stretch* avg_rot_str) +{ + free(avg_rot_str->aver_ang); + free(avg_rot_str->aver_str); + free(avg_rot_str->aver_num_ang); +} + + +static int compute_rot_stretch_for_empty_panels( + struct rg_collection *quadrants, + struct rg_collection *connected, + struct geoptimiser_params *gparams, + struct connected_data *conn_data) +{ + struct avg_rot_and_stretch *avg_rot_str; + int di,i; + + STATUS("Computing rotation and elongation corrections for groups " + "without the required number of measurements.\n"); + + avg_rot_str = initialize_avg_rot_stretch(quadrants->n_rigid_groups); + if ( avg_rot_str == NULL ) { + return 1; } for (di=0; din_rigid_groups; di++) { - if ( conn_data[di].n_peaks_in_conn>=min_num_peaks_per_panel ) { - aver_ang[conn_data[di].num_quad] += conn_data[di].cang; - aver_str[conn_data[di].num_quad] += conn_data[di].cstr; - aver_num_ang[conn_data[di].num_quad]++; + if ( conn_data[di].n_peaks_in_conn >= + gparams->min_num_pix_per_conn_group ) { + avg_rot_str->aver_ang[conn_data[di].num_quad] += + conn_data[di].cang; + avg_rot_str->aver_str[conn_data[di].num_quad] += + conn_data[di].cstr; + avg_rot_str->aver_num_ang[conn_data[di].num_quad]++; } } for ( i=0; in_rigid_groups; i++ ) { - if ( aver_num_ang[i]>0 ) { - aver_ang[i] /= (double)aver_num_ang[i]; - aver_str[i] /= (double)aver_num_ang[i]; + if ( avg_rot_str->aver_num_ang[i] > 0 ) { + avg_rot_str->aver_ang[i] /= + (double)avg_rot_str->aver_num_ang[i]; + avg_rot_str->aver_str[i] /= + (double)avg_rot_str->aver_num_ang[i]; } } for ( di=0; din_rigid_groups; di++ ) { - if ( conn_data[di].n_peaks_in_conn0) { + if ( conn_data[di].n_peaks_in_conn < + gparams->min_num_pix_per_conn_group ) { + if ( avg_rot_str->aver_num_ang[conn_data[di].num_quad] + > 0) { conn_data[di].cang = - aver_ang[conn_data[di].num_quad]; + avg_rot_str->aver_ang[conn_data[di].num_quad]; conn_data[di].cstr = - aver_str[conn_data[di].num_quad]; - STATUS("Connected group %s has not enough peaks " - "(%i). Using average angle: %0.4f degrees\n", - conn_data[di].name, + avg_rot_str->aver_str[conn_data[di].num_quad]; + STATUS("Connected group %s has only %i useful " + "pixels. Using average angle: %0.4f " + "degrees\n", conn_data[di].name, conn_data[di].n_peaks_in_conn, conn_data[di].cang); } else { - STATUS("Connected group %s has not enough peaks " - "(%i). Left unchanged\n", + STATUS("Connected group %s does not have enough " + " peaks (%i). It will not be moved.\n", conn_data[di].name, conn_data[di].n_peaks_in_conn); } } } - free(aver_ang); - free(aver_str); - free(aver_num_ang); - - return 0; + return 0; } -static void correct_angle_and_stretch(struct rg_collection *connected, - struct detector *det, - struct connected_data *conn_data, - double use_clen, double stretch_coeff, - int individual_coffset) +static void correct_corner_coordinates(struct rg_collection *connected, + struct connected_data *conn_data) { int di, ip; @@ -1235,55 +1403,26 @@ static void correct_angle_and_stretch(struct rg_collection *connected, for ( ip=0; iprigid_groups[di]->n_panels; ip++ ) { struct panel *p; - double newx, newy; p = connected->rigid_groups[di]->panels[ip]; - newx = - p->fsx*cos(conn_data[di].cang)- - p->fsy*sin(conn_data[di].cang); - newy = - p->fsx*sin(conn_data[di].cang)+ - p->fsy*cos(conn_data[di].cang); - p->fsx = newx; - p->fsy = newy; - newx = - p->ssx*cos(conn_data[di].cang)- - p->ssy*sin(conn_data[di].cang); - newy = - p->ssx*sin(conn_data[di].cang)+ - p->ssy*cos(conn_data[di].cang); - p->ssx = newx; - p->ssy = newy; - } - } - - - if ( individual_coffset == 0 ) { - - int pi; + if ( ip == 0 ) { - for (pi=0; pin_panels; pi++) { - det->panels[pi].coffset -= use_clen*(1.0-stretch_coeff); - } - STATUS("Using a single offset distance for the whole detector: " - "%f m.\n", det->panels[0].coffset); + p->cnx *= conn_data[di].cstr; + p->cny *= conn_data[di].cstr; - for ( di=0; din_rigid_groups; di++ ) { - conn_data[di].cstr = stretch_coeff; - } + } else { - } else { + struct panel *p0; + double delta_x, delta_y; - STATUS("Using an individual distance for each connected group.\n"); - for ( di=0; din_rigid_groups; di++ ) { - for ( ip=0; iprigid_groups[di]->n_panels; - ip++ ) { + p0 = connected->rigid_groups[di]->panels[0]; - struct panel *p; + delta_x = p->cnx-p0->cnx/conn_data[di].cstr; + delta_y = p->cny-p0->cny/conn_data[di].cstr; - p = connected->rigid_groups[di]->panels[ip]; - p->coffset -= (1.0-conn_data[di].cstr)*p->clen; + p->cnx = p0->cnx + delta_x; + p->cny = p0->cny + delta_y; } } @@ -1291,90 +1430,125 @@ static void correct_angle_and_stretch(struct rg_collection *connected, } -static void shift_panels(struct rg_collection *connected, - struct connected_data *conn_data) +static void correct_rotation_and_stretch(struct rg_collection *connected, + struct detector *det, + struct connected_data *conn_data, + double clen_to_use, + double stretch_coeff, + struct geoptimiser_params *gparams) { int di, ip; + STATUS("Applying rotation and stretch corrections.\n"); for ( di=0; din_rigid_groups; di++ ) { for ( ip=0; iprigid_groups[di]->n_panels; ip++ ) { struct panel *p; + double new_fsx, new_fsy, new_ssx, new_ssy; p = connected->rigid_groups[di]->panels[ip]; - if ( ip == 0 ) { - - p->cnx *= conn_data[di].cstr; - p->cny *= conn_data[di].cstr; + new_fsx = p->fsx*cos(conn_data[di].cang)- + p->fsy*sin(conn_data[di].cang); + new_fsy = p->fsx*sin(conn_data[di].cang)+ + p->fsy*cos(conn_data[di].cang); + new_ssx = p->ssx*cos(conn_data[di].cang)- + p->ssy*sin(conn_data[di].cang); + new_ssy = p->ssx*sin(conn_data[di].cang)+ + p->ssy*cos(conn_data[di].cang); + p->fsx = new_fsx; + p->fsy = new_fsy; + p->ssx = new_ssx; + p->ssy = new_ssy; + } + } - } else { - struct panel *p0; - double delta_x, delta_y; + if ( gparams->individual_coffset ) { - p0 = connected->rigid_groups[di]->panels[0]; + STATUS("Using individual distances for rigid panels.\n"); + for ( di=0; din_rigid_groups; di++ ) { + for ( ip=0; iprigid_groups[di]->n_panels; + ip++ ) { - delta_x = (p->cnx-p0->cnx/conn_data[di].cstr); - delta_y = (p->cny-p0->cny/conn_data[di].cstr); + struct panel *p; - p->cnx = p0->cnx + delta_x * cos(conn_data[di].cang) - - delta_y * sin(conn_data[di].cang); - p->cny = p0->cny + delta_x * sin(conn_data[di].cang) - + delta_y * cos(conn_data[di].cang); + p = connected->rigid_groups[di]->panels[ip]; + p->coffset -= (1.0-conn_data[di].cstr)*clen_to_use; } } - } -} + } else { -static void recompute_panel(struct connected_data *conn_data, int di, int ip, - struct rg_collection *connected, - double *slab_to_x, - double *slab_to_y, - double *recomputed_slab_to_x, - double *recomputed_slab_to_y, - double *displ_x, double *displ_y, - double stretch_coeff, int aw, int *num_pix_displ) -{ + int pi; + for ( di=0; din_rigid_groups; di++ ) { + conn_data[di].cstr = stretch_coeff; + } + + for ( pi=0; pin_panels; pi++) { + det->panels[pi].coffset -= clen_to_use*(1.0-stretch_coeff); + } + STATUS("Using a single offset distance for the whole " + "detector: %f m. Stretch ceofficient: %0.4f\n", + det->panels[0].coffset, stretch_coeff); + } + + correct_corner_coordinates(connected, conn_data); +} + + +static void adjust_panel(struct connected_data *conn_data, int di, int ip, + struct rg_collection *connected, + struct pixel_maps *pixel_maps, + struct pixel_maps *recomputed_pixel_maps, + struct avg_displacements *avg_displ, + double stretch_coeff, + struct enhanced_det *edet, int *num_pix_displ) +{ double c_stretch; struct panel *p; int ifs, iss; c_stretch = conn_data[di].cstr; - if ( fabs(c_stretch)rigid_groups[di]->panels[ip]; for ( ifs=p->min_fs; ifsmax_fs+1; ifs++ ) { for ( iss=p->min_ss; issmax_ss+1; iss++ ) { - recomputed_slab_to_x[ifs+aw*iss] /= c_stretch; - recomputed_slab_to_y[ifs+aw*iss] /= c_stretch; - if ( num_pix_displ[ifs+aw*iss] >= + + int pix_index = ifs+edet->width*iss; + + recomputed_pixel_maps->pix_to_x[pix_index] /= c_stretch; + recomputed_pixel_maps->pix_to_y[pix_index] /= c_stretch; + if ( num_pix_displ[pix_index] >= conn_data[di].num_peaks_per_pixel) { - displ_x[ifs+aw*iss] -= (slab_to_x[ifs+aw*iss]- - recomputed_slab_to_x[ifs+aw*iss]); - displ_y[ifs+aw*iss] -= (slab_to_y[ifs+aw*iss]- - recomputed_slab_to_y[ifs+aw*iss]); + avg_displ->displ_x[pix_index] -= + (pixel_maps->pix_to_x[pix_index]- + recomputed_pixel_maps->pix_to_x[pix_index]); + avg_displ->displ_y[pix_index] -= + (pixel_maps->pix_to_y[pix_index]- + recomputed_pixel_maps->pix_to_y[pix_index]); } } } } -static void recompute_differences(struct rg_collection *connected, - double *slab_to_x, double *slab_to_y, - double *recomputed_slab_to_x, - double *recomputed_slab_to_y, - struct connected_data *conn_data, - double stretch_coeff, int aw, - double *displ_x, double *displ_y, - int *num_pix_displ) +static void adjust_displ_for_stretch(struct rg_collection *connected, + struct pixel_maps *pixel_maps, + struct pixel_maps *recomputed_pixel_maps, + struct connected_data *conn_data, + double stretch_coeff, + struct enhanced_det *edet, + struct avg_displacements *avg_displ, + int *num_pix_displ) { int di, ip; @@ -1382,47 +1556,57 @@ static void recompute_differences(struct rg_collection *connected, for ( di=0; din_rigid_groups; di++ ) { for (ip=0; iprigid_groups[di]->n_panels; ip++) { - - recompute_panel(conn_data, di, ip, connected, - slab_to_x, slab_to_y, - recomputed_slab_to_x, - recomputed_slab_to_y, - displ_x, displ_y, - stretch_coeff, aw, num_pix_displ); + adjust_panel(conn_data, di, ip, connected, + pixel_maps, recomputed_pixel_maps, + avg_displ, stretch_coeff, edet, + num_pix_displ); } } } -static void fill_av_in_panel(struct rg_collection *connected, int di, int ip, - struct connected_data *conn_data, - int *num_pix_displ, int aw, - double *av_in_panel_fs, - double *av_in_panel_ss, - double *displ_x, double *displ_y) +static void fill_av_conn(struct rg_collection *connected, int di, + struct connected_data *conn_data, + int *num_pix_displ, struct enhanced_det *edet, + double *list_displ_in_conn_fs, + double *list_displ_in_conn_ss, + struct avg_displacements *avg_displ) { struct panel *p; - int ifs, iss; + int ifs, iss, ip; - p = connected->rigid_groups[di]->panels[ip]; + int counter = 0; - for ( ifs=p->min_fs; ifsmax_fs+1; ifs++ ) { - for ( iss=p->min_ss; issmax_ss+1; iss++ ) { + for (ip=0; iprigid_groups[di]->n_panels; ip++) { - if (num_pix_displ[ifs+aw*iss]>= - conn_data[di].num_peaks_per_pixel) { - av_in_panel_fs[conn_data[di].n_peaks_in_conn] = - displ_x[ifs+aw*iss]; - av_in_panel_ss[conn_data[di].n_peaks_in_conn] = - displ_y[ifs+aw*iss]; - conn_data[di].n_peaks_in_conn++; + p = connected->rigid_groups[di]->panels[ip]; + + for ( ifs=p->min_fs; ifsmax_fs+1; ifs++ ) { + for ( iss=p->min_ss; issmax_ss+1; iss++ ) { + + int pix_index = ifs+edet->width*iss; + + if ( num_pix_displ[pix_index] >= + conn_data[di].num_peaks_per_pixel ) { + list_displ_in_conn_fs[counter] = + avg_displ->displ_x[pix_index]; + list_displ_in_conn_ss[counter] = + avg_displ->displ_y[pix_index]; + counter++; + } } } } + + if ( counter != conn_data[di].n_peaks_in_conn ) { + printf("counter: %i n_peaks_in_conn: %i\n", + counter, conn_data[di].n_peaks_in_conn); + exit(0); + } } -static void fill_con_data_sh(struct connected_data *conn_data, +static void fill_conn_data_sh(struct connected_data *conn_data, double *av_in_panel_fs, double *av_in_panel_ss, int di, double max_peak_distance) @@ -1431,11 +1615,11 @@ static void fill_con_data_sh(struct connected_data *conn_data, conn_data[di].n_peaks_in_conn); conn_data[di].sh_y = comp_median(av_in_panel_ss, conn_data[di].n_peaks_in_conn); - STATUS("Panel %s, num of considered pixels: %i, shifts (in pixels) x,y: %0.8f, %0.8f\n", + STATUS("Panel %s, num pixels: %i, shifts (in pixels) X,Y: %0.8f, %0.8f\n", conn_data[di].name, conn_data[di].n_peaks_in_conn, conn_data[di].sh_x, conn_data[di].sh_y); - if ( modulus2d(conn_data[di].sh_x, conn_data[di].sh_y) > - 0.8*max_peak_distance ) { + if ( modulus2d(conn_data[di].sh_x, conn_data[di].sh_y ) > + 0.8*max_peak_distance ) { STATUS(" WARNING: absolute shift is: %0.1f > 0.8*%0.1f pixels." " Increase the value of the max_peak_distance parameter!\n", modulus2d(conn_data[di].sh_x, conn_data[di].sh_y), @@ -1444,162 +1628,183 @@ static void fill_con_data_sh(struct connected_data *conn_data, } -static int compute_shifts(struct rg_collection *connected, - struct connected_data *conn_data, - int *num_pix_displ, int aw, - int min_num_peaks_per_panel, - double dfv, double max_peak_distance, - double *displ_x, double *displ_y) +static int compute_shift(struct rg_collection *connected, + struct connected_data *conn_data, + int *num_pix_displ, + struct enhanced_det *edet, + struct geoptimiser_params *gparams, + struct avg_displacements *avg_displ) { - STATUS("Computing average shifts for connected groups.\n"); + STATUS("Computing shift corrections.\n"); - int di, ip; + int di; for ( di=0; din_rigid_groups; di++ ) { - int cmaxfs; - int cmaxss; - int num_all_pixels; - double *av_in_panel_fs; - double *av_in_panel_ss; - - cmaxfs = connected->rigid_groups[di]->panels[0]->max_fs+1- - connected->rigid_groups[di]->panels[0]->min_fs; - cmaxss = connected->rigid_groups[di]->panels[0]->max_ss+1- - connected->rigid_groups[di]->panels[0]->min_ss; + double *list_displ_in_conn_fs; + double *list_displ_in_conn_ss; - num_all_pixels = cmaxfs*cmaxss*connected->rigid_groups[di]->n_panels; - - av_in_panel_fs = malloc(num_all_pixels*sizeof(double)); - if ( av_in_panel_fs == NULL ) { + list_displ_in_conn_fs = malloc(conn_data[di].n_peaks_in_conn* + sizeof(double)); + if ( list_displ_in_conn_fs == NULL ) { ERROR("Failed to allocate memory for computing shifts\n"); return 1; } - av_in_panel_ss = malloc(num_all_pixels*sizeof(double)); - if ( av_in_panel_ss == NULL ) { + list_displ_in_conn_ss = malloc(conn_data[di].n_peaks_in_conn* + sizeof(double)); + if ( list_displ_in_conn_ss == NULL ) { ERROR("Failed to allocate memory for computing shifts\n"); - free(av_in_panel_fs); + free(list_displ_in_conn_fs); return 1; } - conn_data[di].n_peaks_in_conn = 0; - - for (ip=0; iprigid_groups[di]->n_panels; ip++) { - - fill_av_in_panel(connected, di, ip, conn_data, - num_pix_displ, aw, av_in_panel_fs, - av_in_panel_ss, displ_x, displ_y); - - } - - if ( conn_data[di].n_peaks_in_conn>=min_num_peaks_per_panel ) { + fill_av_conn(connected, di, conn_data, num_pix_displ, edet, + list_displ_in_conn_fs, list_displ_in_conn_ss, + avg_displ); - fill_con_data_sh(conn_data, av_in_panel_fs, - av_in_panel_ss, di, - max_peak_distance); + if ( conn_data[di].n_peaks_in_conn >= + gparams->min_num_pix_per_conn_group ) { + fill_conn_data_sh(conn_data, list_displ_in_conn_fs, + list_displ_in_conn_ss, di, + gparams->max_peak_dist); } else { - conn_data[di].sh_x = dfv; - conn_data[di].sh_y = dfv; + conn_data[di].sh_x = -10000.0; + conn_data[di].sh_y = -10000.0; } - free(av_in_panel_fs); - free(av_in_panel_ss); - } + free(list_displ_in_conn_fs); + free(list_displ_in_conn_ss); + } return 0; } -static int compute_shifts_for_empty_panels(struct rg_collection *quadrants, - struct rg_collection *connected, - struct connected_data *conn_data, - int min_num_peaks_per_panel) +struct avg_shift *initialize_avg_shift(int num_rigid_groups) { + struct avg_shift *avg_sh; + int i; - double *aver_x; - double *aver_y; - int *num_aver; - int di, i; + avg_sh = malloc(sizeof(struct avg_shift)); + if ( avg_sh == NULL ) { + ERROR("Failed to allocate memory to compute shifts for " + "empty panels.\n"); + return NULL; + } - // shifts for empty - aver_x = malloc(quadrants->n_rigid_groups*sizeof(double)); - if ( aver_x == NULL ) { + avg_sh->aver_x = malloc(num_rigid_groups*sizeof(double)); + if ( avg_sh->aver_x == NULL ) { ERROR("Failed to allocate memory to compute shifts for " "empty panels.\n"); - return 1; + free(avg_sh); + return NULL; } - aver_y = malloc(quadrants->n_rigid_groups*sizeof(double)); - if ( aver_y == NULL ) { + avg_sh->aver_y = malloc(num_rigid_groups*sizeof(double)); + if ( avg_sh->aver_y == NULL ) { ERROR("Failed to allocate memory to compute shifts for " "empty panels.\n"); - free(aver_x); - return 1; + free(avg_sh->aver_x); + free(avg_sh); + return NULL; } - num_aver = malloc(quadrants->n_rigid_groups*sizeof(int)); - if ( num_aver == NULL ) { + avg_sh->aver_num_sh = malloc(num_rigid_groups*sizeof(int)); + if ( avg_sh->aver_num_sh == NULL ) { ERROR("Failed to allocate memory to compute shifts for " "empty panels.\n"); - free(aver_x); - free(aver_y); - return 1; + free(avg_sh->aver_x); + free(avg_sh->aver_y); + free(avg_sh); + return NULL; } - for ( i=0; in_rigid_groups; i++ ) { - aver_x[i] = 0; - aver_y[i] = 0; - num_aver[i] = 0; + for ( i=0; iaver_x[i] = 0.0; + avg_sh->aver_y[i] = 0.0; + avg_sh->aver_num_sh[i] = 0; } + return avg_sh; +} + + +void free_avg_shift(struct avg_shift *av_sh) { + free(av_sh->aver_x); + free(av_sh->aver_y); + free(av_sh->aver_num_sh); + free(av_sh); +} + + +static int compute_shift_for_empty_panels(struct rg_collection *quadrants, + struct rg_collection *connected, + struct connected_data *conn_data, + struct geoptimiser_params* gparams) +{ + + struct avg_shift *av_sh; + int di, i; + + STATUS("Computing shift corrections for groups without the required " + "number of measurements.\n"); + av_sh = initialize_avg_shift(quadrants->n_rigid_groups); + if ( av_sh == NULL ) return 1; + for ( di=0; din_rigid_groups; di++ ) { - if ( conn_data[di].n_peaks_in_conn>=min_num_peaks_per_panel ) { - aver_x[conn_data[di].num_quad] += conn_data[di].sh_x; - aver_y[conn_data[di].num_quad] += conn_data[di].sh_y; - num_aver[conn_data[di].num_quad]++; + if ( conn_data[di].n_peaks_in_conn >= + gparams->min_num_pix_per_conn_group ) { + av_sh->aver_x[conn_data[di].num_quad] += + conn_data[di].sh_x; + av_sh->aver_y[conn_data[di].num_quad] += + conn_data[di].sh_y; + av_sh->aver_num_sh[conn_data[di].num_quad]++; } } for ( i=0; in_rigid_groups; i++ ) { - if (num_aver[i]>0) { - aver_x[i] /= (double)num_aver[i]; - aver_y[i] /= (double)num_aver[i]; + if (av_sh->aver_num_sh[i]>0) { + av_sh->aver_x[i] /= (double)av_sh->aver_num_sh[i]; + av_sh->aver_y[i] /= (double)av_sh->aver_num_sh[i]; } } for (di=0; din_rigid_groups; di++) { - if ( conn_data[di].n_peaks_in_conn0 ) { - conn_data[di].sh_x = aver_x[conn_data[di].num_quad]; - conn_data[di].sh_y = aver_y[conn_data[di].num_quad]; - STATUS("Panel %s has not enough (%i) peaks. " - "Using average shifts (in pixels) X,Y: " - "%0.2f,%0.2f\n", conn_data[di].name, + if ( conn_data[di].n_peaks_in_conn < + gparams->min_num_pix_per_conn_group ) { + if ( av_sh->aver_num_sh[conn_data[di].num_quad] > 0 ) { + conn_data[di].sh_x = + av_sh->aver_x[conn_data[di].num_quad]; + conn_data[di].sh_y = + av_sh->aver_y[conn_data[di].num_quad]; + STATUS("Panel %s doesn't not have eough (%i) " + "peaks. Using average shifts (in pixels) " + "X,Y: %0.2f,%0.2f\n", conn_data[di].name, conn_data[di].n_peaks_in_conn, conn_data[di].sh_x, conn_data[di].sh_y); } else { STATUS("Panel %s has not enough (%i) peaks. " - "Left unchanged.\n", + "It will not be moved.\n", conn_data[di].name, conn_data[di].n_peaks_in_conn); } } } - free(aver_x); - free(aver_y); - free(num_aver); + free_avg_shift(av_sh); return 0; } -static void correct_shifts(struct rg_collection *connected, - struct connected_data *conn_data, - double dfv, double clen_to_use) +static void correct_shift(struct rg_collection *connected, + struct connected_data *conn_data, + double clen_to_use) { int di; int ip; + STATUS("Applying shift corrections.\n"); + for ( di=0;din_rigid_groups;di++ ) { for (ip=0; iprigid_groups[di]->n_panels; ip++) { @@ -1607,300 +1812,246 @@ static void correct_shifts(struct rg_collection *connected, p = connected->rigid_groups[di]->panels[ip]; - if ( conn_data[di].sh_x>dfv+1.0 && - conn_data[di].sh_y > dfv+1.0 ) { + if ( conn_data[di].sh_x > -9999.0 && + conn_data[di].sh_y > -9999.0 ) { p->cnx -= conn_data[di].sh_x; p->cny -= conn_data[di].sh_y; } else { - STATUS("For some reason connected group %s is " - "empty! Left unchanged.\n", - p->name); + STATUS("For some reason shifts for panel %s has " + "not been computed!\n", p->name); } } } } -static void a_s_counting_loop(int *num_pix_displ, int ifs, int iss, - int di, struct connected_data *conn_data, - int aw, double *slab_to_x, - double *slab_to_y, struct panel *p0, - struct panel *p1, double *displ_x, - double *displ_y, double minrad, - int *num_ang) +static struct connected_stretch_and_angles *initialize_connected_stretch_angles( + struct rg_collection *connected) { - double coX, coY, cdX, cdY; - - if ( num_pix_displ[ifs+aw*iss]>= - conn_data[di].num_peaks_per_pixel ) { - - int ifs1, iss1; - int max_fs1_tmp = p0->max_fs; - int max_ss1_tmp = p0->max_ss; - - coX = slab_to_x[ifs+aw*iss]; - coY = slab_to_y[ifs+aw*iss]; - cdX = coX - displ_x[ifs+aw*iss]; - cdY = coY - displ_y[ifs+aw*iss]; - - for (ifs1=ifs+1; ifs1max_fs; - } - - for (iss1=iss+1; iss1max_ss; - } - - if ( num_pix_displ[ifs1+aw*iss1]>= - conn_data[di].num_peaks_per_pixel ) { - - double dist; - double coX1, coY1, cdX1, cdY1; - double len1, len2; - - dist = modulus2d(ifs-ifs1,iss-iss1); - if ( dist < minrad ) continue; - coX1 = slab_to_x[ifs1+aw*iss1]; - coY1 = slab_to_y[ifs1+aw*iss1]; - cdX1 = - coX1 - displ_x[ifs1+aw*iss1]; - cdY1 = - coY1 - displ_y[ifs1+aw*iss1]; - - len1 = modulus2d(coX1-coX, coY1-coY); - len2 = modulus2d(cdX1-cdX, cdY1-cdY); - if ( len1stretch_coeff = malloc(connected->n_rigid_groups*sizeof(double)); + if ( csaa->stretch_coeff == NULL ) { + ERROR("Failed to allocate memory to compute angles and " + "stretch.\n"); + free(csaa); + return NULL; + } + csaa->num_angles = malloc(connected->n_rigid_groups*sizeof(long)); + if ( csaa->num_angles == NULL ) { + ERROR("Failed to allocate memory to compute angles and " + "stretch.\n"); + free(csaa->stretch_coeff); + free(csaa); + return NULL; } + csaa->num_coeff=0; + return csaa; } -static int a_s_processing_loop(int *num_pix_displ, int ifs, int iss, - int di, struct connected_data *conn_data, - int aw, double *slab_to_x, - double *slab_to_y, struct panel *p0, - struct panel *p1, double *displ_x, - double *displ_y, double minrad, - int max_num_ang, int *num_ang, - double *angles, double *stretches) +static void scan_p1(int ip0, int ip1, struct rg_collection *connected, + struct avg_displacements *avg_displ, + struct connected_data *conn_data, + struct enhanced_det *edet, + struct pixel_maps *pixel_maps, + int* num_pix_displ, + int di, double min_dist, + long *num_ang, int ifs0, int iss0, + double c_x0, double c_y0, double cd_x0, double cd_y0, + int compute, double *angles, double *stretches) { - double coX, coY, cdX, cdY; - if ( num_pix_displ[ifs+aw*iss]>= - conn_data[di].num_peaks_per_pixel ) { + int iss1, ifs1; - int ifs1, iss1; - int max_fs1_tmp = p0->max_fs; - int max_ss1_tmp = p0->max_ss; + struct panel *p1 = connected->rigid_groups[di]->panels[ip1]; - if ( *num_ang>=max_num_ang ) return -2; + int min_ss_p1, min_fs_p1; - coX = slab_to_x[ifs+aw*iss]; - coY = slab_to_y[ifs+aw*iss]; - cdX = coX - displ_x[ifs+aw*iss]; - cdY = coY - displ_y[ifs+aw*iss]; + if ( ip0 == ip1 ) { + min_fs_p1 = ifs0; + min_ss_p1 = iss0; + } else { + min_fs_p1 = p1->min_fs; + min_ss_p1 = p1->min_ss; + } - for (ifs1=ifs+1; ifs1max_ss+1; iss1++ ) { - if ( ifs1 == max_fs1_tmp ) { - max_fs1_tmp = p1->max_fs; - } + for ( ifs1=min_fs_p1; ifs1max_fs+1; ifs1++ ) { - for (iss1=iss+1; iss1width*iss1; - if ( iss1 == max_ss1_tmp ) { - max_ss1_tmp = p1->max_ss; + if ( num_pix_displ[pix_index1]>= + conn_data[di].num_peaks_per_pixel ) { + + double dist; + double c_x1, c_y1, cd_x1, cd_y1; + double d_c_x, d_c_y, d_cd_x, d_cd_y; + double len1, len2; + + c_x1 = pixel_maps->pix_to_x[pix_index1]; + c_y1 = pixel_maps->pix_to_y[pix_index1]; + cd_x1 = c_x1 - avg_displ->displ_x[pix_index1]; + cd_y1 = c_y1 - avg_displ->displ_y[pix_index1]; + d_c_x = c_x1-c_x0; + d_c_y = c_y1-c_y0; + d_cd_x = cd_x1-cd_x0; + d_cd_y = cd_y1-cd_y0; + + dist = modulus2d(d_c_x,d_c_y); + if ( dist < min_dist ) continue; + + len1 = modulus2d(d_c_x, d_c_y); + len2 = modulus2d(d_cd_x, d_cd_y); + if ( len1= - conn_data[di].num_peaks_per_pixel ) { + if (compute) { - double dist; - double coX1, coY1, cdX1, cdY1; - double len1, len2; - double scalM; + double scal_m; double multlen; - if ( *num_ang>=max_num_ang ) return -2; - dist = modulus2d(ifs-ifs1,iss-iss1); - if (dist=multlen ) { + if ( fabs(scal_m)>=multlen ) { angles[*num_ang] = 0.0; } else { - angles[*num_ang] = 1.0; + angles[*num_ang] = acos(scal_m + /multlen); - angles[*num_ang] = - acos(scalM/multlen); - - if ((coY1-coY)*(cdX1-cdX)- - (coX1-coX)*(cdY1-cdY) < 0) { + if (d_c_y * d_cd_x - + d_c_x * d_cd_y < 0) { angles[*num_ang] *= -1.; } } stretches[*num_ang] = len1/len2; - - *num_ang = *num_ang+1; } + + *num_ang = *num_ang+1; } } } - return 0; } +static void scan_p0(int ip0, struct rg_collection *connected, + struct avg_displacements *avg_displ, + struct connected_data *conn_data, + struct enhanced_det *edet, + struct pixel_maps *pixel_maps, + int* num_pix_displ, + int di, double min_dist, + long *num_ang, int compute, + double *angles, double *stretches) +{ + + int iss0, ifs0, ip1; + + struct panel *p0 = connected->rigid_groups[di]->panels[ip0]; + + for ( iss0=p0->min_ss; iss0max_ss+1; iss0++ ) { + + for ( ifs0=p0->min_fs; ifs0max_fs+1; ifs0++ ) { + + int pix_index0 = ifs0+edet->width*iss0; + + if ( num_pix_displ[pix_index0]>= + conn_data[di].num_peaks_per_pixel ) { + + double c_x0, c_y0, cd_x0, cd_y0; + + c_x0 = pixel_maps->pix_to_x[pix_index0]; + c_y0 = pixel_maps->pix_to_y[pix_index0]; + cd_x0 = c_x0 - avg_displ->displ_x[pix_index0]; + cd_y0 = c_y0 - avg_displ->displ_y[pix_index0]; + + for ( ip1=ip0; + ip1rigid_groups[di]->n_panels; + ip1++ ) { + scan_p1(ip0, ip1, connected, avg_displ, + conn_data, edet, pixel_maps, + num_pix_displ, di, min_dist, + num_ang, ifs0, iss0, c_x0, + c_y0, cd_x0, cd_y0, compute, + angles, stretches); + } + } + } + } +} -static int compute_angles_and_stretch(struct rg_collection *connected, - struct connected_data *conn_data, - int *num_pix_displ, - double *slab_to_x, - double *slab_to_y, - double *displ_x, - double *displ_y, - int aw, - int min_num_peaks_per_panel, - double dist_coeff_ang_str, - int num_peaks_per_pixel, - double man_stretching_coeff, - double *stretch_coeff) +static double compute_rotation_and_stretch(struct rg_collection *connected, + struct connected_data *conn_data, + struct enhanced_det *edet, + int *num_pix_displ, + struct pixel_maps *pixel_maps, + struct avg_displacements *avg_displ, + double dist_coeff_for_rot_str, + struct geoptimiser_params *gparams) { int di; - int num_coeff; double stretch_cf; struct connected_stretch_and_angles *csaa; - csaa = malloc(sizeof(struct connected_stretch_and_angles)); + STATUS("Computing rotation and stretch corrections.\n"); + + csaa = initialize_connected_stretch_angles(connected); if ( csaa == NULL ) { - ERROR("Failed to allocate memory to compute angles and " - "stretch.\n"); - return 1; - } - csaa->stretch_coeff = malloc(connected->n_rigid_groups*sizeof(double)); - if ( csaa->stretch_coeff == NULL ) { - ERROR("Failed to allocate memory to compute angles and " - "stretch.\n"); - free(csaa); - return 1; - } - csaa->num_angles = malloc(connected->n_rigid_groups*sizeof(unsigned int)); - if ( csaa->num_angles == NULL ) { - ERROR("Failed to allocate memory to compute angles and " - "stretch.\n"); - free(csaa->stretch_coeff); - free(csaa); - return 1; + return -1.0; } - csaa->num_coeff=0; - for ( di=0; din_rigid_groups; di++ ) { - if ( conn_data[di].n_peaks_in_connmin_num_pix_per_conn_group ) { continue; } - unsigned int max_num_ang = 0; + long max_num_ang = 0; + double min_dist; double* angles; double* stretches; - int cmaxfs; - int cmaxss; - - struct panel *p; - double minrad; - int num_ang = 0; - int ip0, ip1; + struct panel *first_p; + long num_ang = 0; + int ip0; + int num_pix_first_p; - p = connected->rigid_groups[di]->panels[0]; + first_p = connected->rigid_groups[di]->panels[0]; - cmaxfs = p->max_fs+1-p->min_fs; - cmaxss = p->max_ss+1-p->min_ss; + num_pix_first_p = (first_p->max_fs+1-first_p->min_fs)* + (first_p->max_ss+1-first_p->min_ss); // TODO: MINRAD HERE IS NOT UNIVERSAL - minrad = dist_coeff_ang_str*sqrt(cmaxfs*cmaxss* - connected->rigid_groups[di]->n_panels); + min_dist = dist_coeff_for_rot_str*sqrt(num_pix_first_p* + connected->rigid_groups[di]->n_panels); for ( ip0=0; ip0rigid_groups[di]->n_panels; ip0++ ) { - struct panel *p0 = connected->rigid_groups[di]->panels[ip0]; - - for ( ip1=0; ip1rigid_groups[di]->n_panels; - ip1++ ) { - - struct panel *p1 = - connected->rigid_groups [di]->panels[ip1]; - - int ifs, iss; - int min_fs_tmp = p0->min_fs; - int max_fs_tmp = p0->max_fs; - int min_ss_tmp = p0->min_ss; - int max_ss_tmp = p0->max_ss; - - for (ifs=min_fs_tmp; ifsmin_fs; - max_fs_tmp = p1->max_fs; - } - - for (iss=min_ss_tmp; issmin_ss; - max_ss_tmp = p1->max_ss; - } - - a_s_counting_loop(num_pix_displ, - ifs, iss, di, conn_data, - aw, slab_to_x, slab_to_y, - p0, p1, displ_x, - displ_y, minrad, - &num_ang); - - } - } - } + scan_p0(ip0, connected, avg_displ, conn_data, edet, + pixel_maps, num_pix_displ, di, min_dist, + &num_ang, 0, NULL, NULL); } - max_num_ang = conn_data[di].n_peaks_in_conn* - conn_data[di].n_peaks_in_conn; max_num_ang = num_ang+1; angles = malloc(max_num_ang*sizeof(double)); @@ -1910,7 +2061,7 @@ static int compute_angles_and_stretch(struct rg_collection *connected, free(csaa->stretch_coeff); free(csaa->num_angles); free(csaa); - return 1; + return -1.0; } stretches = malloc(max_num_ang*sizeof(double)); if ( stretches == NULL ) { @@ -1920,71 +2071,25 @@ static int compute_angles_and_stretch(struct rg_collection *connected, free(csaa->stretch_coeff); free(csaa->num_angles); free(csaa); - return 1; + return -1.0; } num_ang = 0; for ( ip0=0; ip0rigid_groups[di]->n_panels; ip0++ ) { - struct panel *p0 = connected->rigid_groups[di]->panels[ip0]; - - for ( ip1=0; ip1rigid_groups[di]->n_panels; - ip1++ ) { - - struct panel *p1 = - connected->rigid_groups [di]->panels[ip1]; - - int ifs, iss; - int min_fs_tmp = p0->min_fs; - int max_fs_tmp = p0->max_fs; - int min_ss_tmp = p0->min_ss; - int max_ss_tmp = p0->max_ss; - - for (ifs=min_fs_tmp; ifsmin_fs; - max_fs_tmp = p1->max_fs; - } - - for (iss=min_ss_tmp; issmin_ss; - max_ss_tmp = p1->max_ss; - } - - ret = a_s_processing_loop( - num_pix_displ, - ifs, iss, di, - conn_data, - aw, slab_to_x, - slab_to_y, - p0, p1, displ_x, - displ_y, minrad, - max_num_ang, - &num_ang, angles, - stretches); - - if ( ret == -2 ) break; - } - } - } + scan_p0(ip0, connected, avg_displ, conn_data, edet, + pixel_maps, num_pix_displ, di, min_dist, + &num_ang, 1, angles, stretches); } - if ( num_ang<1 ) continue; - conn_data[di].cang = -comp_median(angles,num_ang); - conn_data[di].cstr = comp_median(stretches,num_ang); + if ( num_ang < 1 ) continue; + conn_data[di].cang = -comp_median(angles, num_ang); + conn_data[di].cstr = comp_median(stretches, num_ang); - STATUS("Connected group %s:\n" - " Number of measurements: %i Angle: %0.4f deg " - "Stretch coefficient: %0.4f\n", - conn_data[di].name, num_ang, conn_data[di].cang, - conn_data[di].cstr); + STATUS("Panel %s, num: %li, angle: %0.4f deg, stretch coeff: " + "%0.4f\n", conn_data[di].name, num_ang, conn_data[di].cang, + conn_data[di].cstr); csaa->stretch_coeff[csaa->num_coeff] = conn_data[di].cstr; csaa->num_angles[csaa->num_coeff] = num_ang; @@ -1994,55 +2099,57 @@ static int compute_angles_and_stretch(struct rg_collection *connected, free(stretches); } + stretch_cf = 1.0; - num_coeff = csaa->num_coeff; + printf("Computing overall stretch coefficient.\n"); - stretch_cf = 1; - if (num_coeff>0) { + if ( csaa->num_coeff>0 ) { - int ipp; + int peaks_per_p; - for ( ipp=num_peaks_per_pixel; ipp>=0; ipp-- ) { + peaks_per_p = gparams->min_num_peaks_per_pix; + + while ( peaks_per_p>=0 ) { double total_num; - int di; + long di; + stretch_cf = 0; total_num = 0; - for ( di=0; di=ipp ) { + for ( di=0; dinum_coeff; di++ ) { + if ( conn_data[di].num_peaks_per_pixel >= + peaks_per_p ) { + stretch_cf += csaa->stretch_coeff[di]* + (double)csaa->num_angles[di]; total_num += csaa->num_angles[di]; } } - if ( total_num>1 ) { - total_num = 1./total_num; - } else { - continue; - } - stretch_cf = 0; - for ( di=0; di=ipp ) { - stretch_cf += - total_num*csaa->stretch_coeff[di]* - (double)csaa->num_angles[di]; - } + + if ( total_num > 0 ) { + + stretch_cf /= total_num; + + printf("(Using only connected groups for which " + "the minimum number of measurements per " + "pixel is %i).\n", peaks_per_p); + break; } - break; + peaks_per_p--; } } - if ( stretch_cfFLT_EPSILON ) { - stretch_cf = man_stretching_coeff; - STATUS("Using manually set stretch coefficient: %0.4f\n", - stretch_cf); - + if ( gparams->nostretch ) { + stretch_cf = 1.0; + STATUS("However, distance offset will not be optimized, so the " + "stretching coefficient has been set to 1.0\n"); for ( di=0; din_rigid_groups; di++ ) { - conn_data[di].cstr = man_stretching_coeff; + conn_data[di].cstr = stretch_cf; } } @@ -2051,9 +2158,7 @@ static int compute_angles_and_stretch(struct rg_collection *connected, free(csaa->num_angles); free(csaa); - *stretch_coeff = stretch_cf; - - return 0; + return stretch_cf; } @@ -2114,19 +2219,18 @@ static int unpack_slab(struct image *image) } for ( ss=0; ssh; ss++ ) { - for ( fs=0; fsw; fs++ ) { - - int idx; - int cfs, css; + for ( fs=0; fsw; fs++ ) { - cfs = fs+p->min_fs; - css = ss+p->min_ss; - idx = cfs + css*image->width; + int idx; + int cfs, css; - image->dp[pi][fs+p->w*ss] = image->data[idx]; - image->bad[pi][fs+p->w*ss] = 0; + cfs = fs+p->min_fs; + css = ss+p->min_ss; + idx = cfs + css*image->width; - } + image->dp[pi][fs+p->w*ss] = image->data[idx]; + image->bad[pi][fs+p->w*ss] = 0; + } } } @@ -2188,8 +2292,7 @@ static int draw_detector(cairo_surface_t *surf, struct image *image, } -static int save_data_to_png(char *filename, struct detector *det, - int max_fs, int max_ss, double default_fill_value, +static int save_data_to_png(char *filename, struct enhanced_det *edet, double *data) { struct image im; @@ -2201,18 +2304,18 @@ static int save_data_to_png(char *filename, struct detector *det, cairo_status_t r; cairo_surface_t *surf; - im.data = malloc((max_fs+1)*(max_ss+1)*sizeof(float)); + im.data = malloc((edet->width)*(edet->height)*sizeof(float)); if ( im.data == NULL ) { ERROR("Failed to allocate memory to save data.\n"); return 1; } - im.det = det; - im.width = max_fs+1; - im.height = max_ss+1; + im.det = edet->det; + im.width = edet->width; + im.height = edet->height; im.flags = NULL; - for ( i=0; i<(max_fs+1)*(max_ss+1); i++) { - if ( data[i] == default_fill_value ) { + for ( i=0; i<(edet->width)*(edet->height); i++) { + if ( data[i] == -10000.0) { im.data[i] = 0.0; } else if ( data[i] > 1.0) { im.data[i] = 1.0; @@ -2264,61 +2367,9 @@ static int save_data_to_png(char *filename, struct detector *det, #endif /* HAVE_SAVE_TO_PNG */ -static void calculate_panel_correction(int di, int ip, int aw, - int *num_pix_displ, - struct rg_collection *connected, - struct connected_data *conn_data) -{ - struct panel *p; - int ifs, iss; - - p = connected->rigid_groups[di]->panels[ip]; - for (ifs=p->min_fs; ifsmax_fs+1; ifs++) { - for (iss=p->min_ss; issmax_ss+1; iss++) { - if ( num_pix_displ[ifs+aw*iss]>= - conn_data[di].num_peaks_per_pixel ) { - conn_data[di].n_peaks_in_conn++; - } - } - } - -} - - -static void compute_abs_displ(struct rg_collection *connected, - struct connected_data *conn_data, - int *num_pix_displ, - double dfv, int di, int ip, int aw, - double *displ_x, - double *displ_y, - double *displ_abs) -{ - struct panel *p; - int ifs, iss; - - if (conn_data[di].sh_x < dfv+1) return; - - p = connected->rigid_groups[di]->panels[ip]; - - for (ifs=p->min_fs; ifsmax_fs+1; ifs++) { - for (iss=p->min_ss; issmax_ss+1; iss++) { - if ( num_pix_displ[ifs+aw*iss]>= - conn_data[di].num_peaks_per_pixel ) { - displ_x[ifs+aw*iss] -= conn_data[di].sh_x; - displ_y[ifs+aw*iss] -= conn_data[di].sh_y; - displ_abs[ifs+aw*iss] = modulus2d( - displ_x[ifs+aw*iss], - displ_y[ifs+aw*iss] - ); - } else { - displ_abs[ifs+aw*iss] = dfv; - } - } - } -} - -int check_and_enforce_cspad_dist(struct detector *det, int enforce) +int check_and_enforce_cspad_dist(struct geoptimiser_params *gparams, + struct detector *det) { int np = 0; int num_errors_found = 0; @@ -2334,7 +2385,7 @@ int check_and_enforce_cspad_dist(struct detector *det, int enforce) struct panel *op = &det->panels[np+1]; dist2 = (( ep->cnx - op->cnx )*( ep->cnx - op->cnx ) + - ( ep->cny - op->cny )*( ep->cny - op->cny )); + ( ep->cny - op->cny )*( ep->cny - op->cny )); if ( dist2 > (dist_to_check+tol)*(dist_to_check+tol) || dist2 < (dist_to_check-tol)*(dist_to_check-tol) ) { @@ -2345,7 +2396,7 @@ int check_and_enforce_cspad_dist(struct detector *det, int enforce) "is outside acceptable margins.\n", ep->name, op->name); - if ( enforce ) { + if ( gparams->enforce_cspad_layout ) { double new_op_cx, new_op_cy; @@ -2368,7 +2419,7 @@ int check_and_enforce_cspad_dist(struct detector *det, int enforce) STATUS("Warning: relative orientation between panels " "%s and %s is incorrect.\n", ep->name, op->name); - if ( enforce ) { + if ( gparams->enforce_cspad_layout ) { STATUS("Enforcing relative orientation....\n"); @@ -2385,67 +2436,313 @@ int check_and_enforce_cspad_dist(struct detector *det, int enforce) } - } + } return num_errors_found; } +static struct connected_data *initialize_conn_data(struct rg_collection *quadrants, + struct rg_collection *connected) +{ + + struct connected_data *conn_data; + + conn_data = malloc(connected->n_rigid_groups* + sizeof(struct connected_data)); + + int di; + + for (di=0; din_rigid_groups; di++) { + + conn_data[di].num_quad = find_quad_for_connected( + connected->rigid_groups[di], + quadrants); + conn_data[di].cang = 0.0; + conn_data[di].cstr = 1.0; + conn_data[di].sh_x = -10000.0; + conn_data[di].sh_y = -10000.0; + conn_data[di].num_peaks_per_pixel = 1; + conn_data[di].name = connected->rigid_groups[di]->name; + conn_data[di].n_peaks_in_conn = 0; + } + + return conn_data; +} + + +static struct pixel_displ_list *initialize_pixel_displacement_list( + struct enhanced_det *edet) +{ + + struct pixel_displ_list *pix_displ_list; + int ipx; + + pix_displ_list = malloc(sizeof(struct pixel_displ_list)); + + pix_displ_list->pix_displ_list = calloc(edet->num_pix, + sizeof(struct single_pixel_displ)); + if ( pix_displ_list->pix_displ_list == NULL ) { + ERROR("Error allocating memory for pixel displacement data.\n"); + free(pix_displ_list); + return NULL; + } + pix_displ_list->curr_pix_displ = calloc(edet->num_pix, + sizeof(struct single_pixel_displ*)); + if ( pix_displ_list->curr_pix_displ == NULL ) { + ERROR("Error allocating memory for pixel displacement data.\n"); + free(pix_displ_list->pix_displ_list); + free(pix_displ_list); + return NULL; + } + pix_displ_list->num_pix_displ = calloc(edet->num_pix, sizeof(int)); + if ( pix_displ_list->num_pix_displ == NULL ) { + ERROR("Error allocating memory for pixel displacement data.\n"); + free(pix_displ_list->curr_pix_displ); + free(pix_displ_list->pix_displ_list); + free(pix_displ_list); + return NULL; + } + + for ( ipx=0; ipxnum_pix; ipx++ ) { + pix_displ_list->pix_displ_list[ipx].dfs = -10000.0; + pix_displ_list->pix_displ_list[ipx].dss = -10000.0; + pix_displ_list->pix_displ_list[ipx].ne = NULL; + pix_displ_list->curr_pix_displ[ipx] = &pix_displ_list->pix_displ_list[ipx]; + pix_displ_list->num_pix_displ[ipx] = 0; + } + + return pix_displ_list; +} + + +static void free_pixel_displacement_list( + struct pixel_displ_list *pix_displ_list, + struct enhanced_det *edet) +{ + int i; + struct single_pixel_displ *curr = NULL; + struct single_pixel_displ *next = NULL; + + for ( i=0; inum_pix; i++ ) { + + curr = &pix_displ_list->pix_displ_list[i]; + + if ( curr->ne != NULL ) { + curr = curr->ne; + while ( curr != NULL ) { + next = curr->ne; + free(curr); + curr = next; + } + } + } + + free(pix_displ_list->curr_pix_displ); + free(pix_displ_list->pix_displ_list); + free(pix_displ_list->num_pix_displ); + free(pix_displ_list); +} + + +static struct avg_displacements *initialize_average_displacement( + struct enhanced_det *edet) + +{ + static struct avg_displacements *avg_displ; + + avg_displ = malloc(sizeof(struct avg_displacements)); + avg_displ->displ_x = calloc(edet->num_pix, sizeof(double)); + if ( avg_displ->displ_x == NULL ) { + ERROR("Error allocating memory for pixel properties.\n"); + return NULL; + } + avg_displ->displ_y = calloc(edet->num_pix, sizeof(double)); + if ( avg_displ->displ_y == NULL ) { + ERROR("Error allocating memory for pixel properties.\n"); + free(avg_displ->displ_x); + free(avg_displ); + return NULL; + } + avg_displ->displ_abs = calloc(edet->num_pix, sizeof(double)); + if ( avg_displ->displ_abs == NULL ) { + ERROR("Error allocating memory for pixel properties.\n"); + free(avg_displ->displ_x); + free(avg_displ->displ_y); + free(avg_displ); + return NULL; + } + + return avg_displ; +} + + +static void free_avg_displacement(struct avg_displacements * avg_displ) +{ + free(avg_displ->displ_x); + free(avg_displ->displ_y); + free(avg_displ->displ_abs); + free(avg_displ); +} + + +static void free_pattern_list(struct pattern_list *pattern_list) +{ + + int pti; + + for ( pti=0; ptin_patterns; pti++ ) { + int nuc; + + image_feature_list_free(pattern_list->patterns[pti]->im_list); + reflist_free(pattern_list->patterns[pti]->ref_list); + for ( nuc=0; nucpatterns[pti]->n_unit_cells; + nuc++) { + cell_free(pattern_list->patterns[pti]->unit_cells[nuc]); + } + free(pattern_list->patterns[pti]->filename); + free(pattern_list->patterns[pti]); + } + free(pattern_list); +} + + +static int *extract_num_pix_free_displ_list(struct pixel_displ_list *pix_displ_list, + struct enhanced_det *edet) +{ + int *num_pix; + + int i; + struct single_pixel_displ *curr = NULL; + struct single_pixel_displ *next = NULL; + + for ( i=0; inum_pix; i++ ) { + + curr = &pix_displ_list->pix_displ_list[i]; + + if ( curr->ne != NULL ) { + curr = curr->ne; + while ( curr != NULL ) { + next = curr->ne; + free(curr); + curr = next; + } + } + } + + num_pix = pix_displ_list->num_pix_displ; + + free(pix_displ_list->curr_pix_displ); + free(pix_displ_list->pix_displ_list); + free(pix_displ_list); + + return num_pix; +} + + +static void recompute_panel_avg_displ(struct rg_collection *connected, + struct connected_data *conn_data, + int *num_pix_displ, int di, int ip, + struct enhanced_det *edet, + struct avg_displacements *avg_displ) +{ + struct panel *p; + int ifs, iss; + + if (conn_data[di].sh_x < -9999.0) return; + + p = connected->rigid_groups[di]->panels[ip]; + + for ( ifs=p->min_fs; ifsmax_fs+1; ifs++ ) { + for ( iss=p->min_ss; issmax_ss+1; iss++ ) { + + int pix_index = ifs+edet->width*iss; + + if ( num_pix_displ[pix_index]>= + conn_data[di].num_peaks_per_pixel ) { + avg_displ->displ_x[pix_index] -= conn_data[di].sh_x; + avg_displ->displ_y[pix_index] -= conn_data[di].sh_y; + avg_displ->displ_abs[pix_index] = modulus2d( + avg_displ->displ_x[pix_index], + avg_displ->displ_y[pix_index] + ); + } else { + avg_displ->displ_abs[pix_index] = -10000.0; + } + } + } +} + + +void recompute_avg_displ(struct rg_collection *connected, + struct connected_data *conn_data, + int *num_pix_displ, + struct enhanced_det *edet, + struct avg_displacements *avg_displ) +{ + + int di, ip; + + for ( di=0;din_rigid_groups;di++ ) { + for (ip=0; iprigid_groups[di]->n_panels; ip++) { + + recompute_panel_avg_displ(connected, conn_data, + num_pix_displ, di, ip, edet, + avg_displ); + + } + } +} + -int optimize_geometry(char *infile, char *outfile, char *geometry_filename, - struct detector *det, struct rg_collection* quadrants, - struct rg_collection* connected, - int min_num_peaks_per_pixel, int min_num_peaks_per_panel, - int only_best_distance, int nostretch, - int individual_coffset, int error_maps, - int enforce_cspad_layout, int no_cspad, - double max_peak_dist, const char *command_line) +int optimize_geometry(struct geoptimiser_params *gparams, + struct detector *det, + struct rg_collection *quadrants, + struct rg_collection *connected) { - int num_pix_in_slab; int max_fs = 0; int max_ss = 0; - int aw = 0; - int pi, di, ip, pti; - int ret1, ret2; + int pi; int ret; int write_ret; int maybe_cspad = 0; - int *num_pix_displ; double res_sum; - double istep; + double avg_res; double clen_to_use; - double man_stretching_coeff = 0.0; - double avc[6] = {0.,0.,0.,0.,0.,0.}; - double dfv = -10000.0; + // for angles and stretch calculation use - // only pixels which are distco*size_panel - // away - double dist_coeff_ang_str = 0.2; - double *displ_x; - double *displ_y; - double *displ_abs; - double totalError; - double min_braggp_dist; + // only pixels which are distco*size_panel + // away + double dist_coeff_for_rot_str = 0.2; + double total_error; - double* slab_to_x; - double* slab_to_y; - double* recomputed_slab_to_x; - double* recomputed_slab_to_y; - double stretch_coeff = 1; + struct pixel_maps *pixel_maps; + struct pixel_maps *recomputed_pixel_maps; + double stretch_coeff = 1.0; - struct single_pix_displ *all_pix_displ; - struct single_pix_displ **curr_pix_displ; struct connected_data *conn_data = NULL; struct pattern_list *pattern_list; + struct cell_params *cell_params; + struct enhanced_det edet; + struct avg_displacements *avg_displ; + struct pixel_displ_list *pix_displ_list; - if ( nostretch ) man_stretching_coeff = 1.0; + STATUS("Maximum distance between peaks: %0.1f pixels.\n", + gparams->max_peak_dist); - if ( det->n_panels == 64 ) { + STATUS("Minimum number of measurements for a pixel to be included in the " + "refinement: %i\n", + gparams->min_num_peaks_per_pix); + STATUS("Minimum number of measurements for connected group for accurate " + "estimation of position/orientation: %i\n", + gparams->min_num_pix_per_conn_group); + + if (det->n_panels == 64 ) { maybe_cspad = 1; } - if ( maybe_cspad && !no_cspad ) { + if ( maybe_cspad && !gparams->no_cspad ) { int num_errors = 0; @@ -2454,14 +2751,13 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, "connected ASICS.\n"); STATUS("If the detector is not a CSPAD, please rerun the " "program with the --no-cspad option.\n"); - if ( enforce_cspad_layout ) { + if ( gparams->enforce_cspad_layout ) { STATUS("Enforcing CSPAD layout...\n"); } - num_errors = check_and_enforce_cspad_dist(det, - enforce_cspad_layout); + num_errors = check_and_enforce_cspad_dist(gparams, det); - if ( enforce_cspad_layout ) { + if ( gparams->enforce_cspad_layout ) { int geom_wr; @@ -2469,9 +2765,10 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, "Please restart geometry optimization using the " "optimized geometry from this run as input geometry " "file.\n"); - geom_wr = write_detector_geometry_2(geometry_filename, - outfile, det, - command_line, 1); + 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; @@ -2480,9 +2777,9 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, return 0; } - if ( !enforce_cspad_layout && num_errors > 0 ) { + if ( !gparams->enforce_cspad_layout && num_errors > 0 ) { ERROR("Relative distance and orientation of connected " - "panels do not respect the CSPAD layout.\n" + "ASICS do not respect the CSPAD layout.\n" "Geometry optimization cannot continue.\n" "Please rerun the program with the " "--enforce-cspad-layout option.\n"); @@ -2490,13 +2787,17 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, } } - pattern_list = read_patterns_from_steam_file(infile, det); + pattern_list = read_patterns_from_steam_file(gparams->infile, det); if ( pattern_list->n_patterns < 1 ) { ERROR("Error reading stream file\n"); return 1; } - compute_avg_cell_parameters(pattern_list, avc); + cell_params = compute_avg_cell_parameters(pattern_list); + if ( cell_params == NULL ) { + free(pattern_list); + return 1; + } res_sum = 0; for ( pi=0; pin_panels; pi++ ) { @@ -2509,362 +2810,224 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, } res_sum += det->panels[pi].res; } + avg_res = res_sum/det->n_panels; - istep = res_sum/det->n_panels; - - aw = max_fs+1; + edet.det = det; + edet.width = max_fs+1; + edet.height = max_ss+1; + edet.num_pix = (max_fs+1)*(max_ss+1); - clen_to_use = compute_clen_to_use(pattern_list, istep, avc, - max_peak_dist, - only_best_distance, - &min_braggp_dist); + clen_to_use = pick_clen_to_use(gparams, pattern_list, avg_res, + cell_params); if ( clen_to_use == -1.0 ) return 1; - if ( max_peak_dist == -1.0 ) { - STATUS ("Setting maximum distance for peak search to 0.5*" - "minimum Bragg peak distance from the data\n"); - max_peak_dist = 0.5*min_braggp_dist; - } - STATUS("Maximum distance for peak search: %0.1f pixels.\n", max_peak_dist); - if ( min_braggp_dist<1.2*max_peak_dist ) { - STATUS("WARNING: The maximum distance for peak search " - "is too small compared to minimum Bragg peak " - "distance: %0.1f < 1.2*%0.1f pixels.\n", - min_braggp_dist, - max_peak_dist); - } - - STATUS("Minimum number of measurements for a pixel to be included in the " - "refinement: %i\n", - min_num_peaks_per_pixel); - STATUS("Minimum number of measurements per connected group for accurate " - "estimation of position/orientation: %i\n", - min_num_peaks_per_panel); - - num_pix_in_slab = (max_fs+1)*(max_ss+1); - displ_x = calloc(num_pix_in_slab,sizeof(double)); - if ( displ_x == NULL ) { - ERROR("Error allocating memory for pixel properties.\n"); - return 1; - } - displ_y = calloc(num_pix_in_slab,sizeof(double)); - if ( displ_y == NULL ) { - ERROR("Error allocating memory for pixel properties.\n"); - free(displ_x); - return 1; - } - displ_abs = calloc(num_pix_in_slab,sizeof(double)); - if ( displ_abs == NULL ) { - ERROR("Error allocating memory for pixel properties.\n"); - free(displ_x); - free(displ_y); + avg_displ = initialize_average_displacement(&edet); + if ( avg_displ == NULL ) { + free(cell_params); + free(pattern_list); return 1; } - slab_to_x = malloc(num_pix_in_slab*sizeof(double)); - slab_to_y = malloc(num_pix_in_slab*sizeof(double)); - if ( slab_to_x == NULL ) { + pixel_maps = initialize_pixel_maps(&edet); + if ( pixel_maps == NULL ) { ERROR("Failed to allocate memory for pixel maps.\n"); - free(displ_x); - free(displ_y); - free(displ_abs); - return 1; - } - slab_to_y = malloc(num_pix_in_slab*sizeof(double)); - if ( slab_to_y == NULL ) { - ERROR("Failed to allocate memory for pixel maps.\n"); - free(displ_x); - free(displ_y); - free(displ_abs); - free(slab_to_x); + free_avg_displacement(avg_displ); + free(cell_params); + free(pattern_list); return 1; } - fill_coordinate_matrices(det, aw, slab_to_x, slab_to_y); - - all_pix_displ = calloc(num_pix_in_slab, - sizeof(struct single_pix_displ)); - if ( all_pix_displ == NULL ) { - ERROR("Error allocating memory for connected structure data.\n"); - free(displ_x); - free(displ_y); - free(displ_abs); - free(slab_to_x); - free(slab_to_y); - return 1; - } - curr_pix_displ = calloc(num_pix_in_slab, - sizeof(struct single_pix_displ*)); - if ( curr_pix_displ == NULL ) { - ERROR("Error allocating memory for connected structure data.\n"); - free(displ_x); - free(displ_y); - free(displ_abs); - free(slab_to_x); - free(slab_to_y); - free(all_pix_displ); - return 1; - } - num_pix_displ = calloc(num_pix_in_slab, sizeof(int)); - if ( num_pix_displ == NULL ) { + pix_displ_list = initialize_pixel_displacement_list(&edet); + if ( pix_displ_list == NULL ) { ERROR("Error allocating memory for connected structure data.\n"); - free(displ_x); - free(displ_y); - free(displ_abs); - free(slab_to_x); - free(slab_to_y); - free(all_pix_displ); - free(curr_pix_displ); + free_avg_displacement(avg_displ); + free_pixel_maps(pixel_maps); + free(pattern_list); return 1; } - conn_data = malloc(connected->n_rigid_groups* - sizeof(struct connected_data)); + conn_data = initialize_conn_data(quadrants, connected); if ( conn_data == NULL ) { ERROR("Error allocating memory for connected structure data.\n"); - free(displ_x); - free(displ_y); - free(displ_abs); - free(slab_to_x); - free(slab_to_y); - free(all_pix_displ); - free(curr_pix_displ); - free(num_pix_displ); + free_avg_displacement(avg_displ); + free_pixel_maps(pixel_maps); + free_pixel_displacement_list(pix_displ_list, &edet); + free(pattern_list); return 1; } - STATUS("Computing pixel statistics.\n"); - ret = compute_pixel_statistics(pattern_list, det, connected, quadrants, - num_pix_in_slab, max_peak_dist, - aw, dfv, - min_num_peaks_per_pixel, - min_num_peaks_per_panel, - only_best_distance, - clen_to_use, slab_to_x, - slab_to_y, conn_data, - displ_x, displ_y, displ_abs, - all_pix_displ, - curr_pix_displ, - num_pix_displ); + + ret = compute_pixel_displacements(pattern_list, &edet, connected, + gparams, clen_to_use, conn_data, + avg_displ, pix_displ_list); if ( ret != 0 ) { - free(displ_x); - free(displ_y); - free(displ_abs); - free(slab_to_x); - free(slab_to_y); - free_all_curr_pix_displ(all_pix_displ, curr_pix_displ, - num_pix_in_slab); - free(num_pix_displ); + free_avg_displacement(avg_displ); + free_pixel_maps(pixel_maps); + free_pixel_displacement_list(pix_displ_list, &edet); free(conn_data); + free(pattern_list); return 1; } - free_all_curr_pix_displ(all_pix_displ, curr_pix_displ, num_pix_in_slab); - for ( pti=0; ptin_patterns; pti++ ) { - int nuc; + free_pattern_list(pattern_list); - image_feature_list_free(pattern_list->patterns[pti]->im_list); - reflist_free(pattern_list->patterns[pti]->ref_list); - for ( nuc=0; nucpatterns[pti]->n_unit_cells; - nuc++) { - cell_free(pattern_list->patterns[pti]->unit_cells[nuc]); - } - free(pattern_list->patterns[pti]->filename); - free(pattern_list->patterns[pti]); + adjust_min_peaks_per_conn(&edet, connected, gparams, conn_data, + pix_displ_list); + + ret = compute_avg_displacements(&edet, connected, pix_displ_list, + conn_data, avg_displ); + if ( ret != 0 ) { + free_avg_displacement(avg_displ); + free_pixel_maps(pixel_maps); + free_pixel_displacement_list(pix_displ_list, &edet); + free(conn_data); + return 1; } - free(pattern_list); - if ( error_maps ) { - STATUS("Saving pixel position errors before correction.\n"); + num_pix_displ = extract_num_pix_free_displ_list(pix_displ_list, &edet); + + STATUS("Computing error before correction.\n"); + total_error = compute_error(connected, &edet, conn_data, + num_pix_displ, avg_displ->displ_abs); + + STATUS("Detector-wide error before correction = %0.4f pixels.\n", + total_error); + + if ( gparams->error_maps ) { + STATUS("Saving error map before correction.\n"); #ifdef HAVE_SAVE_TO_PNG - ret1 = save_data_to_png("error_map_before.png", det, max_fs, max_ss, - dfv, displ_abs); - if ( ret1!=0 ) { + ret = save_data_to_png("error_map_before.png", &edet, + avg_displ->displ_abs); + if ( ret != 0 ) { ERROR("Error while writing data to file.\n"); + + free_avg_displacement(avg_displ); + free_pixel_maps(pixel_maps); free(conn_data); - free(displ_x); - free(displ_y); - free(displ_abs); free(num_pix_displ); - free(slab_to_x); - free(slab_to_y); return 1; } #else /* HAVE_SAVE_TO_PNG */ - STATUS("ERROR: geoptimiser was compiled without GTK and cairo " + ERROR("WARNING: geoptimiser was compiled without GTK and cairo " "support. Error maps will not be saved.\n"); #endif /* HAVE_SAVE_TO_PNG */ } - STATUS("Computing initial error.\n"); - totalError = compute_error(connected, aw, conn_data, - num_pix_displ, displ_abs); - - STATUS("The total initial error = %0.4f pixels.\n", - totalError); - STATUS("Now calculating corrections.\n"); - - for ( di=0;din_rigid_groups;di++ ) { - - conn_data[di].n_peaks_in_conn = 0; - - for (ip=0; iprigid_groups[di]->n_panels; ip++) { - - calculate_panel_correction(di, ip, aw, num_pix_displ, - connected, conn_data); - - } - } - - STATUS("Calculating angles and elongations.\n"); - - ret = compute_angles_and_stretch(connected, conn_data, - num_pix_displ, - slab_to_x, slab_to_y, - displ_x, displ_y, - aw, - min_num_peaks_per_panel, - dist_coeff_ang_str, - min_num_peaks_per_pixel, - man_stretching_coeff, &stretch_coeff); - if ( ret != 0 ) { + stretch_coeff = compute_rotation_and_stretch(connected, conn_data, + &edet, num_pix_displ, + pixel_maps, avg_displ, + dist_coeff_for_rot_str, + gparams); + if ( stretch_coeff < 0.0 ) { + free_avg_displacement(avg_displ); + free_pixel_maps(pixel_maps); free(conn_data); - free(displ_x); - free(displ_y); - free(displ_abs); free(num_pix_displ); - free(slab_to_x); - free(slab_to_y); return 1; } - ret = correct_empty_panels(quadrants, connected, min_num_peaks_per_panel, - conn_data); + ret = compute_rot_stretch_for_empty_panels(quadrants, connected, + gparams, conn_data); if ( ret != 0 ) { + free_avg_displacement(avg_displ); + free_pixel_maps(pixel_maps); free(conn_data); - free(displ_x); - free(displ_y); - free(displ_abs); free(num_pix_displ); - free(slab_to_x); - free(slab_to_y); return 1; } - correct_angle_and_stretch(connected, det, conn_data, clen_to_use, - stretch_coeff, individual_coffset); - shift_panels(connected, conn_data); + correct_rotation_and_stretch(connected, edet.det, conn_data, + clen_to_use, stretch_coeff, + gparams); - recomputed_slab_to_x = malloc(num_pix_in_slab*sizeof(double)); - if ( recomputed_slab_to_x == NULL ) { - ERROR("Failed to allocate memory for pixel maps.\n"); - free(displ_x); - free(displ_y); - free(displ_abs); - free(slab_to_x); - free(slab_to_y); - free(num_pix_displ); + recomputed_pixel_maps = initialize_pixel_maps(&edet); + if ( recomputed_pixel_maps == NULL ) { + free_avg_displacement(avg_displ); + free_pixel_maps(pixel_maps); free(conn_data); - return 1; - } - recomputed_slab_to_y = malloc(num_pix_in_slab*sizeof(double)); - if ( recomputed_slab_to_y == NULL ) { - ERROR("Failed to allocate memory for pixel maps.\n"); - free(displ_x); - free(displ_y); - free(displ_abs); - free(slab_to_x); - free(slab_to_y); free(num_pix_displ); - free(conn_data); - free(recomputed_slab_to_x); return 1; } - fill_coordinate_matrices(det, aw, recomputed_slab_to_x, - recomputed_slab_to_y); - - recompute_differences(connected, slab_to_x, slab_to_y, - recomputed_slab_to_x, - recomputed_slab_to_y, conn_data, - stretch_coeff, aw, displ_x, displ_y, - num_pix_displ); - - ret = compute_shifts(connected, conn_data, num_pix_displ, aw, - min_num_peaks_per_panel, dfv, - max_peak_dist, displ_x, displ_y ); - - if ( ret != 0 ) return 1; - - compute_shifts_for_empty_panels(quadrants, connected, conn_data, - min_num_peaks_per_panel); - - for ( di=0;din_rigid_groups;di++ ) { - for (ip=0; iprigid_groups[di]->n_panels; ip++) { - - compute_abs_displ(connected, conn_data, - num_pix_displ, dfv, di, ip, aw, - displ_x, displ_y, displ_abs); + adjust_displ_for_stretch(connected, pixel_maps, + recomputed_pixel_maps, + conn_data, + stretch_coeff, &edet, avg_displ, + num_pix_displ); - } + ret = compute_shift(connected, conn_data, num_pix_displ, &edet, + gparams, avg_displ); + if ( ret != 0 ) { + free_avg_displacement(avg_displ); + free_pixel_maps(pixel_maps); + free(conn_data); + free(num_pix_displ); + free_pixel_maps(recomputed_pixel_maps); + return 1; } - correct_shifts(connected, conn_data, dfv, clen_to_use); + compute_shift_for_empty_panels(quadrants, connected, conn_data, + gparams); + + correct_shift(connected, conn_data, clen_to_use); - if ( error_maps ) { + recompute_avg_displ(connected, conn_data, + num_pix_displ, &edet, + avg_displ); + if ( gparams->error_maps ) { #ifdef HAVE_SAVE_TO_PNG - STATUS("Saving pixel position errors after correction.\n"); - ret2 = save_data_to_png("error_map_after.png", det, max_fs, max_ss, - dfv, displ_abs); - if ( ret2 !=0 ) { + STATUS("Saving error map after correction.\n"); + + ret = save_data_to_png("error_map_after.png", &edet, + avg_displ->displ_abs); + if ( ret !=0 ) { ERROR("Error while writing data to file.\n"); + + free_avg_displacement(avg_displ); + free_pixel_maps(pixel_maps); free(conn_data); - free(displ_x); - free(displ_y); - free(displ_abs); free(num_pix_displ); - free(slab_to_x); - free(slab_to_y); - free(recomputed_slab_to_x); - free(recomputed_slab_to_y); + free_pixel_maps(recomputed_pixel_maps); return 1; } #else /* HAVE_SAVE_TO_PNG */ STATUS("ERROR: geoptimiser was compiled without GTK and cairo " - "support. Error maps will not be saved.\n"); + "support.\n Error maps will not be saved.\n"); #endif /* HAVE_SAVE_TO_PNG */ } - STATUS("Computing final error.\n"); - totalError = compute_error(connected, aw, conn_data, num_pix_displ, - displ_abs); + STATUS("Computing errors after correction.\n"); + total_error = compute_error(connected, &edet, conn_data, + num_pix_displ, avg_displ->displ_abs); - STATUS("The total final error = %0.4f pixels.\n",totalError); + STATUS("Detector-wide error after correction = %0.4f pixels.\n", + total_error); - write_ret = write_detector_geometry_2(geometry_filename, outfile, det, - command_line, 1); + write_ret = write_detector_geometry_2(gparams->geometry_filename, + gparams->outfile, edet.det, + gparams->command_line, 1); if ( write_ret != 0 ) { ERROR("Error in writing output geometry file.\n"); return 1; } STATUS("All done!\n"); - if ( error_maps ) { + if ( gparams->error_maps ) { #ifdef HAVE_SAVE_TO_PNG @@ -2874,65 +3037,71 @@ int optimize_geometry(char *infile, char *outfile, char *geometry_filename, #endif /* HAVE_SAVE_TO_PNG */ } + free_avg_displacement(avg_displ); + free_pixel_maps(pixel_maps); free(conn_data); - free(displ_x); - free(displ_y); - free(displ_abs); free(num_pix_displ); - free(slab_to_x); - free(slab_to_y); - free(recomputed_slab_to_x); - free(recomputed_slab_to_y); + free_pixel_maps(recomputed_pixel_maps); return 0; } + int main(int argc, char *argv[]) { int c, i; int ret_val; char buffer[256]; char command_line[1024]; - char *outfile = NULL; - char *infile = NULL; - char *geometry_filename = NULL; + char *quadrant_coll_name = NULL; char *connected_coll_name = NULL; - int min_num_peaks_per_pixel = 3; - int min_num_peaks_per_panel = 100; - int only_best_distance = 0; - int enforce_cspad_layout = 0; - int nostretch = 0; - int individual_coffset = 0; - int no_cspad = 0; - int error_maps = 1; - double max_peak_dist = -1.0; + struct geoptimiser_params *gparams; struct detector *det = NULL; struct rg_collection *quadrants; struct rg_collection *connected; struct beam_params beam; + gparams = malloc(sizeof(struct geoptimiser_params)); + + gparams->outfile = NULL; + gparams->infile = NULL; + gparams->geometry_filename = NULL; + gparams->min_num_peaks_per_pix = 3; + gparams->min_num_pix_per_conn_group = 100; + gparams->only_best_distance = 0; + gparams->enforce_cspad_layout = 0; + gparams->nostretch = 0; + gparams->individual_coffset = 0; + gparams->no_cspad = 0; + gparams->error_maps = 1; + gparams->max_peak_dist = 4.0; + const struct option longopts[] = { - /* Options with long and short versions */ - {"help", 0, NULL, 'h'}, - {"version", 0, NULL, 10 }, - {"input", 1, NULL, 'i'}, - {"output", 1, NULL, 'o'}, - {"geometry", 1, NULL, 'g'}, - {"quadrants", 1, NULL, 'q'}, - {"connected", 1, NULL, 'c'}, - {"min-num-peaks-per-pixel", 1, NULL, 'x'}, - {"min-num-peaks-per-panel", 1, NULL, 'p'}, - {"most-few-clen", 0, NULL, 'l'}, - {"max-peak-dist", 1, NULL, 'm'}, - {"individual-dist-offset", 0, NULL, 's'}, - - /* Long-only options with no arguments */ - {"no-stretch", 0, &nostretch, 1}, - {"no-error-maps", 0, &error_maps, 0}, - {"enforce-cspad-layout", 0, &enforce_cspad_layout, 1}, - {"no-cspad", 0, &no_cspad, 1}, + /* Options with long and short versions */ + {"help", 0, NULL, 'h'}, + {"version", 0, NULL, 10 }, + {"input", 1, NULL, 'i'}, + {"output", 1, NULL, 'o'}, + {"geometry", 1, NULL, 'g'}, + {"quadrants", 1, NULL, 'q'}, + {"connected", 1, NULL, 'c'}, + {"min-num-peaks-per-pixel", 1, NULL, 'x'}, + {"min-num-pixels-per-conn-group", 1, NULL, 'p'}, + {"most-few-clen", 0, NULL, 'l'}, + {"max-peak-dist", 1, NULL, 'm'}, + {"individual-dist-offset", 0, NULL, 's'}, + + /* Long-only options with no arguments */ + {"no-stretch", 0, &gparams->nostretch, 1}, + {"no-error-maps", 0, &gparams->error_maps, 0}, + {"enforce-cspad-layout", 0, &gparams->enforce_cspad_layout, 1}, + {"no-cspad", 0, &gparams->no_cspad, 1}, + + /* Long-only options with arguments */ + {"min-num-peaks-per-panel", 1, NULL, 11}, + {0, 0, NULL, 0} }; @@ -2953,16 +3122,17 @@ int main(int argc, char *argv[]) return 0; case 'o' : - outfile = strdup(optarg); + gparams->outfile = strdup(optarg); break; case 'i' : - infile = strdup(optarg); + gparams->infile = strdup(optarg); break; case 'g' : - geometry_filename = strdup(optarg); - det = get_detector_geometry(geometry_filename, &beam); + gparams->geometry_filename = strdup(optarg); + det = get_detector_geometry(gparams->geometry_filename, + &beam); if ( det == NULL ) { ERROR("Failed to read detector geometry from " "'%s'\n", optarg); @@ -2979,39 +3149,48 @@ int main(int argc, char *argv[]) break; case 'x' : - min_num_peaks_per_pixel = atoi(optarg); + gparams->min_num_peaks_per_pix = atoi(optarg); break; case 'p' : - min_num_peaks_per_panel = atoi(optarg); + gparams->min_num_pix_per_conn_group = atoi(optarg); break; case 'l' : - only_best_distance = 1; + gparams->only_best_distance = 1; break; case 'm' : - max_peak_dist = strtof(optarg, NULL); + gparams->max_peak_dist = strtof(optarg, NULL); break; case 's' : - individual_coffset = 1; + gparams->individual_coffset = 1; + break; + + case 11: + ERROR("WARNING: The --min-num-peaks-per-panel option has been " + "renamed to --min-num-pixels-per-conn-group. The " + "old option has been deprecated and will " + "soon be removed. It is currently remapped to " + "--min-num-pixels-per-conn-group\n"); + gparams->min_num_pix_per_conn_group = atoi(optarg); break; } } - if ( geometry_filename == NULL ) { + if ( gparams->geometry_filename == NULL ) { ERROR("You must provide a geometry to optimize.\n"); return 1; } - if ( infile == NULL ) { + if ( gparams->infile == NULL ) { ERROR("You must provide an input stream file.\n"); return 1; } - if ( outfile == NULL ) { + if ( gparams->outfile == NULL ) { ERROR("You must provide an output filename.\n"); return 1; } @@ -3052,12 +3231,7 @@ int main(int argc, char *argv[]) } g_type_init(); - ret_val = optimize_geometry(infile, outfile, geometry_filename, det, - quadrants, connected, min_num_peaks_per_pixel, - min_num_peaks_per_panel, only_best_distance, - nostretch, individual_coffset, error_maps, - enforce_cspad_layout, no_cspad, - max_peak_dist, command_line); + ret_val = optimize_geometry(gparams, det, quadrants, connected); return ret_val; } -- cgit v1.2.3