diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/basis.c | 144 | ||||
-rw-r--r-- | src/cache.c | 2 | ||||
-rw-r--r-- | src/ipr.c | 2 | ||||
-rw-r--r-- | src/mapping.c | 2 | ||||
-rw-r--r-- | src/qdrp.c | 2 | ||||
-rw-r--r-- | src/reflections.c | 87 | ||||
-rw-r--r-- | src/reflections.h | 8 |
7 files changed, 168 insertions, 79 deletions
diff --git a/src/basis.c b/src/basis.c index 939dc93..19004ee 100644 --- a/src/basis.c +++ b/src/basis.c @@ -102,6 +102,7 @@ static int basis_lfom(ControlContext *ctx, double vx, double vy, double vz) { } +/* Select a suitable number of sensible seed vectors */ static ReflectionList *basis_find_seeds(ControlContext *ctx) { double tilt_min; @@ -115,7 +116,7 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) { int i; ReflectionList *seeds; - seeds = reflection_init(); + seeds = reflectionlist_new(); /* Locate the 'plane' in the middle of the "wedge". * This whole procedure assumes there is just one tilt axis. */ @@ -134,7 +135,7 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) { z = z_temp; /* Find the point in the middle of the "wedge" */ - scale = reflection_largest_g(ctx->reflectionlist)/4; + scale = reflection_largest_g(ctx->reflectionlist)/6; x *= scale; y *= scale; z *= scale; @@ -149,21 +150,21 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) { for ( i=1; i<=10; i++ ) { Reflection *vector; - int accept; + Reflection *new_seed; + int accept, lfom; double vx, vy, vz; do { Reflection *check; - int lfom; accept = 1; /* Find a "candidate vector" reflection */ - vector = reflection_find_nearest_longer(ctx->reflectionlist, centre->x, centre->y, centre->z, 1e9); /* 0.5 A^-1 */ + vector = reflection_find_nearest_longer_unknown(ctx->reflectionlist, centre->x, centre->y, centre->z, 1e9); /* 0.5 A^-1 */ if ( !vector ) { - printf("BS: Couldn't find enough seeds\n"); - return NULL; + printf("BS: Only found %i seeds\n", i); + return seeds; } vector->found = 1; @@ -172,8 +173,8 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) { vy = vector->y - centre->y; vz = vector->z - centre->z; - /* Proximity test */ - check = reflection_find_nearest_type(ctx->reflectionlist, vx, vy, vz, REFLECTION_NORMAL); + /* Proximity test: don't duplicate seeds */ + check = reflection_find_nearest(seeds, vx, vy, vz); if ( check ) { if ( distance3d(vx, vy, vz, check->x, check->y, check->z) < 1e9 ) { /* Too close to another seed */ @@ -182,17 +183,16 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) { } } - /* lFOM test */ + /* Record lFOM for later analysis */ lfom = basis_lfom(ctx, vx, vy, vz); - if ( lfom < 1 ) { - accept = 0; - continue; - } - printf("lfom=%i\n", lfom); } while ( !accept ); - reflection_add(seeds, vx, vy, vz, 1.0, REFLECTION_NORMAL); + /* Add to the list of seeds */ + new_seed = reflection_add(seeds, vx, vy, vz, 1.0, REFLECTION_NORMAL); + new_seed->lfom = lfom; + + /* Create a marker in the default list for visualisation */ reflection_add(ctx->reflectionlist, vx, vy, vz, 1.0, REFLECTION_MARKER); } @@ -201,48 +201,98 @@ static ReflectionList *basis_find_seeds(ControlContext *ctx) { } -Basis *basis_find(ControlContext *ctx) { +/* Assemble the most sensible basis from seeds */ +static Basis *basis_assemble_seeds(ReflectionList *seeds) { Basis *basis; - ReflectionList *seeds; Reflection *ref; - int i; - - /* Get the shortlist of seeds */ - seeds = basis_find_seeds(ctx); + int i, j, b; + ReflectionList *seeds_sorted; - /* Assemble the seeds into a basis */ - basis = malloc(sizeof(Basis)); - ref = seeds->reflections; - for ( i=1; i<=3; i++ ) { + seeds_sorted = reflection_sort_lfom(seeds); - double vx, vy, vz; + basis = malloc(10*sizeof(Basis)); - vx = ref->x; - vy = ref->y; - vz = ref->z; + for ( b=0; b<10; b++ ) { - switch ( i ) { - case 1 : { - basis->a.x = vx; - basis->a.y = vy; - basis->a.z = vz; - } - case 2 : { - basis->b.x = vx; - basis->b.y = vy; - basis->b.z = vz; - } - case 3 : { - basis->c.x = vx; - basis->c.y = vy; - basis->c.z = vz; + ref = seeds_sorted->reflections; + j = 0; /* Number of basis components already found */ + for ( i=1; i<=seeds_sorted->n_reflections; i++ ) { + + double vx, vy, vz; + + vx = ref->x; + vy = ref->y; + vz = ref->z; + + printf("Seed %2i: lFOM=%6i", i, ref->lfom); + + switch ( j ) { + case 0 : { + basis[b].a.x = vx; + basis[b].a.y = vy; + basis[b].a.z = vz; + j++; + printf(" *"); + break; + } + case 1 : { + if ( (angle_between(vx, vy, vz, basis[b].a.x, basis[b].a.y, basis[b].a.z) > M_PI/6) + && (angle_between(vx, vy, vz, basis[b].a.x, basis[b].a.y, basis[b].a.z) < M_PI-M_PI/6) ) { + basis[b].b.x = vx; + basis[b].b.y = vy; + basis[b].b.z = vz; + j++; + printf(" * (%4.1f deg)", rad2deg(angle_between(vx, vy, vz, basis[b].a.x, basis[b].a.y, basis[b].a.z))); + break; + } + } + case 2 : { + double cx, cy, cz; + cx = basis[b].a.y*basis[b].b.z - basis[b].a.z*basis[b].b.y; + cy = - basis[b].a.x*basis[b].b.z + basis[b].a.z*basis[b].b.x; + cz = basis[b].a.x*basis[b].b.y - basis[b].a.y*basis[b].b.x; + if ( (angle_between(vx, vy, vz, basis[b].a.x, basis[b].a.y, basis[b].a.z) > M_PI/6) + && (angle_between(vx, vy, vz, basis[b].b.x, basis[b].b.y, basis[b].b.z) > M_PI/6) + && (angle_between(vx, vy, vz, cx, cy, cz) < M_PI/2-M_PI/6) ) { + basis[b].c.x = vx; + basis[b].c.y = vy; + basis[b].c.z = vz; + j++; + printf(" * (%4.1f deg)", rad2deg(angle_between(vx, vy, vz, cx, cy, cz))); + break; + } + } } + + printf("\n"); + if ( j >= 3 ) break; + ref = ref->next; + } - - ref = ref->next; + if ( j < 3 ) { + break; + } + } + return basis; + +} + +Basis *basis_find(ControlContext *ctx) { + + Basis *basis; + ReflectionList *seeds; + + /* Get the shortlist of seeds */ + seeds = basis_find_seeds(ctx); + if ( seeds->n_reflections < 3 ) { + printf("BS: Not enough seeds to form a basis\n"); + return NULL; } + printf("BS: Found %i seeds\n", seeds->n_reflections); + + basis = basis_assemble_seeds(seeds); printf("BS: eFOM = %7.3f %%\n", basis_efom(ctx->reflectionlist, basis)*100); diff --git a/src/cache.c b/src/cache.c index eb164aa..078fa72 100644 --- a/src/cache.c +++ b/src/cache.c @@ -38,7 +38,7 @@ ReflectionList *cache_load(const char *filename) { cachedreflection_size = sizeof(Reflection) - sizeof(Reflection *); - reflectionlist = reflection_init(); + reflectionlist = reflectionlist_new(); f = fopen(filename, "rb"); if ( !f ) { fprintf(stderr, "Couldn't read reflections from cache\n"); @@ -51,7 +51,7 @@ static ReflectionList *ipr_generate(ControlContext *ctx, Basis *basis) { max_res = 1e10; //works until I put some more in the reflect.cache header } - ordered = reflection_init(); + ordered = reflectionlist_new(); max_order_a = max_res/modulus(basis->a.x, basis->a.y, basis->a.z); max_order_b = max_res/modulus(basis->b.x, basis->b.y, basis->b.z); diff --git a/src/mapping.c b/src/mapping.c index 8570966..c2fed73 100644 --- a/src/mapping.c +++ b/src/mapping.c @@ -23,7 +23,7 @@ ReflectionList *mapping_create(ControlContext *ctx) { int i; /* Create reflection context */ - reflectionlist = reflection_init(); + reflectionlist = reflectionlist_new(); /* Pass all images through itrans * (let itrans add the reflections to reflectionlist for now) */ @@ -107,7 +107,7 @@ static int qdrp_parseline(ControlContext *ctx, const char *line) { } ctx->started = 1; - ctx->reflectionlist = reflection_init(); + ctx->reflectionlist = reflectionlist_new(); } diff --git a/src/reflections.c b/src/reflections.c index a0ff76d..56871f0 100644 --- a/src/reflections.c +++ b/src/reflections.c @@ -22,26 +22,18 @@ #include "reflections.h" #include "utils.h" -static void reflection_addfirst(ReflectionList *reflectionlist) { - /* Create first items on lists - saves faffing later. Corresponds to a central marker. - Reflections are only stored if they have non-zero value. */ - reflectionlist->reflections = malloc(sizeof(Reflection)); - reflectionlist->reflections->next = NULL; - reflectionlist->reflections->x = 0; - reflectionlist->reflections->y = 0; - reflectionlist->reflections->z = 0; - reflectionlist->reflections->type = REFLECTION_CENTRAL; - reflectionlist->last_reflection = reflectionlist->reflections; - reflectionlist->n_reflections = 1; +static void reflectionlist_init(ReflectionList *reflectionlist) { + reflectionlist->n_reflections = 0; reflectionlist->list_capped = 0; + reflectionlist->reflections = NULL; + reflectionlist->last_reflection = NULL; } -ReflectionList *reflection_init() { +ReflectionList *reflectionlist_new() { ReflectionList *reflectionlist = malloc(sizeof(ReflectionList)); - reflectionlist->n_reflections = 0; - reflection_addfirst(reflectionlist); - reflectionlist->list_capped = 0; + + reflectionlist_init(reflectionlist); return reflectionlist; @@ -87,15 +79,12 @@ void reflection_clear(ReflectionList *reflectionlist) { reflection = next; } while ( reflection ); - reflectionlist->n_reflections = 0; - reflectionlist->list_capped = 0; - reflection_addfirst(reflectionlist); + reflectionlist_init(reflectionlist); } void reflection_free(ReflectionList *reflectionlist) { reflection_clear(reflectionlist); - free(reflectionlist->reflections); free(reflectionlist); } @@ -115,8 +104,6 @@ Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, d nearest = reflection_find_nearest_type(reflectionlist, x, y, z, type); if ( nearest && distance3d(x, y, z, nearest->x, nearest->y, nearest->z) < 0.1e9 ) return NULL; - reflectionlist->n_reflections++; - new_reflection = malloc(sizeof(Reflection)); new_reflection->next = NULL; new_reflection->x = x; @@ -126,8 +113,14 @@ Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, d new_reflection->type = type; new_reflection->found = 0; - reflectionlist->last_reflection->next = new_reflection; - reflectionlist->last_reflection = new_reflection; + if ( reflectionlist->last_reflection ) { + reflectionlist->last_reflection->next = new_reflection; + reflectionlist->last_reflection = new_reflection; + } else { + reflectionlist->reflections = new_reflection; + reflectionlist->last_reflection = new_reflection; + } + reflectionlist->n_reflections++; return new_reflection; @@ -222,9 +215,15 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord } void reflection_add_from_reflection(ReflectionList *reflectionlist, Reflection *r) { + if ( reflectionlist->last_reflection ) { + reflectionlist->last_reflection->next = r; + reflectionlist->last_reflection = r; + } else { + reflectionlist->reflections = r; + reflectionlist->last_reflection = r; + } r->next = NULL; - reflectionlist->last_reflection->next = r; - reflectionlist->last_reflection = r; + reflectionlist->n_reflections++; } double reflection_largest_g(ReflectionList *reflectionlist) { @@ -269,7 +268,7 @@ Reflection *reflection_find_nearest(ReflectionList *reflectionlist, double x, do } -Reflection *reflection_find_nearest_longer(ReflectionList *reflectionlist, double x, double y, double z, double min_distance) { +Reflection *reflection_find_nearest_longer_unknown(ReflectionList *reflectionlist, double x, double y, double z, double min_distance) { double max = +INFINITY; Reflection *reflection; @@ -316,3 +315,39 @@ Reflection *reflection_find_nearest_type(ReflectionList *reflectionlist, double } +/* This destroys the lfom values in the input list */ +ReflectionList *reflection_sort_lfom(ReflectionList *in) { + + ReflectionList *out; + Reflection *best; + + out = reflectionlist_new(); + + do { + + Reflection *reflection; + int lfom = 0; + + reflection = in->reflections; + best = NULL; + while ( reflection ) { + if ( reflection->lfom > lfom ) { + best = reflection; + lfom = reflection->lfom; + } + reflection = reflection->next; + }; + + if ( best ) { + Reflection *new; + new = reflection_add(out, best->x, best->y, best->z, best->intensity, best->type); + new->lfom = best->lfom; + best->lfom = 0; + } + + } while ( best ); + + return out; + +} + diff --git a/src/reflections.h b/src/reflections.h index 9ac0cec..c753430 100644 --- a/src/reflections.h +++ b/src/reflections.h @@ -39,7 +39,10 @@ typedef struct reflection_struct { signed int l; ReflectionType type; + + /* Stuff used when finding bases */ int found; /* This reflection has been used in the seed-finding process */ + int lfom; /* Line FoM for this seed */ struct reflection_struct *next; /* MUST BE LAST in order for caching to work */ @@ -55,7 +58,7 @@ typedef struct reflectionlist_struct { } ReflectionList; -extern ReflectionList *reflection_init(void); +extern ReflectionList *reflectionlist_new(void); extern void reflection_clear(ReflectionList *reflectionlist); extern void reflection_clear_markers(ReflectionList *reflectionlist); extern void reflection_free(ReflectionList *reflectionlist); @@ -69,8 +72,9 @@ extern void reflection_add_from_reflection(ReflectionList *reflectionlist, Refle extern double reflection_largest_g(ReflectionList *reflectionlist); extern int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta); extern Reflection *reflection_find_nearest(ReflectionList *reflectionlist, double x, double y, double z); -extern Reflection *reflection_find_nearest_longer(ReflectionList *reflectionlist, double x, double y, double z, double min_distance); +extern Reflection *reflection_find_nearest_longer_unknown(ReflectionList *reflectionlist, double x, double y, double z, double min_distance); extern Reflection *reflection_find_nearest_type(ReflectionList *reflectionlist, double x, double y, double z, ReflectionType type); +extern ReflectionList *reflection_sort_lfom(ReflectionList *in); #endif /* REFLECTION_H */ |