aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--libcrystfel/src/image.h3
-rw-r--r--src/geoptimiser.c1782
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;
}