aboutsummaryrefslogtreecommitdiff
path: root/src/basis.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/basis.c')
-rw-r--r--src/basis.c144
1 files changed, 97 insertions, 47 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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 02:09:50 +0000