diff options
| -rw-r--r-- | libcrystfel/src/image.h | 3 | ||||
| -rw-r--r-- | src/geoptimiser.c | 1782 |
2 files changed, 646 insertions, 1139 deletions
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h index 89c01785..9fd9b495 100644 --- a/libcrystfel/src/image.h +++ b/libcrystfel/src/image.h @@ -174,6 +174,9 @@ struct image { const struct copy_hdf5_field *copyme; struct stuff_from_stream *stuff_from_stream; + double avg_clen; /* Average camera length extracted + * from stuff_from_stream */ + int id; /* ID number of the thread * handling this image */ int serial; /* Monotonically ascending serial diff --git a/src/geoptimiser.c b/src/geoptimiser.c index 5bb046a4..0a6885d0 100644 --- a/src/geoptimiser.c +++ b/src/geoptimiser.c @@ -3,13 +3,13 @@ * * Refine detector geometry * - * Copyright © 2014-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2014-2016 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: * 2014-2015 Oleksandr Yefanov * 2014-2015 Valerio Mariani - * 2014-2015 Thomas White <taw@physics.org> + * 2014-2016 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -51,13 +51,13 @@ #endif /* HAVE_GTK */ #endif /* HAVE_CAIRO */ -#include <detector.h> -#include <stream.h> -#include <version.h> -#include <crystal.h> -#include <image.h> -#include <utils.h> -#include <render.h> +#include "detector.h" +#include "stream.h" +#include "version.h" +#include "crystal.h" +#include "image.h" +#include "utils.h" +#include "render.h" #include "hdfsee-render.h" @@ -133,36 +133,11 @@ struct connected_data }; -struct pattern { - ImageFeatureList *im_list; - RefList *ref_list; - double clen; - UnitCell **unit_cells; - int n_unit_cells; - double lambda; - char *filename; -}; - - -struct pattern_list { - struct pattern **patterns; - int n_patterns; -}; - - struct single_pixel_displ { double dfs; double dss; - 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 single_pixel_displ *ne; }; @@ -173,24 +148,6 @@ struct connected_stretch_and_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 { @@ -207,302 +164,178 @@ struct avg_shift int *aver_num_sh; }; -struct pixel_maps + +struct gpanel { - double *pix_to_x; - double *pix_to_y; -}; + struct panel *p; + /* Individual pixel displacements */ + struct single_pixel_displ *pix_displ_list; + struct single_pixel_displ **curr_pix_displ; + int *num_pix_displ; -struct enhanced_det -{ - struct detector *det; - int width; - int height; - int num_pix; + struct avg_shift avg_shift; + struct avg_rot_and_stretch avg_rot_and_stretch; + + /* Average displacements for each pixel */ + double *avg_displ_x; + double *avg_displ_y; + double *avg_displ_abs; }; -static void compute_x_y(struct detector *det, double fs, double ss, - double * x, double *y) +static void compute_x_y(double fs, double ss, struct panel *p, + double *x, double *y) { - struct panel *p; double xs, ys; - double dfs, dss; - - p = find_panel(det, fs, ss); - - dss = ss-p->min_ss; - dfs = fs-p->min_fs; - xs = dfs*p->fsx + dss*p->ssx; - ys = dfs*p->fsy + dss*p->ssy; + xs = fs*p->fsx + ss*p->ssx; + ys = fs*p->fsy + ss*p->ssy; *x = xs + p->cnx; *y = ys + p->cny; } -static Reflection *find_closest_reflection(RefList *rlist, - double fx, double fy, - struct detector *det, - double *d) +static Reflection *find_closest_reflection(struct image *image, + double fx, double fy, + double *d) { - double dmin = HUGE_VAL; Reflection *closest = NULL; - Reflection *refl; - RefListIterator *iter; - - - for ( refl = first_refl(rlist, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) { - double ds; - double rfs, rss; - double rx, ry; - - get_detector_pos(refl, &rfs, &rss); - - compute_x_y(det, rfs, rss, &rx, &ry); - - ds = distance(rx, ry, fx, fy); - - if ( ds < dmin ) { - dmin = ds; - closest = refl; - } - - } - - if ( dmin < HUGE_VAL ) { - *d = dmin; - return closest; - } + int i; - *d = +INFINITY; - return NULL; -} + for ( i=0; i<image->n_crystals; i++ ) { + Reflection *refl; + RefListIterator *iter; + RefList *rlist = crystal_get_reflections(image->crystals[i]); -static double compute_average_clen (struct detector *det, char **clen_from, - double *offset) -{ + for ( refl = first_refl(rlist, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double ds; + double rfs, rss; + double rx, ry; - int np, num_pan; - double sum_clen; + get_detector_pos(refl, &rfs, &rss); - sum_clen = 0; - num_pan = 0; + compute_x_y(rfs, rss, get_panel(refl), &rx, &ry); - for ( np=0; np<det->n_panels; np++ ) { + ds = distance(rx, ry, fx, fy); - struct panel p = det->panels[np]; + if ( ds < dmin ) { + dmin = ds; + closest = refl; + } - if ( p.clen_from != NULL ) { - *clen_from = strdup(p.clen_from); - *offset = p.coffset; - return -1; - } else { - sum_clen += p.clen+p.coffset; - num_pan += 1; } } - return sum_clen/num_pan; - + if ( closest == NULL ) { + *d = +INFINITY; + } else { + *d = dmin; + } + return closest; } -static double extract_f_from_stuff(const char *field_name, - struct stuff_from_stream* stuff) +static double get_average_clen(struct image *image) { int i; + struct stuff_from_stream *stuff = image->stuff_from_stream; - char field_name_plus_equal[256]; - snprintf(field_name_plus_equal, 256, "hdf5%s = ", field_name); + if ( stuff == NULL ) { + ERROR("No 'stuff' from stream!\n"); + return -1.0; + } for ( i=0; i<stuff->n_fields; i++ ) { - if ( strncmp(stuff->fields[i], field_name_plus_equal, - strlen(field_name_plus_equal)) == 0 ) { - return atof(stuff->fields[i]+ - strlen(field_name_plus_equal)); + if ( strncmp(stuff->fields[i], "average_camera_length = ", + 24) == 0 ) + { + return atof(stuff->fields[i]+24); } } - ERROR("Failed to recover camera length from stream file\n"); - return -1; + ERROR("Failed to recover average camera length from stream file\n"); + return -1.0; } -static struct pattern_list *read_patterns_from_steam_file(const char *infile, - struct detector *det) +static struct image *read_patterns_from_stream(const char *infile, + struct detector *det, int *n) { Stream *st; - - struct pattern_list *pattern_list; - - int max_patterns, n_chunks; - - n_chunks = 0; - max_patterns = 0; - - pattern_list = malloc(sizeof(struct pattern_list)); - if ( pattern_list == NULL ) { - ERROR("Failed to allocate memory for loaded patterns.\n"); - return NULL; - } - pattern_list->n_patterns =0; - - pattern_list->patterns = malloc(1024*sizeof(struct pattern*)); - if ( pattern_list->patterns == NULL ) { - ERROR("Failed to allocate memory for loaded patterns.\n"); - free(pattern_list); + struct image *images; + int n_chunks = 0; + int max_images = 1024; + int n_images = 0; + + images = malloc(max_images * sizeof(struct image)); + if ( images == NULL ) { + ERROR("Failed to allocate memory for images.\n"); return NULL; } - pattern_list->n_patterns = 0; - max_patterns = 1024; st = open_stream_for_read(infile); if ( st == NULL ) { ERROR("Failed to open input stream '%s'\n", infile); - free(pattern_list->patterns); - free(pattern_list); + free(images); return NULL; } do { - struct image cur; - int i; + images[n_images].det = det; - cur.det = det; - cur.stuff_from_stream = NULL; + if ( read_chunk_2(st, &images[n_images], + STREAM_READ_REFLECTIONS + | STREAM_READ_PEAKS + | STREAM_READ_UNITCELL) != 0 ) break; - if ( read_chunk_2(st, &cur, STREAM_READ_REFLECTIONS - | STREAM_READ_PEAKS | STREAM_READ_UNITCELL) != 0 ) { - break; - } + n_chunks++; /* Number of chunks processed */ - n_chunks +=1; + /* Reject if there are no crystals (not indexed) */ + if ( images[n_images].n_crystals == 0 ) continue; - if ( cur.n_crystals !=0 ) { + images[n_images].avg_clen = get_average_clen(&images[n_images]); - struct pattern *patn; - double avg_clen = 0.0; - double offset = 0.0; - char *clen_from; + n_images++; /* Number of images accepted */ - if ( pattern_list->n_patterns == max_patterns ) { + if ( n_images == max_images ) { - struct pattern **patterns_new; + struct image *images_new; - patterns_new = realloc(pattern_list->patterns, - (max_patterns+1024)* - sizeof(struct pattern *)); - if ( patterns_new == NULL ) { - ERROR("Failed to allocate " - "memory for loaded patterns.\n"); - free(pattern_list->patterns); - free(pattern_list); - return NULL; - } - - max_patterns += 1024; - pattern_list->patterns = patterns_new; - } - - patn = malloc(sizeof(struct pattern)); - if ( patn == NULL ) { - ERROR("Failed to allocate memory for loaded " + images_new = realloc(images, + (max_images+1024)*sizeof(struct image)); + if ( images_new == NULL ) { + ERROR("Failed to allocate memory for " "patterns.\n"); - free(pattern_list->patterns); - free(pattern_list); + free(images); return NULL; } - patn->filename = cur.filename; - patn->unit_cells = NULL; - patn->n_unit_cells = 0; - patn->im_list = cur.features; - patn->ref_list = reflist_new(); - - clen_from = NULL; - avg_clen = compute_average_clen(det, &clen_from, &offset); - if ( avg_clen == -1 ) { - avg_clen = extract_f_from_stuff(clen_from, - cur.stuff_from_stream)*1e-3; - avg_clen += offset; - } - patn->clen = avg_clen; - free(clen_from); - - patn->lambda = cur.lambda; - - for ( i=0; i<cur.n_crystals; i++ ) { - - RefList *crystal_reflist; - Reflection *refl; - RefListIterator *iter; - UnitCell *cell; - UnitCell **new_unit_cells; - - cell = crystal_get_cell(cur.crystals[i]); - - new_unit_cells = realloc(patn->unit_cells, - (patn->n_unit_cells+1)* - sizeof(UnitCell *)); - if ( new_unit_cells == NULL ) { - ERROR("Failed to allocate memory for " - "loaded patterns.\n"); - free(pattern_list->patterns); - free(pattern_list); - free(patn); - return NULL; - } - - new_unit_cells[patn->n_unit_cells] = cell; - patn->n_unit_cells++; - patn->unit_cells = new_unit_cells; - - crystal_reflist = crystal_get_reflections( - cur.crystals[i]); - - for ( refl = first_refl(crystal_reflist, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - Reflection *n; - int h, k, l; - - get_indices(refl, &h, &k, &l); - n = add_refl(patn->ref_list, h, k, l); - copy_data(n, refl); - } - - } - - pattern_list->patterns[pattern_list->n_patterns] = patn; - pattern_list->n_patterns++; - - if ( pattern_list->n_patterns%1000 == 0 ) { - STATUS("Loaded %i indexed patterns from %i " - "total patterns.\n", - pattern_list->n_patterns, ++n_chunks); - } + max_images += 1024; + images = images_new; + } + if ( n_chunks % 1000 == 0 ) { + STATUS("Loaded %i indexed patterns from %i " + "total patterns.\n", + n_images, n_chunks); } } while ( 1 ); close_stream(st); + *n = n_images; - STATUS("Found %d indexed patterns in file %s (from a total of %d).\n", - pattern_list->n_patterns, infile, n_chunks ); - - return pattern_list; + return images; } @@ -514,39 +347,42 @@ static struct rvec get_q_from_xyz(double rx, double ry, double dist, double l) double twotheta = atan2(r, dist); double az = atan2(ry, rx); - q.u = 1.0/(l*1e9) * sin(twotheta)*cos(az); - q.v = 1.0/(l*1e9) * sin(twotheta)*sin(az); - q.w = 1.0/(l*1e9) * (cos(twotheta) - 1.0); + q.u = 1.0/l * sin(twotheta)*cos(az); + q.v = 1.0/l * sin(twotheta)*sin(az); + q.w = 1.0/l * (cos(twotheta) - 1.0); return q; } -static struct cell_params *compute_avg_cell_parameters(struct pattern_list - *pattern_list) +static UnitCell *compute_avg_cell_parameters(struct image *images, int n) { int numavc; int j, i; double minc[6]; double maxc[6]; - double avg_cpar[6]; + double avg_cpar[6] = {0, 0, 0, 0, 0, 0}; + UnitCell *avg; - for (j=0; j<6; j++) { - minc[j] = 1e10; - maxc[j] = -1e10; + for ( j=0; j<6; j++ ) { + minc[j] = 1e100; + maxc[j] = -1e100; } + numavc = 0; - for (i=0; i<pattern_list->n_patterns; i++) { + for ( i=0; i<n; i++ ) { - struct pattern *ptn; + struct image *image; double cpar[6]; int j, cri; - ptn = pattern_list->patterns[i]; + image = &images[i]; - for ( cri=0; cri<ptn->n_unit_cells; cri++ ) { + for ( cri=0; cri<images->n_crystals; cri++ ) { - cell_get_parameters(ptn->unit_cells[cri], + UnitCell *cell = crystal_get_cell(image->crystals[cri]); + + cell_get_parameters(cell, &cpar[0], // a &cpar[1], // b &cpar[2], // c @@ -554,17 +390,10 @@ static struct cell_params *compute_avg_cell_parameters(struct pattern_list &cpar[4], // beta &cpar[5]); // gamma - cpar[0] *= 1e9; - cpar[1] *= 1e9; - cpar[2] *= 1e9; - cpar[3] = rad2deg(cpar[3]); - cpar[4] = rad2deg(cpar[4]); - cpar[5] = rad2deg(cpar[5]); - for ( j=0; j<6; j++ ) { avg_cpar[j] += cpar[j]; - if (cpar[j]<minc[j]) minc[j] = cpar[j]; - if (cpar[j]>maxc[j]) maxc[j] = cpar[j]; + if ( cpar[j]<minc[j] ) minc[j] = cpar[j]; + if ( cpar[j]>maxc[j] ) maxc[j] = cpar[j]; } numavc++; @@ -572,46 +401,44 @@ static struct cell_params *compute_avg_cell_parameters(struct pattern_list } - if ( numavc>0 ) { + if ( numavc > 0 ) { for ( j=0; j<6; j++ ) avg_cpar[j] /= numavc; } - struct cell_params *cparams = malloc(sizeof(struct cell_params)); + avg = cell_new(); - 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]; + cell_set_parameters(avg, avg_cpar[0], avg_cpar[1], avg_cpar[2], + avg_cpar[3], avg_cpar[4], avg_cpar[5]); - STATUS("Average cell coordinates:\n"); - STATUS("Average a, b, c (in nm): %6.3f, %6.3f, %6.3f\n", - cparams->a, cparams->b, cparams->c); + STATUS("Average cell parameters:\n"); + STATUS("Average a, b, c (in A): %6.3f, %6.3f, %6.3f\n", + avg_cpar[0]*1e10, avg_cpar[1]*1e10, avg_cpar[2]*1e10); 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]); + "\t%6.3f - %6.3f A,\n" + "\t%6.3f - %6.3f A,\n" + "\t%6.3f - %6.3f A\n", + minc[0]*1e10, maxc[0]*1e10, minc[1]*1e10, + maxc[1]*1e10, minc[2]*1e10, maxc[2]*1e10); STATUS("Average alpha,beta,gamma in degrees: %6.3f, %6.3f, %6.3f\n", - cparams->alpha, cparams->beta, cparams->gamma); + rad2deg(avg_cpar[3]), rad2deg(avg_cpar[4]), rad2deg(avg_cpar[5])); 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]); + rad2deg(minc[3]), rad2deg(maxc[3]), + rad2deg(minc[4]), rad2deg(maxc[4]), + rad2deg(minc[5]), rad2deg(maxc[5])); - return cparams; + return avg; } static double pick_clen_to_use(struct geoptimiser_params *gparams, - struct pattern_list *pattern_list, - double avg_res, struct cell_params *cparams) + struct image *images, int n, + double avg_res, UnitCell *avg) { int cp, i, u; int num_clens; - int max_clens; int best_clen; int *clens_population; double *clens; @@ -619,23 +446,18 @@ static double pick_clen_to_use(struct geoptimiser_params *gparams, double min_braggp_dist; double clen_to_use; struct rvec cqu; - - max_clens = 1024; - - clens = calloc(max_clens,sizeof(double)); - if ( clens == NULL ) { - ERROR("Failed to allocate memory for clen calculation.\n"); - return -1.0; - } - clens_population = calloc(max_clens,sizeof(int)); - if ( clens_population == NULL ) { - ERROR("Failed to allocate memory for clen calculation.\n"); - free(clens); - return -1.0; - } - lambdas = calloc(max_clens,sizeof(double)); - if ( lambdas == NULL ) { + double a, b, c, al, be, ga; + + /* These need to be big enough for the number of different camera + * lengths in the data set. There are probably only a few, but assume + * the worst case here - a unique camera length for each frame */ + clens = calloc(n, sizeof(double)); + clens_population = calloc(n, sizeof(int)); + lambdas = calloc(n, sizeof(double)); + if ((lambdas == NULL) || (clens == NULL) || (clens_population == NULL)) + { ERROR("Failed to allocate memory for clen calculation.\n"); + free(lambdas); free(clens); free(clens_population); return -1.0; @@ -643,60 +465,29 @@ static double pick_clen_to_use(struct geoptimiser_params *gparams, num_clens = 0; - for ( cp=0; cp<pattern_list->n_patterns; cp++ ) { + for ( cp=0; cp<n; cp++ ) { int i; int found = 0; for ( i=0; i<num_clens; i++ ) { - if ( fabs(pattern_list->patterns[cp]->clen-clens[i]) - <0.0001 ) { + if ( fabs(images[cp].avg_clen - clens[i]) <0.0001 ) { clens_population[i]++; - lambdas[i] += pattern_list->patterns[cp]->lambda; + lambdas[i] += images[cp].lambda; found = 1; break; } } - if ( found == 1) continue; - - if ( num_clens == max_clens ) { - - int *clens_population_new; - double *clens_new; - double *lambdas_new; - - clens_population_new = realloc(clens_population, - (max_clens+1024)*sizeof(int)); - clens_new = realloc(clens_population, - (max_clens+1024)*sizeof(double)); - lambdas_new = realloc(clens_population, - (max_clens+1024)*sizeof(double)); + if ( found ) continue; - 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; - } - - max_clens += 1024; - clens_population_new = clens_population; - clens = clens_new; - lambdas = lambdas_new; - } - - clens[num_clens] = pattern_list->patterns[cp]->clen; + clens[num_clens] = images[cp].avg_clen; clens_population[num_clens] = 1; - lambdas[num_clens] = pattern_list->patterns[cp]->lambda; + lambdas[num_clens] = images[cp].lambda; num_clens++; } - for ( u=0; u<num_clens; u++ ) { lambdas[u] /= clens_population[u]; } @@ -711,35 +502,38 @@ static double pick_clen_to_use(struct geoptimiser_params *gparams, best_clen = 0; clen_to_use = clens[0]; - for ( i=0; i<num_clens; i++) { - if ( clens_population[i] >0) { - 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], - 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]; - } + cell_get_parameters(avg, &a, &b, &c, &al, &be, &ga); + for ( i=0; i<num_clens; i++ ) { + + assert(clens_population[i] > 0); + + cqu = get_q_from_xyz(1.0/avg_res, 0, clens[i], lambdas[i]); + + min_braggp_dist = fmin(fmin((1.0/cqu.u)/a, + (1.0/cqu.u)/b), + (1.0/cqu.u)/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], 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]; } + } 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); + STATUS("Only %i patterns with camera length %0.4f m will be " + "used.\n", clens_population[best_clen], clen_to_use); } free(clens); @@ -763,8 +557,8 @@ static int find_quad_for_connected(struct rigid_group *rg, struct panel *p; int qi; - // The quadrant for a group of connected panels is the quadrant to which - // the first panel in the connected set belong + /* The quadrant for a group of connected panels is the quadrant to which + * the first panel in the connected set belongs */ p = rg->panels[0]; for ( qi=0; qi<quadrants->n_rigid_groups; qi++ ) { @@ -773,51 +567,46 @@ static int find_quad_for_connected(struct rigid_group *rg, } } - // Hopefully never reached - return -1; + /* Hopefully never reached */ + ERROR("Couldn't find quadrant for connected group!\n"); + abort(); } -static int fill_avg_pixel_displ(struct pixel_displ_list *pix_displ_list, - int pix_index, - struct avg_displacements *avg_displ) +/* Take all the (valid) displacements for pixel "i" in panel "gp", calculate + * the median displacements in each direction and the modulus */ +static int fill_avg_pixel_displ(struct gpanel *gp, int i) { 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 ) { + list_fs_displ = calloc(gp->num_pix_displ[i], sizeof(double)); + list_ss_displ = calloc(gp->num_pix_displ[i], sizeof(double)); + if ( (list_fs_displ == NULL) || (list_ss_displ == NULL) ) { ERROR("Failed to allocate memory for pixel statistics.\n"); free(list_fs_displ); + free(list_ss_displ); return 1; } - pix_displ_list->curr_pix_displ[pix_index] = - &pix_displ_list->pix_displ_list[pix_index]; + gp->curr_pix_displ[i] = &gp->pix_displ_list[i]; - for ( ei=0; ei<pix_displ_list->num_pix_displ[pix_index]; ei++ ) { + for ( ei=0; ei<gp->num_pix_displ[i]; ei++ ) { - struct single_pixel_displ * pixel_displ_entry; - pixel_displ_entry = pix_displ_list->curr_pix_displ[pix_index]; + struct single_pixel_displ *pix; - if ( pixel_displ_entry->dfs == -10000.0 ) break; - list_fs_displ[count] = pixel_displ_entry->dfs; - list_ss_displ[count] = pixel_displ_entry->dss; + pix = gp->curr_pix_displ[i]; + + if ( pix->dfs == -10000.0 ) break; + list_fs_displ[count] = pix->dfs; + list_ss_displ[count] = pix->dss; count++; - if ( pixel_displ_entry->ne == NULL ) { + if ( pix->ne == NULL ) { break; } else { - pix_displ_list->curr_pix_displ[pix_index] = - pix_displ_list->curr_pix_displ[pix_index]->ne; + gp->curr_pix_displ[i] = gp->curr_pix_displ[i]->ne; } } @@ -827,19 +616,20 @@ static int fill_avg_pixel_displ(struct pixel_displ_list *pix_displ_list, 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]); + /* FIXME: Is this fs/ss or x/y ??? */ + gp->avg_displ_x[i] = comp_median(list_fs_displ, count); + gp->avg_displ_y[i] = comp_median(list_ss_displ, count); + gp->avg_displ_abs[i] = modulus2d(gp->avg_displ_x[i], + gp->avg_displ_y[i]); + free(list_fs_displ); free(list_ss_displ); - return 0; } -static int allocate_next_element(struct single_pixel_displ** curr_pix_displ, - int pix_index) +static int allocate_next_element(struct single_pixel_displ **curr_pix_displ, + int pix_index) { curr_pix_displ[pix_index]->ne = malloc(sizeof(struct single_pixel_displ)); if ( curr_pix_displ[pix_index]->ne == NULL ) { @@ -853,167 +643,153 @@ static int allocate_next_element(struct single_pixel_displ** curr_pix_displ, } -static int add_distance_to_list(struct enhanced_det *edet, +static int add_distance_to_list(struct gpanel *gp, 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)); + pix_index = ((int)rint(imfe->fs) + gp->p->w*(int)rint(imfe->ss)); - if ( pix_displ_list->num_pix_displ[pix_index]>0 ) { + if ( gp->num_pix_displ[pix_index] > 0 ) { int ret; - ret = allocate_next_element(pix_displ_list->curr_pix_displ, - pix_index); + ret = allocate_next_element(gp->curr_pix_displ, pix_index); - if ( ret != 0) return ret; + 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]++; + compute_x_y(rfs, rss, get_panel(refl), &crx, &cry); + gp->curr_pix_displ[pix_index]->dfs = fx - crx; + gp->curr_pix_displ[pix_index]->dss = fy - cry; + gp->curr_pix_displ[pix_index]->ne = NULL; + gp->num_pix_displ[pix_index]++; return 0; } -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) +static int count_pixels_with_min_peaks(struct gpanel *gp, int min_num_peaks) { - int pix_index; int pixel_count = 0; int ifs, iss; - for ( ifs=p->min_fs; ifs<p->max_fs+1; ifs++ ) { - for ( iss=p->min_ss; iss<p->max_ss+1; iss++ ) { - - pix_index = ifs+edet->width*iss; - - if ( pix_displ_list->num_pix_displ[pix_index] >= - min_num_peaks ) { - pixel_count += 1; - } + for ( iss=0; iss<gp->p->h; iss++ ) { + for ( ifs=0; ifs<gp->p->w; ifs++ ) { + int idx = ifs+gp->p->w*iss; + if ( gp->num_pix_displ[idx] >= min_num_peaks ) { + pixel_count += 1; } + + } } return pixel_count; } -static void adjust_min_peaks_per_conn(struct enhanced_det *edet, - struct rg_collection *connected, +static void adjust_min_peaks_per_conn(struct rg_collection *connected, + struct gpanel *gpanels, + struct detector *det, struct geoptimiser_params *gparams, - struct connected_data *conn_data, - struct pixel_displ_list *pix_displ_list) + struct connected_data *conn_data) { - int min_num_peaks, di, ip; - struct panel *p; - STATUS("Adjusting the minimum number of measurements per pixel in " "order to have enough measurements for each connected group.\n"); + for ( di=0; di<connected->n_rigid_groups; di++ ) { for ( min_num_peaks=gparams->min_num_peaks_per_pix; - min_num_peaks>0; min_num_peaks-- ) { - + min_num_peaks>0; min_num_peaks-- ) + { int di_count = 0; - for (ip=0; ip<connected->rigid_groups[di]->n_panels; - ip++) { - + for ( ip=0; ip<connected->rigid_groups[di]->n_panels; + ip++ ) + { int pix_count; + struct panel *p; + struct gpanel *gp; p = connected->rigid_groups[di]->panels[ip]; + gp = &gpanels[panel_number(det, p)]; - - pix_count = count_pixels_with_min_peaks(p, - min_num_peaks, - pix_displ_list, - edet); + pix_count = count_pixels_with_min_peaks(gp, + min_num_peaks); 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; + if ( di_count >= gparams->min_num_pix_per_conn_group ) { + conn_data[di].num_peaks_per_pixel = min_num_peaks; break; } } - STATUS("Minimum number of measurement " - "per pixel for connected group " - "%s has been set to %i\n", - conn_data[di].name, + STATUS("Minimum number of measurements per pixel for connected " + "group %s has been set to %i\n", conn_data[di].name, conn_data[di].num_peaks_per_pixel); } } -static int compute_pixel_displacements(struct pattern_list *pattern_list, - struct enhanced_det *edet, +static int compute_pixel_displacements(struct image *images, int n_images, + struct gpanel *gpanels, + struct detector *det, struct rg_collection *connected, struct geoptimiser_params *gparams, double clen_to_use, - struct connected_data *conn_data, - struct avg_displacements *avg_displ, - struct pixel_displ_list *pix_displ_list) + struct connected_data *conn_data) { - int cp, ich; - + int cp; STATUS("Computing pixel displacements.\n"); - for ( cp=0; cp<pattern_list->n_patterns; cp++ ) { - ImageFeatureList *flist = pattern_list->patterns[cp]->im_list; + for ( cp=0; cp<n_images; cp++ ) { + + int fi; + ImageFeatureList *flist = images[cp].features; if ( gparams->only_best_distance ) { - if ( fabs(pattern_list->patterns[cp]->clen-clen_to_use) > - 0.0001 ) { + if ( fabs(images[cp].avg_clen - clen_to_use) > 0.0001 ) { continue; } } - for ( ich=0; - ich<image_feature_count(pattern_list->patterns[cp]->im_list); - ich++ ) { + for ( fi=0; fi<image_feature_count(images[cp].features); fi++ ) { double min_dist; double fx, fy; Reflection *refl; - int ret; + struct imagefeature *imfe; - RefList *rlist = pattern_list->patterns[cp]->ref_list; + imfe = image_get_feature(flist, fi); + if ( imfe == NULL ) continue; - struct imagefeature *imfe = image_get_feature(flist, ich); - compute_x_y(edet->det, imfe->fs, imfe->ss, &fx, &fy); + compute_x_y(imfe->fs, imfe->ss, imfe->p, &fx, &fy); - refl = find_closest_reflection(rlist, fx, fy, edet->det, + /* Find the closest reflection (from all crystals) */ + refl = find_closest_reflection(&images[cp], fx, fy, &min_dist); - if ( refl == NULL ) continue; if ( min_dist < gparams->max_peak_dist ) { + struct gpanel *gp; + int r; + gp = &gpanels[panel_number(det, imfe->p)]; - ret = add_distance_to_list(edet, imfe, - pix_displ_list, - refl, fx, fy); - - if ( ret != 0 ) return ret; + r = add_distance_to_list(gp, imfe, refl, fx, fy); + if ( r ) return r; } } @@ -1023,23 +799,22 @@ static int compute_pixel_displacements(struct pattern_list *pattern_list, } -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) +static int compute_avg_pix_displ(struct gpanel *gp, int idx, + int num_peaks_per_pixel) { int ret; - if ( pix_displ_list->num_pix_displ[pix_index] >= - num_peaks_per_pixel ) { + if ( gp->num_pix_displ[idx] >= num_peaks_per_pixel ) { + + ret = fill_avg_pixel_displ(gp, idx); + if ( ret != 0 ) return ret; - 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; + + gp->avg_displ_x[idx] = -10000.0; + gp->avg_displ_y[idx] = -10000.0; + gp->avg_displ_abs[idx] = -10000.0; + } return 0; @@ -1047,11 +822,10 @@ static int compute_avg_pix_displ(struct pixel_displ_list *pix_displ_list, } -static int compute_avg_displacements(struct enhanced_det *edet, +static int compute_avg_displacements(struct detector *det, struct rg_collection *connected, - struct pixel_displ_list *pix_displ_list, struct connected_data *conn_data, - struct avg_displacements *avg_displ) + struct gpanel *gpanels) { int di, ip, ifs, iss; int pix_index, ret; @@ -1059,23 +833,26 @@ static int compute_avg_displacements(struct enhanced_det *edet, for ( di=0; di<connected->n_rigid_groups; di++ ) { - for (ip=0; ip<connected->rigid_groups[di]->n_panels; - ip++) { + for ( ip=0; ip<connected->rigid_groups[di]->n_panels; ip++ ) { + + int pp; + struct gpanel *gp; p = connected->rigid_groups[di]->panels[ip]; + pp = panel_number(det, p); + gp = &gpanels[pp]; - for ( ifs=p->min_fs; ifs<p->max_fs+1; ifs++ ) { - for ( iss=p->min_ss; iss<p->max_ss+1; iss++ ) { + for ( iss=0; iss<p->h; iss++ ) { + for ( ifs=0; ifs<p->w; ifs++ ) { - pix_index = ifs+edet->width*iss; + pix_index = ifs+p->w*iss; - ret = compute_avg_pix_displ(pix_displ_list, - conn_data[di].num_peaks_per_pixel, - pix_index, avg_displ); + ret = compute_avg_pix_displ(gp, pix_index, + conn_data[di].num_peaks_per_pixel); - if ( ret != 0 ) return ret; + if ( ret != 0 ) return ret; - } + } } } @@ -1085,10 +862,9 @@ static int compute_avg_displacements(struct enhanced_det *edet, static double compute_error(struct rg_collection *connected, - struct enhanced_det* edet, + struct detector *det, struct connected_data *conn_data, - int *num_pix_displ, - double *displ_abs) + struct gpanel *gpanels) { double total_error = 0; int num_total_error = 0; @@ -1101,30 +877,35 @@ static double compute_error(struct rg_collection *connected, int num_connected_error = 0; int ifs, iss; - for (ip=0; ip<connected->rigid_groups[di]->n_panels; ip++) { + for ( ip=0; ip<connected->rigid_groups[di]->n_panels; ip++ ) { p = connected->rigid_groups[di]->panels[ip]; - for (ifs=p->min_fs; ifs<p->max_fs+1; ifs++) { - for (iss=p->min_ss; iss<p->max_ss+1; iss++) { - - int pix_index; - pix_index = ifs+edet->width*iss; + for ( iss=0; iss<p->h; iss++ ) { + for ( ifs=0; ifs<p->w; ifs++ ) { - if ( num_pix_displ[pix_index]>= - conn_data[di].num_peaks_per_pixel ) { + int pix_index; + struct gpanel *gp; + int pp = panel_number(det, p); - double cer; + gp = &gpanels[pp]; + pix_index = ifs+p->w*iss; - cer = displ_abs[pix_index]* - displ_abs[pix_index]; - connected_error += cer; - num_connected_error++; - total_error += cer; - num_total_error++; - } + if ( gp->num_pix_displ[pix_index] + >= conn_data[di].num_peaks_per_pixel ) + { + double cer; + + cer = gp->avg_displ_abs[pix_index] + * gp->avg_displ_abs[pix_index]; + connected_error += cer; + num_connected_error++; + total_error += cer; + num_total_error++; } } + } + } if ( num_connected_error > 0 ) { @@ -1151,64 +932,7 @@ static double compute_error(struct rg_collection *connected, } -static struct pixel_maps *initialize_pixel_maps(struct enhanced_det *edet) -{ - - int pi; - struct pixel_maps *pixel_maps; - - 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; pi<edet->det->n_panels; pi++ ) { - - struct panel *p; - int iss, ifs; - - 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++) { - - double xs, ys; - int pix_index; - - 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; -} - - -void free_pixel_maps(struct pixel_maps* pixel_maps) -{ - free(pixel_maps->pix_to_x); - free(pixel_maps->pix_to_y); - free(pixel_maps); -} - - - -struct avg_rot_and_stretch* initialize_avg_rot_stretch( - int num_rigid_groups) +struct avg_rot_and_stretch *initialize_avg_rot_stretch(int num_rigid_groups) { int i; @@ -1433,67 +1157,54 @@ static void correct_rotation_and_stretch(struct rg_collection *connected, } -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) +static void adjust_panel(struct connected_data *conn_data, + struct rg_collection *connected, double c_stretch, + double stretch_coeff, int num_peaks_per_pixel, + struct panel *p, struct gpanel *gp) { - double c_stretch; - struct panel *p; int ifs, iss; - c_stretch = conn_data[di].cstr; - + /* FIXME: What does "TODO" mean? */ //TODO if ( fabs(c_stretch)<FLT_EPSILON ) c_stretch = stretch_coeff; - p = connected->rigid_groups[di]->panels[ip]; + for ( iss=0; iss<p->h; iss++ ) { + for ( ifs=0; ifs<p->w; ifs++ ) { - for ( ifs=p->min_fs; ifs<p->max_fs+1; ifs++ ) { - for ( iss=p->min_ss; iss<p->max_ss+1; iss++ ) { + double x, y; + int idx = ifs+gp->p->w*iss; - int pix_index = ifs+edet->width*iss; + if ( gp->num_pix_displ[idx] < num_peaks_per_pixel) continue; - 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) { + compute_x_y(ifs, iss, p, &x, &y); + gp->avg_displ_x[idx] -= x - x/c_stretch; + gp->avg_displ_y[idx] -= y - y/c_stretch; - 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 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) + double stretch_coeff, struct detector *det, + struct gpanel *gpanels) { - int di, ip; for ( di=0; di<connected->n_rigid_groups; di++ ) { for (ip=0; ip<connected->rigid_groups[di]->n_panels; ip++) { - adjust_panel(conn_data, di, ip, connected, - pixel_maps, recomputed_pixel_maps, - avg_displ, stretch_coeff, edet, - num_pix_displ); + struct panel *p; + struct gpanel *gp; + + p = connected->rigid_groups[di]->panels[ip]; + gp = &gpanels[panel_number(det, p)]; + + adjust_panel(conn_data, connected, conn_data[di].cstr, + stretch_coeff, + conn_data[di].num_peaks_per_pixel, p, gp); } } } @@ -1501,103 +1212,101 @@ static void adjust_displ_for_stretch(struct rg_collection *connected, static void fill_av_conn(struct rg_collection *connected, int di, struct connected_data *conn_data, - int *num_pix_displ, struct enhanced_det *edet, + struct detector *det, struct gpanel *gpanels, double *list_displ_in_conn_fs, - double *list_displ_in_conn_ss, - struct avg_displacements *avg_displ) + double *list_displ_in_conn_ss) { - struct panel *p; - int ifs, iss, ip; - + int ip; int counter = 0; - for (ip=0; ip<connected->rigid_groups[di]->n_panels; ip++) { - - p = connected->rigid_groups[di]->panels[ip]; - - for ( ifs=p->min_fs; ifs<p->max_fs+1; ifs++ ) { - for ( iss=p->min_ss; iss<p->max_ss+1; iss++ ) { + for ( ip=0; ip<connected->rigid_groups[di]->n_panels; ip++ ) { - 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++; - } + int ifs, iss; + struct panel *p = connected->rigid_groups[di]->panels[ip]; + struct gpanel *gp = &gpanels[panel_number(det, p)]; + + for ( iss=0; iss<p->h; iss++ ) { + for ( ifs=0; ifs<p->w; ifs++ ) { + + int pix_index = ifs+p->w*iss; + + if ( gp->num_pix_displ[pix_index] + >= conn_data[di].num_peaks_per_pixel ) + { + list_displ_in_conn_fs[counter] = + gp->avg_displ_x[pix_index]; + list_displ_in_conn_ss[counter] = + gp->avg_displ_y[pix_index]; + counter++; } + + } } } if ( counter != conn_data[di].n_peaks_in_conn ) { - printf("counter: %i n_peaks_in_conn: %i\n", + ERROR("counter: %i n_peaks_in_conn: %i\n", counter, conn_data[di].n_peaks_in_conn); - exit(0); + abort(); } } 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) + double *av_in_panel_fs, double *av_in_panel_ss, + int di, double mpd) { conn_data[di].sh_x = comp_median(av_in_panel_fs, 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 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 ) { - 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), - max_peak_distance); + + if ( modulus2d(conn_data[di].sh_x, conn_data[di].sh_y ) > 0.8*mpd ) { + STATUS("WARNING: absolute shift is: %0.1f > 0.8*%0.1f pixels. " + "Increase the value of max_peak_distance!\n", + modulus2d(conn_data[di].sh_x, conn_data[di].sh_y), mpd); } } static int compute_shift(struct rg_collection *connected, struct connected_data *conn_data, - int *num_pix_displ, - struct enhanced_det *edet, + struct detector *det, struct geoptimiser_params *gparams, - struct avg_displacements *avg_displ) + struct gpanel *gpanels) { - STATUS("Computing shift corrections.\n"); - int di; + STATUS("Computing shift corrections.\n"); + for ( di=0; di<connected->n_rigid_groups; di++ ) { double *list_displ_in_conn_fs; double *list_displ_in_conn_ss; - 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; - } - list_displ_in_conn_ss = malloc(conn_data[di].n_peaks_in_conn* - sizeof(double)); - if ( list_displ_in_conn_ss == NULL ) { + list_displ_in_conn_fs = malloc(conn_data[di].n_peaks_in_conn + * sizeof(double)); + list_displ_in_conn_ss = malloc(conn_data[di].n_peaks_in_conn + * sizeof(double)); + if ( (list_displ_in_conn_fs == NULL) + || (list_displ_in_conn_ss == NULL) ) + { ERROR("Failed to allocate memory for computing shifts\n"); free(list_displ_in_conn_fs); + free(list_displ_in_conn_ss); return 1; } - fill_av_conn(connected, di, conn_data, num_pix_displ, edet, - list_displ_in_conn_fs, list_displ_in_conn_ss, - avg_displ); + fill_av_conn(connected, di, conn_data, det, gpanels, + list_displ_in_conn_fs, list_displ_in_conn_ss); - if ( conn_data[di].n_peaks_in_conn >= - gparams->min_num_pix_per_conn_group ) { + 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, @@ -1610,6 +1319,7 @@ static int compute_shift(struct rg_collection *connected, free(list_displ_in_conn_fs); free(list_displ_in_conn_ss); } + return 0; } @@ -1794,11 +1504,8 @@ static struct connected_stretch_and_angles *initialize_connected_stretch_angles( 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, + struct detector *det, struct gpanel *gpanels, int di, double min_dist, long *num_ang, int ifs0, int iss0, double c_x0, double c_y0, double cd_x0, double cd_y0, @@ -1806,8 +1513,11 @@ static void scan_p1(int ip0, int ip1, struct rg_collection *connected, { int iss1, ifs1; + struct panel *p1; + struct gpanel *gp1; - struct panel *p1 = connected->rigid_groups[di]->panels[ip1]; + p1 = connected->rigid_groups[di]->panels[ip1]; + gp1 = &gpanels[panel_number(det, p1)]; int min_ss_p1, min_fs_p1; @@ -1815,129 +1525,115 @@ static void scan_p1(int ip0, int ip1, struct rg_collection *connected, min_fs_p1 = ifs0; min_ss_p1 = iss0; } else { - min_fs_p1 = p1->min_fs; - min_ss_p1 = p1->min_ss; + min_fs_p1 = 0; + min_ss_p1 = 0; } - for ( iss1=min_ss_p1; iss1<p1->max_ss+1; iss1++ ) { - - for ( ifs1=min_fs_p1; ifs1<p1->max_fs+1; ifs1++ ) { - - int pix_index1 = ifs1+edet->width*iss1; + for ( iss1=min_ss_p1; iss1<p1->h; iss1++ ) { + for ( ifs1=min_fs_p1; ifs1<p1->w; ifs1++ ) { - 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; + int pix_index1 = ifs1+p1->w*iss1; - 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; + if ( gp1->num_pix_displ[pix_index1] + < conn_data[di].num_peaks_per_pixel ) continue; - 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<FLT_EPSILON || len2<FLT_EPSILON ) { - continue; - } + compute_x_y(ifs1, iss1, p1, &c_x1, &c_y1); + cd_x1 = c_x1 - gp1->avg_displ_x[pix_index1]; + cd_y1 = c_y1 - gp1->avg_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; - if (compute) { + dist = modulus2d(d_c_x,d_c_y); + if ( dist < min_dist ) continue; - double scal_m; - double multlen; + len1 = modulus2d(d_c_x, d_c_y); + len2 = modulus2d(d_cd_x, d_cd_y); + if ( len1<FLT_EPSILON || len2<FLT_EPSILON ) continue; - scal_m = d_c_x * d_cd_x+ - d_c_y * d_cd_y - - FLT_EPSILON; + if ( compute ) { - multlen = len1*len2; - if ( fabs(scal_m)>=multlen ) { - angles[*num_ang] = 0.0; - } else { + double scal_m; + double multlen; - angles[*num_ang] = acos(scal_m - /multlen); + scal_m = d_c_x * d_cd_x+ + d_c_y * d_cd_y - + FLT_EPSILON; - if (d_c_y * d_cd_x - - d_c_x * d_cd_y < 0) { - angles[*num_ang] *= -1.; - } + multlen = len1*len2; + if ( fabs(scal_m)>=multlen ) { + angles[*num_ang] = 0.0; + } else { - } + angles[*num_ang] = acos(scal_m/multlen); - stretches[*num_ang] = len1/len2; + if (d_c_y * d_cd_x - d_c_x * d_cd_y < 0) { + angles[*num_ang] *= -1.; } - *num_ang = *num_ang+1; } + + stretches[*num_ang] = len1/len2; } + + *num_ang = *num_ang+1; + } } } +/* Executed for each panel in the connected group */ 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, + struct detector *det, struct gpanel *gpanels, int di, double min_dist, long *num_ang, int compute, double *angles, double *stretches) { - int iss0, ifs0, ip1; + struct gpanel *gp; + struct panel *p0; - struct panel *p0 = connected->rigid_groups[di]->panels[ip0]; + p0 = connected->rigid_groups[di]->panels[ip0]; + gp = &gpanels[panel_number(det, p0)]; - for ( iss0=p0->min_ss; iss0<p0->max_ss+1; iss0++ ) { + for ( iss0=0; iss0<p0->h; iss0++ ) { + for ( ifs0=0; ifs0<p0->w; ifs0++ ) { - for ( ifs0=p0->min_fs; ifs0<p0->max_fs+1; ifs0++ ) { + double c_x0, c_y0, cd_x0, cd_y0; + int pix_index0 = ifs0+p0->w*iss0; - int pix_index0 = ifs0+edet->width*iss0; + if ( gp->num_pix_displ[pix_index0] + < conn_data[di].num_peaks_per_pixel ) continue; - if ( num_pix_displ[pix_index0]>= - conn_data[di].num_peaks_per_pixel ) { + compute_x_y(ifs0, iss0, p0, &c_x0, &c_y0); + cd_x0 = c_x0 - gp->avg_displ_x[pix_index0]; + cd_y0 = c_y0 - gp->avg_displ_y[pix_index0]; - 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; - ip1<connected->rigid_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); - } - } + for ( ip1=ip0; ip1<connected->rigid_groups[di]->n_panels; + ip1++ ) + { + scan_p1(ip0, ip1, connected, + conn_data, det, gpanels, di, min_dist, + num_ang, ifs0, iss0, c_x0, + c_y0, cd_x0, cd_y0, compute, + angles, stretches); } + + } } } 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, + struct detector *det, + struct gpanel *gpanels, double dist_coeff_for_rot_str, struct geoptimiser_params *gparams) { @@ -1949,15 +1645,9 @@ static double compute_rotation_and_stretch(struct rg_collection *connected, STATUS("Computing rotation and stretch corrections.\n"); csaa = initialize_connected_stretch_angles(connected); - if ( csaa == NULL ) { - return -1.0; - } + if ( csaa == NULL ) return -1.0; for ( di=0; di<connected->n_rigid_groups; di++ ) { - if ( conn_data[di].n_peaks_in_conn < - gparams->min_num_pix_per_conn_group ) { - continue; - } long max_num_ang = 0; @@ -1970,27 +1660,28 @@ static double compute_rotation_and_stretch(struct rg_collection *connected, int ip0; int num_pix_first_p; + /* Enough peaks in this connected group? */ + if ( conn_data[di].n_peaks_in_conn + < gparams->min_num_pix_per_conn_group ) continue; + first_p = connected->rigid_groups[di]->panels[0]; - num_pix_first_p = (first_p->max_fs+1-first_p->min_fs)* - (first_p->max_ss+1-first_p->min_ss); + num_pix_first_p = first_p->w * first_p->h; - // TODO: MINRAD HERE IS NOT UNIVERSAL - min_dist = dist_coeff_for_rot_str*sqrt(num_pix_first_p* - connected->rigid_groups[di]->n_panels); + /* FIXME: minrad here is not universal */ + min_dist = dist_coeff_for_rot_str * sqrt(num_pix_first_p + * connected->rigid_groups[di]->n_panels); for ( ip0=0; ip0<connected->rigid_groups[di]->n_panels; ip0++ ) { - - scan_p0(ip0, connected, avg_displ, conn_data, edet, - pixel_maps, num_pix_displ, di, min_dist, - &num_ang, 0, NULL, NULL); + scan_p0(ip0, connected, conn_data, det, gpanels, + di, min_dist, &num_ang, 0, NULL, NULL); } max_num_ang = num_ang+1; angles = malloc(max_num_ang*sizeof(double)); if ( angles == NULL ) { - ERROR("Error in allocating memory for angle " + ERROR("Error allocating memory for angle " "optimization\n"); free(csaa->stretch_coeff); free(csaa->num_angles); @@ -1999,7 +1690,7 @@ static double compute_rotation_and_stretch(struct rg_collection *connected, } stretches = malloc(max_num_ang*sizeof(double)); if ( stretches == NULL ) { - ERROR("Error in allocating memory for stretch " + ERROR("Error allocating memory for stretch " "optimization\n"); free(angles); free(csaa->stretch_coeff); @@ -2011,19 +1702,18 @@ static double compute_rotation_and_stretch(struct rg_collection *connected, num_ang = 0; for ( ip0=0; ip0<connected->rigid_groups[di]->n_panels; ip0++ ) { - - scan_p0(ip0, connected, avg_displ, conn_data, edet, - pixel_maps, num_pix_displ, di, min_dist, - &num_ang, 1, angles, stretches); + scan_p0(ip0, connected, conn_data, det, gpanels, + 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); 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); + "%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; @@ -2178,8 +1868,8 @@ static int draw_detector(cairo_surface_t *surf, struct image *image, } -static int save_data_to_png(char *filename, struct enhanced_det *edet, - double *data) +static int save_data_to_png(char *filename, struct detector *det, + struct gpanel *gpanels) { struct image im; int i; @@ -2190,21 +1880,19 @@ static int save_data_to_png(char *filename, struct enhanced_det *edet, cairo_status_t r; cairo_surface_t *surf; - im.det = edet->det; - im.width = edet->width; - im.height = edet->height; + im.det = det; im.bad = NULL; - im.dp = malloc(edet->det->n_panels*sizeof(float *)); + im.dp = malloc(det->n_panels*sizeof(float *)); if ( im.dp == NULL ) { ERROR("Failed to allocate data\n"); return 1; } - for ( i=0; i<edet->det->n_panels; i++ ) { + for ( i=0; i<det->n_panels; i++ ) { int fs, ss; struct panel *p; - p = &edet->det->panels[i]; + p = &det->panels[i]; im.dp[i] = calloc(p->w * p->h, sizeof(float)); if ( im.dp[i] == NULL ) { @@ -2216,19 +1904,16 @@ static int save_data_to_png(char *filename, struct enhanced_det *edet, for ( fs=0; fs<p->w; fs++ ) { int idx; - int cfs, css; float val; - cfs = fs+p->min_fs; - css = ss+p->min_ss; - idx = cfs + css*edet->width; + idx = fs + ss*p->w; - if ( data[idx] == -10000.0) { + if ( gpanels[i].avg_displ_abs[idx] == -10000.0) { val = 0.0; - } else if ( data[idx] > 1.0) { + } else if ( gpanels[i].avg_displ_abs[idx] > 1.0) { val = 1.0; } else { - val = (float)data[idx]; + val = (float)gpanels[i].avg_displ_abs[idx]; } val *= 10.0; /* render_panels sets this as max */ @@ -2267,7 +1952,7 @@ static int save_data_to_png(char *filename, struct enhanced_det *edet, cairo_fill(cr); cairo_destroy(cr); - for ( i=0; i<edet->det->n_panels; i++ ) { + for ( i=0; i<det->n_panels; i++ ) { free(im.dp[i]); } free(im.dp); @@ -2387,224 +2072,128 @@ static struct connected_data *initialize_conn_data(struct rg_collection *quadran } -static struct pixel_displ_list *initialize_pixel_displacement_list( - struct enhanced_det *edet) +static int initialize_pixel_displacement_list(struct gpanel *gp) { - - 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 ) { + gp->pix_displ_list = calloc(gp->p->w*gp->p->h, + sizeof(struct single_pixel_displ)); + if ( gp->pix_displ_list == NULL ) { ERROR("Error allocating memory for pixel displacement data.\n"); - free(pix_displ_list); - return NULL; + return 1; } - pix_displ_list->curr_pix_displ = calloc(edet->num_pix, - sizeof(struct single_pixel_displ*)); - if ( pix_displ_list->curr_pix_displ == NULL ) { + + gp->curr_pix_displ = calloc(gp->p->w*gp->p->h, + sizeof(struct single_pixel_displ *)); + if ( gp->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; + free(gp->pix_displ_list); + return 1; } - pix_displ_list->num_pix_displ = calloc(edet->num_pix, sizeof(int)); - if ( pix_displ_list->num_pix_displ == NULL ) { + gp->num_pix_displ = calloc(gp->p->w*gp->p->h, sizeof(int)); + if ( gp->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; + free(gp->pix_displ_list); + free(gp->curr_pix_displ); + return 1; } - for ( ipx=0; ipx<edet->num_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; + for ( ipx=0; ipx<gp->p->w*gp->p->h; ipx++ ) { + gp->pix_displ_list[ipx].dfs = -10000.0; + gp->pix_displ_list[ipx].dss = -10000.0; + gp->pix_displ_list[ipx].ne = NULL; + gp->curr_pix_displ[ipx] = &gp->pix_displ_list[ipx]; + gp->num_pix_displ[ipx] = 0; } - return pix_displ_list; + return 0; } -static void free_pixel_displacement_list( - struct pixel_displ_list *pix_displ_list, - struct enhanced_det *edet) +static void free_displ_lists(struct gpanel *gpanels, int n) { - int i; + int j; struct single_pixel_displ *curr = NULL; struct single_pixel_displ *next = NULL; - for ( i=0; i<edet->num_pix; i++ ) { - - curr = &pix_displ_list->pix_displ_list[i]; + for ( j=0; j<n; j++ ) { - 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; pti<pattern_list->n_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; nuc<pattern_list->patterns[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; + int i; + struct gpanel *gp = &gpanels[j]; - for ( i=0; i<edet->num_pix; i++ ) { + for ( i=0; i<gp->p->w*gp->p->h; i++ ) { - curr = &pix_displ_list->pix_displ_list[i]; + curr = &gp->pix_displ_list[i]; - if ( curr->ne != NULL ) { - curr = curr->ne; - while ( curr != NULL ) { - next = curr->ne; - free(curr); - curr = next; + 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(gp->curr_pix_displ); + free(gp->pix_displ_list); - 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, struct gpanel *gp, + int num_peaks_per_pixel, + double sh_x, double sh_y) { - struct panel *p; int ifs, iss; - if (conn_data[di].sh_x < -9999.0) return; + for ( iss=0; iss<p->h; iss++ ) { + for ( ifs=0; ifs<p->w; ifs++ ) { - p = connected->rigid_groups[di]->panels[ip]; + int pix_index = ifs+p->w*iss; - for ( ifs=p->min_fs; ifs<p->max_fs+1; ifs++ ) { - for ( iss=p->min_ss; iss<p->max_ss+1; iss++ ) { + if ( gp->num_pix_displ[pix_index] >= num_peaks_per_pixel ) { - int pix_index = ifs+edet->width*iss; + gp->avg_displ_x[pix_index] -= sh_x; + gp->avg_displ_y[pix_index] -= sh_y; + gp->avg_displ_abs[pix_index] = modulus2d( + gp->avg_displ_x[pix_index], + gp->avg_displ_y[pix_index]); - 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; - } + } else { + gp->avg_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) + struct detector *det, struct gpanel *gpanels) { int di, ip; for ( di=0;di<connected->n_rigid_groups;di++ ) { + + if (conn_data[di].sh_x < -9999.0) continue; + for (ip=0; ip<connected->rigid_groups[di]->n_panels; ip++) { - recompute_panel_avg_displ(connected, conn_data, - num_pix_displ, di, ip, edet, - avg_displ); + struct panel *p; + struct gpanel *gp; + + p = connected->rigid_groups[di]->panels[ip]; + gp = &gpanels[panel_number(det, p)]; + recompute_panel_avg_displ(connected, conn_data, p, gp, + conn_data[di].num_peaks_per_pixel, + conn_data[di].sh_x, + conn_data[di].sh_y); } } @@ -2616,13 +2205,9 @@ int optimize_geometry(struct geoptimiser_params *gparams, struct rg_collection *quadrants, struct rg_collection *connected) { - int max_fs = 0; - int max_ss = 0; int pi; int ret; int write_ret; - int maybe_cspad = 0; - int *num_pix_displ; double res_sum; double avg_res; @@ -2634,32 +2219,24 @@ int optimize_geometry(struct geoptimiser_params *gparams, double dist_coeff_for_rot_str = 0.2; double total_error; - struct pixel_maps *pixel_maps; - struct pixel_maps *recomputed_pixel_maps; double stretch_coeff = 1.0; 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; + struct image *images; + int n_images = 0; + UnitCell *avg_cell; + struct gpanel *gpanels; STATUS("Maximum distance between peaks: %0.1f pixels.\n", gparams->max_peak_dist); - 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", + 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 && !gparams->no_cspad ) { + if ( (det->n_panels == 64) && !gparams->no_cspad ) { int num_errors = 0; @@ -2668,10 +2245,8 @@ int optimize_geometry(struct geoptimiser_params *gparams, "connected ASICS.\n"); STATUS("If the detector is not a CSPAD, please rerun the " "program with the --no-cspad option.\n"); - if ( gparams->enforce_cspad_layout ) { - STATUS("Enforcing CSPAD layout...\n"); - } + STATUS("Enforcing CSPAD layout...\n"); num_errors = check_and_enforce_cspad_dist(gparams, det); if ( gparams->enforce_cspad_layout ) { @@ -2680,8 +2255,8 @@ int optimize_geometry(struct geoptimiser_params *gparams, STATUS("Saving geometry with enforced CSPAD layout.\n" "Please restart geometry optimization using the " - "optimized geometry from this run as input geometry " - "file.\n"); + "optimized geometry from this run as input " + "geometry file.\n"); geom_wr = write_detector_geometry_2( gparams->geometry_filename, gparams->outfile, det, @@ -2704,127 +2279,95 @@ int optimize_geometry(struct geoptimiser_params *gparams, } } - pattern_list = read_patterns_from_steam_file(gparams->infile, det); - if ( pattern_list->n_patterns < 1 ) { + images = read_patterns_from_stream(gparams->infile, det, &n_images); + if ( (n_images < 1) || (images == NULL) ) { ERROR("Error reading stream file\n"); return 1; } - cell_params = compute_avg_cell_parameters(pattern_list); - if ( cell_params == NULL ) { - free(pattern_list); + avg_cell = compute_avg_cell_parameters(images, n_images); + if ( avg_cell == NULL ) { + free(images); return 1; } res_sum = 0; for ( pi=0; pi<det->n_panels; pi++ ) { - - if ( det->panels[pi].max_fs > max_fs ) { - max_fs = det->panels[pi].max_fs; - } - if ( det->panels[pi].max_ss > max_ss ) { - max_ss = det->panels[pi].max_ss; - } res_sum += det->panels[pi].res; } avg_res = res_sum/det->n_panels; - 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 = pick_clen_to_use(gparams, pattern_list, avg_res, - cell_params); - - if ( clen_to_use == -1.0 ) return 1; + clen_to_use = pick_clen_to_use(gparams, images, n_images, avg_res, + avg_cell); + if ( clen_to_use < 0.0 ) return 1; - avg_displ = initialize_average_displacement(&edet); - if ( avg_displ == NULL ) { - free(cell_params); - free(pattern_list); + gpanels = calloc(det->n_panels, sizeof(struct gpanel)); + if ( gpanels == NULL ) { + ERROR("Failed to allocate panels.\n"); return 1; } + for ( pi=0; pi<det->n_panels; pi++ ) { - pixel_maps = initialize_pixel_maps(&edet); - if ( pixel_maps == NULL ) { - ERROR("Failed to allocate memory for pixel maps.\n"); - free_avg_displacement(avg_displ); - free(cell_params); - free(pattern_list); - return 1; - } + struct gpanel *gp = &gpanels[pi]; + int npx = det->panels[pi].w * det->panels[pi].h; + + gp->p = &det->panels[pi]; + + gp->avg_displ_x = malloc(npx*sizeof(double)); + gp->avg_displ_y = malloc(npx*sizeof(double)); + gp->avg_displ_abs = malloc(npx*sizeof(double)); + + if ( (gp->avg_displ_x == NULL) || (gp->avg_displ_y == NULL) + || (gp->avg_displ_abs == NULL) ) + { + ERROR("Failed to allocate displacements.\n"); + return 1; + } + + if ( initialize_pixel_displacement_list(gp) ) { + ERROR("Failed to allocate lists.\n"); + return 1; + } - pix_displ_list = initialize_pixel_displacement_list(&edet); - if ( pix_displ_list == NULL ) { - ERROR("Error allocating memory for connected structure data.\n"); - free_avg_displacement(avg_displ); - free_pixel_maps(pixel_maps); - free(pattern_list); - return 1; } conn_data = initialize_conn_data(quadrants, connected); if ( conn_data == NULL ) { ERROR("Error allocating memory for connected structure data.\n"); - free_avg_displacement(avg_displ); - free_pixel_maps(pixel_maps); - free_pixel_displacement_list(pix_displ_list, &edet); - free(pattern_list); - return 1; - } - - - ret = compute_pixel_displacements(pattern_list, &edet, connected, - gparams, clen_to_use, conn_data, - avg_displ, pix_displ_list); - if ( ret != 0 ) { - 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_pattern_list(pattern_list); + if ( compute_pixel_displacements(images, n_images, gpanels, det, + connected, gparams, clen_to_use, + conn_data) ) return 1; - adjust_min_peaks_per_conn(&edet, connected, gparams, conn_data, - pix_displ_list); + adjust_min_peaks_per_conn(connected, gpanels, det, gparams, conn_data); - 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); + if ( compute_avg_displacements(det, connected, conn_data, gpanels) ) { free(conn_data); + free(images); return 1; } - num_pix_displ = extract_num_pix_free_displ_list(pix_displ_list, &edet); + free_displ_lists(gpanels, det->n_panels); + /* gpanels[].num_pix_displ is still there */ STATUS("Computing error before correction.\n"); - total_error = compute_error(connected, &edet, conn_data, - num_pix_displ, avg_displ->displ_abs); + total_error = compute_error(connected, det, conn_data, gpanels); STATUS("Detector-wide error before correction: RMSD = %0.4f pixels.\n", total_error); if ( gparams->error_maps ) { + STATUS("Saving error map before correction.\n"); #ifdef HAVE_SAVE_TO_PNG - ret = save_data_to_png("error_map_before.png", &edet, - avg_displ->displ_abs); - if ( ret != 0 ) { + if ( save_data_to_png("error_map_before.png", det, gpanels) ) { ERROR("Error while writing data to file.\n"); - - free_avg_displacement(avg_displ); - free_pixel_maps(pixel_maps); free(conn_data); - free(num_pix_displ); + free(images); return 1; } @@ -2837,58 +2380,32 @@ int optimize_geometry(struct geoptimiser_params *gparams, } - stretch_coeff = compute_rotation_and_stretch(connected, conn_data, - &edet, num_pix_displ, - pixel_maps, avg_displ, + det, gpanels, 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(num_pix_displ); return 1; } 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); + if ( ret ) { free(conn_data); - free(num_pix_displ); return 1; } - - correct_rotation_and_stretch(connected, edet.det, conn_data, + correct_rotation_and_stretch(connected, det, conn_data, clen_to_use, stretch_coeff, gparams); - 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); - free(num_pix_displ); - return 1; - } + adjust_displ_for_stretch(connected, conn_data, stretch_coeff, + det, gpanels); - 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); + ret = compute_shift(connected, conn_data, det, gparams, gpanels); 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; } @@ -2897,9 +2414,7 @@ int optimize_geometry(struct geoptimiser_params *gparams, correct_shift(connected, conn_data, clen_to_use); - recompute_avg_displ(connected, conn_data, - num_pix_displ, &edet, - avg_displ); + recompute_avg_displ(connected, conn_data, det, gpanels); if ( gparams->error_maps ) { @@ -2907,16 +2422,10 @@ int optimize_geometry(struct geoptimiser_params *gparams, STATUS("Saving error map after correction.\n"); - ret = save_data_to_png("error_map_after.png", &edet, - avg_displ->displ_abs); - if ( ret !=0 ) { + ret = save_data_to_png("error_map_after.png", det, gpanels); + if ( ret ) { ERROR("Error while writing data to file.\n"); - - 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; } @@ -2930,14 +2439,13 @@ int optimize_geometry(struct geoptimiser_params *gparams, } STATUS("Computing errors after correction.\n"); - total_error = compute_error(connected, &edet, conn_data, - num_pix_displ, avg_displ->displ_abs); + total_error = compute_error(connected, det, conn_data, gpanels); STATUS("Detector-wide error after correction: RMSD = %0.4f pixels.\n", total_error); write_ret = write_detector_geometry_2(gparams->geometry_filename, - gparams->outfile, edet.det, + gparams->outfile, det, gparams->command_line, 1); if ( write_ret != 0 ) { ERROR("Error in writing output geometry file.\n"); @@ -2954,11 +2462,7 @@ int optimize_geometry(struct geoptimiser_params *gparams, #endif /* HAVE_SAVE_TO_PNG */ } - free_avg_displacement(avg_displ); - free_pixel_maps(pixel_maps); free(conn_data); - free(num_pix_displ); - free_pixel_maps(recomputed_pixel_maps); return 0; } |