diff options
Diffstat (limited to 'libcrystfel/src/reflist.c')
| -rw-r--r-- | libcrystfel/src/reflist.c | 1048 |
1 files changed, 1048 insertions, 0 deletions
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c new file mode 100644 index 00000000..18856c67 --- /dev/null +++ b/libcrystfel/src/reflist.c @@ -0,0 +1,1048 @@ +/* + * reflist.c + * + * Fast reflection/peak list + * + * (c) 2011 Thomas White <taw@physics.org> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#include <stdlib.h> +#include <assert.h> +#include <stdio.h> +#include <pthread.h> + +#include "reflist.h" +#include "utils.h" + +/** + * SECTION:reflist + * @short_description: The fast reflection list + * @title: RefList + * @section_id: + * @see_also: + * @include: "reflist.h" + * @Image: + * + * The fast reflection list stores reflections in an RB-tree indexed + * by the Miller indices h, k and l. Any reflection can be found in a + * length of time which scales logarithmically with the number of reflections in + * the list. + * + * A RefList can contain any number of reflections, and can store more than + * one reflection with a given set of indices, for example when two distinct + * reflections are to be stored according to their asymmetric indices. + * + * There are getters and setters which can be used to get and set values for an + * individual reflection. The reflection list does not calculate any values, + * only stores what it was given earlier. As such, you will need to carefully + * examine which fields your prior processing steps have filled in. + */ + + +struct _refldata { + + /* Symmetric indices (i.e. the "real" indices) */ + signed int hs; + signed int ks; + signed int ls; + + /* Partiality and related geometrical stuff */ + 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 */ + double fs; + double ss; + + /* The distance from the exact Bragg position to the coordinates + * given above. */ + double excitation_error; + + /* Non-zero if this reflection can be used for scaling */ + int scalable; + + /* Non-zero if this reflection should be used as a "guide star" for + * post refinement */ + int refinable; + + /* Intensity */ + double intensity; + double esd_i; + + /* Phase */ + double phase; + int have_phase; + + /* Redundancy */ + int redundancy; + + /* User-specified temporary values */ + double temp1; + double temp2; +}; + + +enum _nodecol { + RED, + BLACK +}; + + +struct _reflection { + + /* Listy stuff */ + unsigned int serial; /* Unique serial number, key */ + struct _reflection *child[2]; /* Child nodes */ + struct _reflection *next; /* Next and previous in doubly linked */ + struct _reflection *prev; /* list of duplicate reflections */ + enum _nodecol col; /* Colour (red or black) */ + + /* Payload */ + pthread_mutex_t lock; /* Protects the contents of "data" */ + struct _refldata data; +}; + + +struct _reflist { + + struct _reflection *head; + struct _reflection *tail; + +}; + + +#define SERIAL(h, k, l) ((((h)+256)<<18) + (((k)+256)<<9) + ((l)+256)) +#define GET_H(serial) ((((serial) & 0xfffc0000)>>18)-256) +#define GET_K(serial) ((((serial) & 0x0003fe00)>>9)-256) +#define GET_L(serial) (((serial) & 0x000001ff)-256) + +/**************************** Creation / deletion *****************************/ + +static Reflection *new_node(unsigned int serial) +{ + Reflection *new; + + new = calloc(1, sizeof(struct _reflection)); + new->serial = serial; + new->next = NULL; + new->prev = NULL; + new->child[0] = NULL; + new->child[1] = NULL; + new->col = RED; + pthread_mutex_init(&new->lock, NULL); + + return new; +} + + +/** + * reflist_new: + * + * Creates a new reflection list. + * + * Returns: the new reflection list, or NULL on error. + */ +RefList *reflist_new() +{ + RefList *new; + + new = malloc(sizeof(struct _reflist)); + if ( new == NULL ) return NULL; + + new->head = NULL; + + return new; +} + + +/** + * reflection_new: + * @h: The h index of the new reflection + * @k: The k index of the new reflection + * @l: The l index of the new reflection + * + * Creates a new individual reflection. You'll probably want to use + * add_refl_to_list() at some later point. + */ +Reflection *reflection_new(signed int h, signed int k, signed int l) +{ + return new_node(SERIAL(h, k, l)); +} + + +/** + * reflection_free: + * @refl: The reflection to free. + * + * Destroys an individual reflection. + */ +void reflection_free(Reflection *refl) +{ + pthread_mutex_destroy(&refl->lock); + free(refl); +} + + +static void recursive_free(Reflection *refl) +{ + if ( refl->child[0] != NULL ) recursive_free(refl->child[0]); + if ( refl->child[1] != NULL ) recursive_free(refl->child[1]); + + while ( refl != NULL ) { + Reflection *next = refl->next; + reflection_free(refl); + refl = next; + } +} + + +/** + * reflist_free: + * @list: The reflection list to free. + * + * Destroys a reflection list. + */ +void reflist_free(RefList *list) +{ + if ( list == NULL ) return; + if ( list->head != NULL ) { + recursive_free(list->head); + } /* else empty list */ + free(list); +} + + +/********************************** Search ************************************/ + +/** + * find_refl: + * @list: The reflection list to search in + * @h: The 'h' index to search for + * @k: The 'k' index to search for + * @l: The 'l' index to search for + * + * This function finds the first reflection in 'list' with the given indices. + * + * Since a %RefList can contain multiple reflections with the same indices, you + * may need to use next_found_refl() to get the other reflections. + * + * Returns: The found reflection, or NULL if no reflection with the given + * indices could be found. + **/ +Reflection *find_refl(const RefList *list, + signed int h, signed int k, signed int l) +{ + unsigned int search = SERIAL(h, k, l); + Reflection *refl; + + if ( list->head == NULL ) return NULL; + + /* Indices greater than or equal to 256 are filtered out when + * reflections are added, so don't even bother looking. + * (also, looking for such reflections causes trouble because the search + * serial number would be invalid) */ + if ( abs(h) >= 256 ) return NULL; + if ( abs(k) >= 256 ) return NULL; + if ( abs(l) >= 256 ) return NULL; + + refl = list->head; + + while ( refl != NULL ) { + + if ( refl->serial == search ) { + + assert(search == refl->serial); + assert(h == GET_H(refl->serial)); + assert(k == GET_K(refl->serial)); + assert(l == GET_L(refl->serial)); + return refl; + + } else { + + int dir = search > refl->serial; + if ( refl->child[dir] != NULL ) { + refl = refl->child[dir]; + } else { + /* Hit the bottom of the tree */ + return NULL; + } + + } + + } + + return NULL; +} + + +/** + * next_found_refl: + * @refl: A reflection returned by find_refl() or next_found_refl() + * + * This function returns the next reflection in @refl's list with the same + * indices. + * + * Returns: The found reflection, or NULL if there are no more reflections with + * the same indices. + **/ +Reflection *next_found_refl(Reflection *refl) +{ + if ( refl->next != NULL ) assert(refl->serial == refl->next->serial); + + return refl->next; /* Well, that was easy... */ +} + + +/********************************** Getters ***********************************/ + +/** + * get_excitation_error: + * @refl: A %Reflection + * + * Returns: The excitation error for the reflection. + **/ +double get_excitation_error(const Reflection *refl) +{ + return refl->data.excitation_error; +} + + +/** + * get_detector_pos: + * @refl: A %Reflection + * @fs: Location at which to store the fast scan offset of the reflection + * @ss: Location at which to store the slow scan offset of the reflection + * + **/ +void get_detector_pos(const Reflection *refl, double *fs, double *ss) +{ + *fs = refl->data.fs; + *ss = refl->data.ss; +} + + +/** + * get_indices: + * @refl: A %Reflection + * @h: Location at which to store the 'h' index of the reflection + * @k: Location at which to store the 'k' index of the reflection + * @l: Location at which to store the 'l' index of the reflection + * + **/ +void get_indices(const Reflection *refl, + signed int *h, signed int *k, signed int *l) +{ + *h = GET_H(refl->serial); + *k = GET_K(refl->serial); + *l = GET_L(refl->serial); +} + + +/** + * get_symmetric_indices: + * @refl: A %Reflection + * @hs: Location at which to store the 'h' index of the reflection + * @ks: Location at which to store the 'k' index of the reflection + * @ls: Location at which to store the 'l' index of the reflection + * + * This function gives the symmetric indices, that is, the "real" indices before + * squashing down to the asymmetric reciprocal unit. This may be useful if the + * list is indexed according to the asymmetric indices, but you still need + * access to the symmetric version. This happens during post-refinement. + * + **/ +void get_symmetric_indices(const Reflection *refl, + signed int *hs, signed int *ks, + signed int *ls) +{ + *hs = refl->data.hs; + *ks = refl->data.ks; + *ls = refl->data.ls; +} + + +/** + * get_partiality: + * @refl: A %Reflection + * + * Returns: The partiality of the reflection. + **/ +double get_partiality(const Reflection *refl) +{ + return refl->data.p; +} + + +/** + * get_intensity: + * @refl: A %Reflection + * + * Returns: The intensity of the reflection. + **/ +double get_intensity(const Reflection *refl) +{ + return refl->data.intensity; +} + + +/** + * get_partial: + * @refl: A %Reflection + * @r1: Location at which to store the first excitation error + * @r2: Location at which to store the second excitation error + * @p: Location at which to store the partiality + * @clamp_low: Location at which to store the first clamp status + * @clamp_high: Location at which to store the second clamp status + * + * This function is used during post refinement (in conjunction with + * set_partial()) to get access to the details of the partiality calculation. + * + **/ +void get_partial(const Reflection *refl, double *r1, double *r2, double *p, + int *clamp_low, int *clamp_high) +{ + *r1 = refl->data.r1; + *r2 = refl->data.r2; + *p = get_partiality(refl); + *clamp_low = refl->data.clamp1; + *clamp_high = refl->data.clamp2; +} + + +/** + * get_scalable: + * @refl: A %Reflection + * + * Returns: non-zero if this reflection can be scaled. + * + **/ +int get_scalable(const Reflection *refl) +{ + return refl->data.scalable; +} + + +/** + * get_refinable: + * @refl: A %Reflection + * + * Returns: non-zero if this reflection can be used for post refinement. + * + **/ +int get_refinable(const Reflection *refl) +{ + return refl->data.refinable; +} + + +/** + * get_redundancy: + * @refl: A %Reflection + * + * The redundancy of the reflection is the number of measurements that have been + * made of it. Note that a redundancy of zero may have a special meaning, such + * as that the reflection was impossible to integrate. Note further that each + * reflection in the list has its own redundancy, even if there are multiple + * copies of the reflection in the list. The total number of reflection + * measurements should always be the sum of the redundancies in the entire list. + * + * Returns: the number of measurements of this reflection. + * + **/ +int get_redundancy(const Reflection *refl) +{ + return refl->data.redundancy; +} + + +/** + * get_esd_intensity: + * @refl: A %Reflection + * + * Returns: the standard error in the intensity measurement (as returned by + * get_intensity()) for this reflection. + * + **/ +double get_esd_intensity(const Reflection *refl) +{ + return refl->data.esd_i; +} + + +/** + * get_phase: + * @refl: A %Reflection + * @have_phase: Place to store a non-zero value if the phase is set, or NULL. + * + * Returns: the phase for this reflection. + * + **/ +double get_phase(const Reflection *refl, int *have_phase) +{ + if ( have_phase != NULL ) *have_phase = refl->data.have_phase; + return refl->data.phase; +} + + +/** + * get_temp1: + * @refl: A %Reflection + * + * The temporary values can be used according to the needs of the calling + * program. + * + * Returns: the first temporary value for this reflection. + * + **/ +double get_temp1(const Reflection *refl) +{ + return refl->data.temp1; +} + + +/** + * get_temp2: + * @refl: A %Reflection + * + * The temporary values can be used according to the needs of the calling + * program. + * + * Returns: the second temporary value for this reflection. + * + **/ +double get_temp2(const Reflection *refl) +{ + return refl->data.temp2; +} + + +/********************************** Setters ***********************************/ + +/** + * copy_data: + * @to: %Reflection to copy data into + * @from: %Reflection to copy data from + * + * This function is used to copy the data (which is everything listed above in + * the list of getters and setters, apart from the indices themselves) from one + * reflection to another. This might be used when creating a new list from an + * old one, perhaps using the asymmetric indices instead of the raw indices for + * the new list. + * + **/ +void copy_data(Reflection *to, const Reflection *from) +{ + memcpy(&to->data, &from->data, sizeof(struct _refldata)); +} + + +/** + * set_detector_pos: + * @refl: A %Reflection + * @exerr: The excitation error for this reflection + * @fs: The fast scan offset of the reflection + * @ss: The slow scan offset of the reflection + * + **/ +void set_detector_pos(Reflection *refl, double exerr, double fs, double ss) +{ + refl->data.excitation_error = exerr; + refl->data.fs = fs; + refl->data.ss = ss; +} + + +/** + * set_partial: + * @refl: A %Reflection + * @r1: The first excitation error + * @r2: The second excitation error + * @p: The partiality + * @clamp_low: The first clamp status + * @clamp_high: The second clamp status + * + * This function is used during post refinement (in conjunction with + * get_partial()) to get access to the details of the partiality calculation. + * + **/ +void set_partial(Reflection *refl, double r1, double r2, double p, + double clamp_low, double clamp_high) +{ + refl->data.r1 = r1; + refl->data.r2 = r2; + refl->data.p = p; + refl->data.clamp1 = clamp_low; + refl->data.clamp2 = clamp_high; +} + + +/** + * set_int: + * @refl: A %Reflection + * @intensity: The intensity for the reflection. + * + * Set the intensity for the reflection. Note that retrieval is done with + * get_intensity(). + **/ +void set_int(Reflection *refl, double intensity) +{ + refl->data.intensity = intensity; +} + + +/** + * set_scalable: + * @refl: A %Reflection + * @scalable: Non-zero if this reflection should be scaled. + * + **/ +void set_scalable(Reflection *refl, int scalable) +{ + refl->data.scalable = scalable; +} + + +/** + * set_refinable: + * @refl: A %Reflection + * @refinable: Non-zero if this reflection can be used for post refinement. + * + **/ +void set_refinable(Reflection *refl, int refinable) +{ + refl->data.refinable = refinable; +} + + +/** + * set_redundancy: + * @refl: A %Reflection + * @red: New redundancy for the reflection + * + * The redundancy of the reflection is the number of measurements that have been + * made of it. Note that a redundancy of zero may have a special meaning, such + * as that the reflection was impossible to integrate. Note further that each + * reflection in the list has its own redundancy, even if there are multiple + * copies of the reflection in the list. The total number of reflection + * measurements should always be the sum of the redundancies in the entire list. + * + **/ +void set_redundancy(Reflection *refl, int red) +{ + refl->data.redundancy = red; +} + + +/** + * set_esd_intensity: + * @refl: A %Reflection + * @esd: New standard error for this reflection's intensity measurement + * + **/ +void set_esd_intensity(Reflection *refl, double esd) +{ + refl->data.esd_i = esd; +} + + +/** + * set_ph: + * @refl: A %Reflection + * @phase: New phase for the reflection + * + **/ +void set_ph(Reflection *refl, double phase) +{ + refl->data.phase = phase; + refl->data.have_phase = 1; +} + + +/** + * set_symmetric_indices: + * @refl: A %Reflection + * @hs: The 'h' index of the reflection + * @ks: The 'k' index of the reflection + * @ls: The 'l' index of the reflection + * + * This function gives the symmetric indices, that is, the "real" indices before + * squashing down to the asymmetric reciprocal unit. This may be useful if the + * list is indexed according to the asymmetric indices, but you still need + * access to the symmetric version. This happens during post-refinement. + * + **/ +void set_symmetric_indices(Reflection *refl, + signed int hs, signed int ks, signed int ls) +{ + refl->data.hs = hs; + refl->data.ks = ks; + refl->data.ls = ls; +} + + +/** + * set_temp1 + * @refl: A %Reflection + * @temp: New temporary value for the reflection + * + * The temporary values can be used according to the needs of the calling + * program. + * + **/ +void set_temp1(Reflection *refl, double temp) +{ + refl->data.temp1 = temp; +} + + +/** + * set_temp2 + * @refl: A %Reflection + * @temp: New temporary value for the reflection + * + * The temporary values can be used according to the needs of the calling + * program. + * + **/ +void set_temp2(Reflection *refl, double temp) +{ + refl->data.temp2 = temp; +} + + +/********************************* Insertion **********************************/ + +static Reflection *rotate_once(Reflection *refl, int dir) +{ + Reflection *s = refl->child[!dir]; + + refl->child[!dir] = s->child[dir]; + s->child[dir] = refl; + + refl->col = RED; + s->col = BLACK; + + return s; +} + + +static Reflection *rotate_twice(Reflection *refl, int dir) +{ + refl->child[!dir] = rotate_once(refl->child[!dir], !dir); + return rotate_once(refl, dir); +} + + +static int is_red(Reflection *refl) +{ + return (refl != NULL) && (refl->col == RED); +} + + +static Reflection *insert_node(Reflection *refl, Reflection *new) +{ + if ( refl == NULL ) { + + refl = new; + + } else { + + int dir; + Reflection *rcd; + + assert(new->serial != refl->serial); + dir = new->serial > refl->serial; + refl->child[dir] = insert_node(refl->child[dir], new); + + rcd = refl->child[dir]; + if ( is_red(rcd) ) { + + if ( is_red(refl->child[!dir]) ) { + + refl->col = RED; + refl->child[0]->col = BLACK; + refl->child[1]->col = BLACK; + + } else { + + if ( is_red(rcd->child[dir] ) ) { + refl = rotate_once(refl, !dir); + } else if ( is_red(rcd->child[!dir] ) ) { + refl = rotate_twice(refl, !dir); + } + + } + } + + } + + return refl; +} + + +/** + * add_refl + * @list: A %RefList + * @h: The 'h' index of the reflection + * @k: The 'k' index of the reflection + * @l: The 'l' index of the reflection + * + * Adds a new reflection to @list. Note that the implementation allows there to + * be multiple reflections with the same indices in the list, so this function + * should succeed even if the given indices already feature in the list. + * + * Returns: The newly created reflection, or NULL on failure. + * + **/ +Reflection *add_refl(RefList *list, signed int h, signed int k, signed int l) +{ + Reflection *new; + Reflection *f; + + assert(abs(h)<256); + assert(abs(k)<256); + assert(abs(l)<256); + + new = new_node(SERIAL(h, k, l)); + if ( new == NULL ) return NULL; + + f = find_refl(list, h, k, l); + if ( f == NULL ) { + + list->head = insert_node(list->head, new); + list->head->col = BLACK; + + } else { + + /* New reflection is identical to a previous one */ + while ( f->next != NULL ) { + f = f->next; + } + f->next = new; + new->prev = f; + } + + return new; +} + + +/** + * add_refl_to_list + * @refl: A %Reflection + * @list: A %RefList + * + * Adds a reflection to @list. The reflection that actually gets added will be + * a newly created one, and all the data will be copied across. The original + * reflection will be destroyed and the new reflection returned. + * + * Returns: The newly created reflection, or NULL on failure. + * + **/ +Reflection *add_refl_to_list(Reflection *refl, RefList *list) +{ + signed int h, k, l; + Reflection *r_added; + + get_indices(refl, &h, &k, &l); + r_added = add_refl(list, h, k, l); + if ( r_added == NULL ) return NULL; + + copy_data(r_added, refl); + reflection_free(refl); + + return r_added; +} + + +/********************************* Iteration **********************************/ + +struct _reflistiterator { + + int stack_size; + int stack_ptr; + Reflection **stack; + +}; + + +/** + * first_refl: + * @list: A %RefList to iterate over + * @piter: Address at which to store a %RefListIterator + * + * This function sets up the state required for iteration over the entire list, + * and then returns the first reflection in the list. An iterator object will + * be created and its address stored at the location given in piter. + * + * Returns: the first reflection in the list. + * + **/ +Reflection *first_refl(RefList *list, RefListIterator **piter) +{ + RefListIterator *iter; + + iter = malloc(sizeof(struct _reflistiterator)); + iter->stack_size = 32; + iter->stack = malloc(iter->stack_size*sizeof(Reflection *)); + iter->stack_ptr = 0; + *piter = iter; + + Reflection *refl = list->head; + + do { + + if ( refl != NULL ) { + iter->stack[iter->stack_ptr++] = refl; + if ( iter->stack_ptr == iter->stack_size ) { + iter->stack_size += 32; + iter->stack = realloc(iter->stack, + iter->stack_size*sizeof(Reflection *)); + } + refl = refl->child[0]; + continue; + } + + if ( iter->stack_ptr == 0 ) { + free(iter->stack); + free(iter); + return NULL; + } + + refl = iter->stack[--iter->stack_ptr]; + + return refl; + + } while ( 1 ); +} + + +/** + * next_refl: + * @refl: A reflection + * @iter: A %RefListIterator + * + * This function looks up the next reflection in the list that was given earlier + * to first_refl(). + * + * Returns: the next reflection in the list, or NULL if no more. + * + **/ +Reflection *next_refl(Reflection *refl, RefListIterator *iter) +{ + int returned = 1; + + do { + + if ( returned ) refl = refl->child[1]; + returned = 0; + + if ( refl != NULL ) { + + iter->stack[iter->stack_ptr++] = refl; + if ( iter->stack_ptr == iter->stack_size ) { + iter->stack_size += 32; + iter->stack = realloc(iter->stack, + iter->stack_size*sizeof(Reflection *)); + } + refl = refl->child[0]; + continue; + + } + if ( iter->stack_ptr == 0 ) { + free(iter->stack); + free(iter); + return NULL; + } + + return iter->stack[--iter->stack_ptr]; + + } while ( 1 ); +} + + +/*********************************** Voodoo ***********************************/ + +static int recursive_depth(Reflection *refl) +{ + int depth_left, depth_right; + + if ( refl == NULL ) return 0; + + depth_left = recursive_depth(refl->child[0]); + depth_right = recursive_depth(refl->child[1]); + + return 1 + biggest(depth_left, depth_right); +} + + +static int recursive_count(Reflection *refl) +{ + int count_left, count_right; + + if ( refl == NULL ) return 0; + + count_left = recursive_count(refl->child[0]); + count_right = recursive_count(refl->child[1]); + + return 1 + count_left + count_right; +} + + +/** + * num_reflections: + * @list: A %RefList + * + * Returns: the number of reflections in @list. + * + **/ +int num_reflections(RefList *list) +{ + return recursive_count(list->head); +} + + +/** + * tree_depth: + * @list: A %RefList + * + * If the depth of the tree is more than about 20, access to the list will be + * slow. This should never happen. + * + * Returns: the depth of the RB-tree used internally to represent @list. + * + **/ +int tree_depth(RefList *list) +{ + return recursive_depth(list->head); +} + + +/** + * lock_reflection: + * @refl: A %Reflection + * + * Acquires a lock on the reflection. + */ +void lock_reflection(Reflection *refl) +{ + pthread_mutex_lock(&refl->lock); +} + + +/** + * unlock_reflection: + * @refl: A %Reflection + * + * Releases a lock on the reflection. + */ +void unlock_reflection(Reflection *refl) +{ + pthread_mutex_unlock(&refl->lock); +} |