From 1a206036c00ca86270fe59bfc475dd059bf38969 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sun, 6 Feb 2011 21:20:01 +0100 Subject: Start work on binary tree --- src/calibrate_detector.c | 3 +- src/cubeit.c | 3 +- src/geometry.c | 52 +++++++----------- src/geometry.h | 8 ++- src/image.h | 33 +---------- src/indexamajig.c | 8 +-- src/partialator.c | 81 ++++++++++++++------------- src/pattern_sim.c | 2 +- src/peaks.c | 140 ++++++++++++----------------------------------- src/post-refinement.c | 127 ++++++++++++++++++++++-------------------- src/post-refinement.h | 5 +- src/reflist.c | 139 ++++++++++++++++++++++++++++++++++++++++++++++ src/reflist.h | 69 +++++++++++++++++++++++ src/reintegrate.c | 3 +- src/templates.c | 52 ++++++++++-------- 15 files changed, 419 insertions(+), 306 deletions(-) create mode 100644 src/reflist.c create mode 100644 src/reflist.h (limited to 'src') diff --git a/src/calibrate_detector.c b/src/calibrate_detector.c index 2252842f..9aeae803 100644 --- a/src/calibrate_detector.c +++ b/src/calibrate_detector.c @@ -168,8 +168,7 @@ static void add_image(void *args, int cookie) image.flags = NULL; image.indexed_cell = NULL; image.filename = pargs->filename; - image.cpeaks = NULL; - image.n_cpeaks = 0; + image.reflections = NULL; image.det = NULL; STATUS("%3i: Processing '%s'\n", cookie, pargs->filename); diff --git a/src/cubeit.c b/src/cubeit.c index 808e2e23..7b40f8e5 100644 --- a/src/cubeit.c +++ b/src/cubeit.c @@ -218,8 +218,7 @@ static void sum_image(void *pg, int cookie) image.flags = NULL; image.indexed_cell = NULL; image.filename = apargs->filename; - image.cpeaks = NULL; - image.n_cpeaks = 0; + image.reflections = NULL; image.det = pargs->det; STATUS("Processing '%s'\n", apargs->filename); diff --git a/src/geometry.c b/src/geometry.c index 9bb41d12..9f9b953c 100644 --- a/src/geometry.c +++ b/src/geometry.c @@ -22,6 +22,7 @@ #include "image.h" #include "peaks.h" #include "beam-parameters.h" +#include "reflist.h" #define MAX_CPEAKS (256 * 256) @@ -117,7 +118,7 @@ static double partiality(double r1, double r2, double r) static int check_reflection(struct image *image, double mres, int output, - struct cpeak *cpeaks, int np, + RefList *reflections, signed int h, signed int k, signed int l, double asx, double asy, double asz, double bsx, double bsy, double bsz, @@ -136,6 +137,7 @@ static int check_reflection(struct image *image, double mres, int output, double divergence = image->div; double lambda = image->lambda; double klow, kcen, khigh; /* Wavenumber */ + Reflection *refl; /* "low" gives the largest Ewald sphere, * "high" gives the smallest Ewald sphere. */ @@ -207,17 +209,9 @@ static int check_reflection(struct image *image, double mres, int output, if ( p == -1 ) return 0; /* Add peak to list */ - cpeaks[np].h = h; - cpeaks[np].k = k; - cpeaks[np].l = l; - cpeaks[np].x = xda; - cpeaks[np].y = yda; - cpeaks[np].r1 = rlow; - cpeaks[np].r2 = rhigh; - cpeaks[np].p = part; - cpeaks[np].clamp1 = clamp_low; - cpeaks[np].clamp2 = clamp_high; - np++; + refl = add_refl(reflections, h, k, l); + set_detector_pos(refl, 0.0, xda, yda); + set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high); if ( output ) { printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %5.2f\n", @@ -228,23 +222,18 @@ static int check_reflection(struct image *image, double mres, int output, } -struct cpeak *find_intersections(struct image *image, UnitCell *cell, - int *n, int output) +RefList *find_intersections(struct image *image, UnitCell *cell, + int output) { double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; - struct cpeak *cpeaks; - int np = 0; + RefList *reflections; int hmax, kmax, lmax; double mres; signed int h, k, l; - cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS); - if ( cpeaks == NULL ) { - *n = 0; - return NULL; - } + reflections = reflist_new(); cell_get_reciprocal(cell, &asx, &asy, &asz, &bsx, &bsy, &bsz, @@ -260,29 +249,30 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell, for ( l=-lmax; l #endif -extern struct cpeak *find_intersections(struct image *image, UnitCell *cell, - int *n, int output); +#include "reflist.h" -extern double integrate_all(struct image *image, struct cpeak *cpeaks, int n); +extern RefList *find_intersections(struct image *image, UnitCell *cell, + int output); + +extern double integrate_all(struct image *image, RefList *reflections); #endif /* GEOMETRY_H */ diff --git a/src/image.h b/src/image.h index 10fb2406..ce0e7343 100644 --- a/src/image.h +++ b/src/image.h @@ -24,6 +24,7 @@ #include "utils.h" #include "cell.h" #include "detector.h" +#include "reflist.h" #define MAX_CELL_CANDIDATES (32) @@ -52,35 +53,6 @@ struct imagefeature { typedef struct _imagefeaturelist ImageFeatureList; -/* This structure represents a predicted peak in an image */ -struct cpeak -{ - /* Indices */ - signed int h; - signed int k; - signed int l; - - double min_distance; - int valid; - - /* Partiality */ - double r1; /* First excitation error */ - double r2; /* Second excitation error */ - double p; /* Partiality */ - int clamp1; /* Clamp status for r1 */ - int clamp2; /* Clamp status for r2 */ - - /* Location in image */ - int x; - int y; - - int scalable; - - /* Intensity */ - double intensity; -}; - - /* Structure describing an image */ struct image { @@ -93,8 +65,7 @@ struct image { struct detector *det; struct beam_params *beam; /* The nominal beam parameters */ char *filename; - struct cpeak *cpeaks; - int n_cpeaks; + RefList *reflections; int id; /* ID number of the thread * handling this image */ diff --git a/src/indexamajig.c b/src/indexamajig.c index 6c01f95e..40d7cbd8 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -263,8 +263,7 @@ static struct image *get_simage(struct image *template, int alternate) image->f0 = template->f0; /* Prevent muppetry */ - image->cpeaks = NULL; - image->n_cpeaks = 0; + image->reflections = NULL; return image; } @@ -326,8 +325,7 @@ static void process_image(void *pp, int cookie) image.indexed_cell = NULL; image.id = cookie; image.filename = filename; - image.cpeaks = NULL; - image.n_cpeaks = 0; + image.reflections = NULL; image.det = pargs->static_args.det; STATUS("Processing '%s'\n", image.filename); @@ -434,7 +432,7 @@ done: free(image.data); free(image.flags); image_feature_list_free(image.features); - free(image.cpeaks); + reflist_free(image.reflections); hdfile_close(hdfile); } diff --git a/src/partialator.c b/src/partialator.c index d1ad3768..a4f3381c 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -34,6 +34,7 @@ #include "beam-parameters.h" #include "post-refinement.h" #include "hrs-scaling.h" +#include "reflist.h" /* Maximum number of iterations of NLSq to do for each image per macrocycle. */ @@ -84,9 +85,9 @@ static void refine_image(int mytask, void *tasks) struct image *image = pargs->image; double nominal_photon_energy = pargs->image->beam->photon_energy; struct hdfile *hdfile; - struct cpeak *spots; - int n, i; + int i; double dev, last_dev; + RefList *reflections; hdfile = hdfile_open(image->filename); if ( hdfile == NULL ) { @@ -110,16 +111,17 @@ static void refine_image(int mytask, void *tasks) a/1.0e-9, b/1.0e-9, c/1.0e-9, rad2deg(al), rad2deg(be), rad2deg(ga)); - spots = find_intersections(image, image->indexed_cell, &n, 0); + /* FIXME: Don't do this each time */ + reflections = find_intersections(image, image->indexed_cell, 0); dev = +INFINITY; i = 0; do { last_dev = dev; - dev = pr_iterate(image, pargs->i_full, pargs->sym, &spots, &n); + dev = pr_iterate(image, pargs->i_full, pargs->sym, reflections); STATUS("Iteration %2i: mean dev = %5.2f\n", i, dev); i++; } while ( (fabs(last_dev - dev) > 1.0) && (i < MAX_CYCLES) ); - mean_partial_dev(image, spots, n, pargs->sym, + mean_partial_dev(image, reflections, pargs->sym, pargs->i_full, pargs->graph); if ( pargs->pgraph ) { fprintf(pargs->pgraph, "%5i %5.2f\n", mytask, dev); @@ -128,7 +130,7 @@ static void refine_image(int mytask, void *tasks) free(image->data); if ( image->flags != NULL ) free(image->flags); hdfile_close(hdfile); - free(spots); + reflist_free(reflections); /* Muppet proofing */ image->data = NULL; @@ -163,22 +165,23 @@ static void refine_all(struct image *images, int n_total_patterns, } -static void uniquify(struct cpeak *spot, const char *sym) +static void uniquify(Reflection *refl, const char *sym) { + signed int h, k, l; signed int ha, ka, la; - get_asymm(spot->h, spot->k, spot->l, &ha, &ka, &la, sym); - spot->h = ha; - spot->k = ka; - spot->l = la; + get_indices(refl, &h, &k, &l); + get_asymm(h, k, l, &ha, &ka, &la, sym); + set_indices(refl, h, k, l); } +/* FIXME: Get rid of this */ static void integrate_image(struct image *image, ReflItemList *obs, const char *sym) { - struct cpeak *spots; - int j, n; + RefList *reflections; + Reflection *refl; struct hdfile *hdfile; double nominal_photon_energy = image->beam->photon_energy; @@ -199,42 +202,39 @@ static void integrate_image(struct image *image, ReflItemList *obs, } /* Figure out which spots should appear in this pattern */ - spots = find_intersections(image, image->indexed_cell, &n, 0); + reflections = find_intersections(image, image->indexed_cell, 0); /* For each reflection, estimate the partiality */ - for ( j=0; jcpeaks = spots; - image->n_cpeaks = n; + image->reflections = reflections; free(image->data); if ( image->flags != NULL ) free(image->flags); @@ -253,21 +253,20 @@ static void select_scalable_reflections(struct image *images, int n) for ( m=0; mlambda); @@ -535,44 +531,26 @@ int find_projected_peaks(struct image *image, UnitCell *cell, if ( dist > domain_r ) continue; } - for ( j=0; jcpeaks = cpeaks; - image->n_cpeaks = n_cpeaks; + optimise_reflist(reflections); + + STATUS("Found %i reflections\n", n_reflections); + image->reflections = reflections; - return n_cpeaks; + return n_reflections; } @@ -685,20 +663,14 @@ void output_intensities(struct image *image, UnitCell *cell, int use_closer, FILE *ofh, int circular_domain, double domain_r) { - int i; - int n_found; - int n_indclose = 0; - int n_foundclose = 0; - int n_veto = 0; - int n_veto_second = 0; double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; + Reflection *refl; - if ( image->n_cpeaks == 0 ) { + if ( image->reflections == NULL ) { find_projected_peaks(image, cell, circular_domain, domain_r); } - if ( image->n_cpeaks == 0 ) return; /* Get exclusive access to the output stream if necessary */ if ( mutex != NULL ) pthread_mutex_lock(mutex); @@ -709,16 +681,20 @@ void output_intensities(struct image *image, UnitCell *cell, &bsx, &bsy, &bsz, &csx, &csy, &csz); - for ( i=0; in_cpeaks; i++ ) { + for ( refl = first_refl(image->reflections); + refl != NULL; + refl = next_refl(refl) ) { float x, y, intensity; double d; int idx; double bg, max; struct panel *p; + double px, py; + signed int h, k, l; - p = find_panel(image->det, image->cpeaks[i].x, - image->cpeaks[i].y); + get_detector_pos(refl, &px, &py); + p = find_panel(image->det, px, py); if ( p == NULL ) continue; if ( p->no_index ) continue; @@ -729,9 +705,7 @@ void output_intensities(struct image *image, UnitCell *cell, if ( image->features != NULL ) { f = image_feature_closest(image->features, - image->cpeaks[i].x, - image->cpeaks[i].y, - &d, &idx); + px, py, &d, &idx); } else { f = NULL; } @@ -754,7 +728,6 @@ void output_intensities(struct image *image, UnitCell *cell, * region. Integration of the same * peak included a bad region this time. */ - n_veto_second++; continue; } intensity = f->intensity; @@ -763,81 +736,36 @@ void output_intensities(struct image *image, UnitCell *cell, int r; - r = integrate_peak(image, - image->cpeaks[i].x, - image->cpeaks[i].y, - &x, &y, &intensity, &bg, &max, + r = integrate_peak(image, px, py, &x, &y, + &intensity, &bg, &max, polar, 1); if ( r ) { /* Plain old ordinary peak veto */ - n_veto++; continue; } } - if ( (f != NULL) && (d < PEAK_CLOSE) ) { - n_indclose++; - } - } else { int r; - r = integrate_peak(image, - image->cpeaks[i].x, - image->cpeaks[i].y, - &x, &y, &intensity, &bg, &max, polar, - 0); + r = integrate_peak(image, px, py, &x, &y, + &intensity, &bg, &max, polar, 0); if ( r ) { /* Plain old ordinary peak veto */ - n_veto++; continue; } } /* Write h,k,l, integrated intensity and centroid coordinates */ + get_indices(refl, &h, &k, &l); fprintf(ofh, "%3i %3i %3i %6f (at %5.2f,%5.2f) max=%6f bg=%6f\n", - image->cpeaks[i].h, image->cpeaks[i].k, image->cpeaks[i].l, - intensity, x, y, max, bg); - - image->cpeaks[i].x = x; - image->cpeaks[i].y = y; - - } - - n_found = image_feature_count(image->features); - for ( i=0; ifeatures, i); - - if ( f == NULL ) continue; - - for ( j=0; jn_cpeaks; j++ ) { - - double d; - - d = pow(image->cpeaks[j].x-f->x, 2.0) - + pow(image->cpeaks[j].y-f->y, 2.0); - - if ( d < PEAK_CLOSE ) n_foundclose++; - - } + h, l, l, intensity, x, y, max, bg); } - fprintf(ofh, "Peak statistics: %i peaks found by the peak search out of " - "%i were close to indexed positions. " - "%i indexed positions out of %i were close to detected peaks.\n", - n_foundclose, n_found, n_indclose, image->n_cpeaks); - fprintf(ofh, "%i integrations using indexed locations were aborted because " - "they hit one or more bad pixels.\n", n_veto); - fprintf(ofh, "%i integrations using peak search locations were aborted " - "because they hit one or more bad pixels.\n", n_veto_second); /* Blank line at end */ fprintf(ofh, "\n"); diff --git a/src/post-refinement.c b/src/post-refinement.c index b860bb42..25664e69 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -65,7 +65,7 @@ static double partiality_rgradient(double r, double profile_radius) /* Return the gradient of parameter 'k' given the current status of 'image'. */ -double gradient(struct image *image, int k, struct cpeak spot, double r) +double gradient(struct image *image, int k, Reflection *refl, double r) { double ds, tt, azi; double nom, den; @@ -74,24 +74,32 @@ double gradient(struct image *image, int k, struct cpeak spot, double r) double bsx, bsy, bsz; double csx, csy, csz; double xl, yl, zl; + signed int hi, ki, li; + double r1, r2, p; + int clamp_low, clamp_high; + + get_indices(refl, &hi, &ki, &li); cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); - xl = spot.h*asx + spot.k*bsx + spot.l*csx; - yl = spot.h*asy + spot.k*bsy + spot.l*csy; - zl = spot.h*asz + spot.k*bsz + spot.l*csz; + xl = hi*asx + ki*bsx + li*csx; + yl = hi*asy + ki*bsy + li*csy; + zl = hi*asz + ki*bsz + li*csz; + - ds = 2.0 * resolution(image->indexed_cell, spot.h, spot.k, spot.l); + ds = 2.0 * resolution(image->indexed_cell, hi, ki, li); tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+k); azi = angle_between(1.0, 0.0, 0.0, xl, yl, 0.0); + get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high); + /* Calculate the gradient of partiality wrt excitation error. */ - if ( spot.clamp1 == 0 ) { - g += partiality_gradient(spot.r1, r); + if ( clamp_low == 0 ) { + g += partiality_gradient(r1, r); } - if ( spot.clamp2 == 0 ) { - g += partiality_gradient(spot.r2, r); + if ( clamp_high == 0 ) { + g += partiality_gradient(r2, r); } /* For many gradients, just multiply the above number by the gradient @@ -99,7 +107,7 @@ double gradient(struct image *image, int k, struct cpeak spot, double r) switch ( k ) { case REF_SCALE : - return -spot.p*pow(image->osf, -2.0); + return -p*pow(image->osf, -2.0); case REF_DIV : nom = sqrt(2.0) * ds * sin(image->div/2.0); @@ -107,33 +115,33 @@ double gradient(struct image *image, int k, struct cpeak spot, double r) return (nom/den) * g; case REF_R : - if ( spot.clamp1 == 0 ) { - g += partiality_rgradient(spot.r1, r); + if ( clamp_low == 0 ) { + g += partiality_rgradient(r1, r); } - if ( spot.clamp2 == 0 ) { - g += partiality_rgradient(spot.r2, r); + if ( clamp_high == 0 ) { + g += partiality_rgradient(r2, r); } return g; /* Cell parameters and orientation */ case REF_ASX : - return spot.h * pow(sin(tt), -1.0) * cos(azi) * g; + return hi * pow(sin(tt), -1.0) * cos(azi) * g; case REF_BSX : - return spot.k * pow(sin(tt), -1.0) * cos(azi) * g; + return ki * pow(sin(tt), -1.0) * cos(azi) * g; case REF_CSX : - return spot.l * pow(sin(tt), -1.0) * cos(azi) * g; + return li * pow(sin(tt), -1.0) * cos(azi) * g; case REF_ASY : - return spot.h * pow(sin(tt), -1.0) * sin(azi) * g; + return hi * pow(sin(tt), -1.0) * sin(azi) * g; case REF_BSY : - return spot.k * pow(sin(tt), -1.0) * sin(azi) * g; + return ki * pow(sin(tt), -1.0) * sin(azi) * g; case REF_CSY : - return spot.l * pow(sin(tt), -1.0) * sin(azi) * g; + return li * pow(sin(tt), -1.0) * sin(azi) * g; case REF_ASZ : - return spot.h * pow(cos(tt), -1.0) * g; + return hi * pow(cos(tt), -1.0) * g; case REF_BSZ : - return spot.k * pow(cos(tt), -1.0) * g; + return ki * pow(cos(tt), -1.0) * g; case REF_CSZ : - return spot.l * pow(cos(tt), -1.0) * g; + return li * pow(cos(tt), -1.0) * g; } @@ -218,31 +226,34 @@ void apply_shift(struct image *image, int k, double shift) } -double mean_partial_dev(struct image *image, struct cpeak *spots, int n, +double mean_partial_dev(struct image *image, RefList *reflections, const char *sym, double *i_full, FILE *graph) { - int h, n_used; + int n_used; double delta_I = 0.0; + Reflection *refl; n_used = 0; - for ( h=0; hosf)); + p = get_partiality(refl); + delta_I += fabs(I_partial - (p * I_full / image->osf)); n_used++; if ( graph != NULL ) { fprintf(graph, "%3i %3i %3i %5.2f (at %5.2f,%5.2f)" " %5.2f %5.2f\n", - hind, kind, lind, I_partial/spots[h].p, xc, yc, - spots[h].p, I_partial / I_full); + hind, kind, lind, I_partial/p, xc, yc, + p, I_partial / I_full); } } @@ -268,19 +280,21 @@ double mean_partial_dev(struct image *image, struct cpeak *spots, int n, /* Perform one cycle of post refinement on 'image' against 'i_full' */ double pr_iterate(struct image *image, double *i_full, const char *sym, - struct cpeak **pspots, int *n) + RefList *reflections) { gsl_matrix *M; gsl_vector *v; gsl_vector *shifts; - int h, param; - struct cpeak *spots = *pspots; + int param; + Reflection *refl; M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS); v = gsl_vector_calloc(NUM_PARAMS); /* Construct the equations, one per reflection in this image */ - for ( h=0; h<*n; h++ ) { + for ( refl = first_refl(reflections); + refl != NULL; + refl = next_refl(refl) ) { signed int hind, kind, lind; signed int ha, ka, la; @@ -288,17 +302,17 @@ double pr_iterate(struct image *image, double *i_full, const char *sym, float I_partial; float xc, yc; int k; + double x, y; + double p; - hind = spots[h].h; - kind = spots[h].k; - lind = spots[h].l; + get_indices(refl, &hind, &kind, &lind); - if ( !spots[h].scalable ) continue; + if ( !get_scalable(refl) ) continue; - /* Actual measurement of this reflection from this - * pattern? */ - /* FIXME: Coordinates aren't whole numbers */ - if ( integrate_peak(image, spots[h].x, spots[h].y, + /* Actual measurement of this reflection from this pattern? */ + get_detector_pos(refl, &x, &y); + /* FIXME: Get rid of this */ + if ( integrate_peak(image, x, y, &xc, &yc, &I_partial, NULL, NULL, 1, 0) ) { continue; @@ -306,7 +320,8 @@ double pr_iterate(struct image *image, double *i_full, const char *sym, get_asymm(hind, kind, lind, &ha, &ka, &la, sym); I_full = lookup_intensity(i_full, ha, ka, la); - delta_I = I_partial - (spots[h].p * I_full / image->osf); + p = get_partiality(refl); + delta_I = I_partial - (p * I_full / image->osf); for ( k=0; kprofile_radius) - * gradient(image, k, spots[h], + * gradient(image, k, refl, image->profile_radius); M_c *= pow(I_full, 2.0); @@ -329,8 +344,7 @@ double pr_iterate(struct image *image, double *i_full, const char *sym, } - gr = gradient(image, k, spots[h], - image->profile_radius); + gr = gradient(image, k, refl, image->profile_radius); v_c = delta_I * I_full * gr; gsl_vector_set(v, k, v_c); @@ -350,8 +364,5 @@ double pr_iterate(struct image *image, double *i_full, const char *sym, gsl_vector_free(v); gsl_vector_free(shifts); - free(spots); - spots = find_intersections(image, image->indexed_cell, n, 0); - *pspots = spots; - return mean_partial_dev(image, spots, *n, sym, i_full, NULL); + return mean_partial_dev(image, reflections, sym, i_full, NULL); } diff --git a/src/post-refinement.h b/src/post-refinement.h index 983ecc28..b3c6c074 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -21,6 +21,7 @@ #include #include "image.h" +#include "reflist.h" /* Refineable parameters */ @@ -44,12 +45,12 @@ enum { void apply_shift(struct image *image, int k, double shift); -double mean_partial_dev(struct image *image, struct cpeak *spots, int n, +double mean_partial_dev(struct image *image, RefList *reflections, const char *sym, double *i_full, FILE *graph); double pr_iterate(struct image *image, double *i_full, const char *sym, - struct cpeak **pspots, int *n); + RefList *reflections); #endif /* POST_REFINEMENT_H */ diff --git a/src/reflist.c b/src/reflist.c new file mode 100644 index 00000000..990bc535 --- /dev/null +++ b/src/reflist.c @@ -0,0 +1,139 @@ +/* + * reflist.c + * + * Fast reflection/peak list + * + * (c) 2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#include + +#include "reflist.h" + + +struct _reflection { + + unsigned int serial; /* Serial number */ + + signed int h; + signed int k; + signed int l; + + double excitation_error; + + /* Partiality */ + double r1; /* First excitation error */ + double r2; /* Second excitation error */ + double p; /* Partiality */ + int clamp1; /* Clamp status for r1 */ + int clamp2; /* Clamp status for r2 */ + + /* Location in image */ + int x; + int y; + + int scalable; + + /* Intensity */ + double intensity; + + +}; + + +struct _reflist { + + struct ref *head; + +}; + + +/**************************** Creation / deletion *****************************/ + +/* Create a reflection list */ +RefList *reflist_new() +{ + RefList *new; + + new = malloc(sizeof(struct _reflist)); + + return new; +} + + +void reflist_free(RefList *list) +{ +} + + +/********************************** Search ************************************/ + +Reflection *find_refl(RefList *list, INDICES) +{ +} + + +/********************************** Getters ***********************************/ + +double get_excitation_error(Reflection *refl) +{ +} + + +void get_detector_pos(Reflection *refl, double *x, double *y) +{ +} + + +void get_indices(Reflection *refl, signed int *h, signed int *k, signed int *l) +{ +} + + +/********************************** Setters ***********************************/ + +void set_detector_pos(Reflection *refl, double exerr, double x, double y) +{ +} + + +void set_partial(Reflection *refl, double r1, double r2, double p, + double clamp_low, double clamp_high) +{ +} + + +/********************************* Insertion **********************************/ + +Reflection *add_refl(RefList *list, INDICES) +{ +} + + +Reflection *add_refl_with_det_pos(RefList *refl, INDICES, double exerr, + double x, double y) +{ +} + + +/********************************* Iteration **********************************/ + +Reflection *first_refl(RefList *list) +{ +} + + +Reflection *next_refl(Reflection *refl) +{ +} + + +/*********************************** Voodoo ***********************************/ + +void optimise_reflist(RefList *list) +{ +} diff --git a/src/reflist.h b/src/reflist.h new file mode 100644 index 00000000..0ba09a3b --- /dev/null +++ b/src/reflist.h @@ -0,0 +1,69 @@ +/* + * reflist.h + * + * Fast reflection/peak list + * + * (c) 2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifndef REFLIST_H +#define REFLIST_H + +#ifdef HAVE_CONFIG_H +#include +#endif + + +typedef struct _reflist RefList; +typedef struct _reflection Reflection; + +#define INDICES signed int h, signed int k, signed int l + +/* Creation/deletion */ +extern RefList *reflist_new(void); +extern void reflist_free(RefList *list); + +/* Search */ +extern Reflection *find_refl(RefList *list, INDICES); + +/* Get */ +extern double get_excitation_error(Reflection *refl); +extern void get_detector_pos(Reflection *refl, double *x, double *y); +extern void get_indices(Reflection *refl, + signed int *h, signed int *k, signed int *l); +extern double get_partiality(Reflection *refl); +extern double get_intensity(Reflection *refl); +extern void get_partial(Reflection *refl, double *r1, double *r2, double *p, + int *clamp_low, int *clamp_high); +extern int get_scalable(Reflection *refl); + +/* Set */ +extern void set_detector_pos(Reflection *refl, double exerr, + double x, double y); +extern void set_partial(Reflection *refl, double r1, double r2, double p, + double clamp_low, double clamp_high); +extern void set_indices(Reflection *refl, + signed int h, signed int k, signed int l); +extern void set_int(Reflection *refl, double intensity); +extern void set_scalable(Reflection *refl, int scalable); + +/* Insertion */ +extern Reflection *add_refl(RefList *list, INDICES); +extern Reflection *add_refl_with_det_pos(RefList *list, INDICES, double exerr, + double x, double y); + + +/* Deletion */ +extern void delete_refl(Reflection *refl); + +/* Iteration */ +extern Reflection *first_refl(RefList *list); +extern Reflection *next_refl(Reflection *refl); + +/* Voodoo */ +extern void optimise_reflist(RefList *list); + +#endif /* REFLIST_H */ diff --git a/src/reintegrate.c b/src/reintegrate.c index 9d882d4b..09634cd7 100644 --- a/src/reintegrate.c +++ b/src/reintegrate.c @@ -110,8 +110,7 @@ static void process_image(void *pg, int cookie) image.flags = NULL; image.indexed_cell = NULL; image.filename = apargs->filename; - image.cpeaks = NULL; - image.n_cpeaks = 0; + image.reflections = NULL; image.det = pargs->det; STATUS("Processing '%s'\n", apargs->filename); diff --git a/src/templates.c b/src/templates.c index 4bd85a76..9abbcb17 100644 --- a/src/templates.c +++ b/src/templates.c @@ -21,6 +21,7 @@ #include "geometry.h" #include "hdf5-file.h" #include "peaks.h" +#include "reflist.h" #include @@ -42,8 +43,7 @@ struct _indexingprivate_template struct template { double omega; double phi; - int n; - struct cpeak *spots; + RefList *spots; }; @@ -164,24 +164,22 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename, for ( omega = 0.0; omega < omega_max-omega_step; omega += omega_step ) { for ( phi = 0.0; phi < phi_max-phi_step; phi += phi_step ) { - int n; - struct cpeak *cpeaks; + RefList *reflections; UnitCell *cell_rot; assert(i < n_templates); cell_rot = rotate_cell(cell, omega, phi, 0.0); - cpeaks = find_intersections(&image, cell_rot, &n, 0); - if ( cpeaks == NULL ) { + reflections = find_intersections(&image, cell_rot, 0); + if ( reflections == NULL ) { ERROR("Template calculation failed.\n"); return NULL; } priv->templates[i].omega = omega; priv->templates[i].phi = phi; - priv->templates[i].n = n; - priv->templates[i].spots = cpeaks; + priv->templates[i].spots = reflections; i++; free(cell_rot); @@ -227,23 +225,28 @@ static int fast_integrate_peak(struct image *image, int xp, int yp) } -static double integrate_all_rot(struct image *image, struct cpeak *cpeaks, - int n, double rot) +static double integrate_all_rot(struct image *image, RefList *reflections, + double rot) { double itot = 0.0; - int i; double cosr, sinr; int num_int = 0; + Reflection *refl; cosr = cos(rot); sinr = sin(rot); - for ( i=0; itemplates[i].spots, - priv->templates[i].n, rot); + rot); best = val > max; } else { val = mean_distance(image, priv->templates[i].spots, - priv->templates[i].n, rot); + rot); best = val < max; } -- cgit v1.2.3