aboutsummaryrefslogtreecommitdiff
path: root/src/basis.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/basis.c')
-rw-r--r--src/basis.c234
1 files changed, 1 insertions, 233 deletions
diff --git a/src/basis.c b/src/basis.c
index 8a282e3..841ccc6 100644
--- a/src/basis.c
+++ b/src/basis.c
@@ -1,7 +1,7 @@
/*
* basis.c
*
- * Find approximate lattices to feed various procedures
+ * Handle basis structures
*
* (c) 2007 Thomas White <taw27@cam.ac.uk>
*
@@ -14,10 +14,7 @@
#endif
#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include "utils.h"
#include "reflections.h"
#include "basis.h"
@@ -71,232 +68,3 @@ static double basis_efom(ReflectionList *reflectionlist, Basis *basis) {
}
-static int basis_lfom(ControlContext *ctx, double vx, double vy, double vz) {
-
- Reflection *tcentre;
- int lfom;
- double tol;
- int j;
-
- lfom = 0;
- tol = modulus(vx, vy, vz)/10.0;
- tcentre = ctx->reflectionlist->reflections;
- do {
-
- for ( j=-20; j<=20; j++ ) {
-
- Reflection *check;
-
- check = reflectionlist_find_nearest(ctx->reflectionlist, tcentre->x+vx*j, tcentre->y+vy*j, tcentre->z+vz*j);
- if ( check && (distance3d(check->x, check->y, check->z, tcentre->x+vx*j, tcentre->y+vy*j, tcentre->z+vz*j) < tol) ) {
- lfom++;
- }
-
- }
-
- tcentre = tcentre->next;
-
- } while ( tcentre );
-
- return lfom;
-
-}
-
-/* Select a suitable number of sensible seed vectors */
-static ReflectionList *basis_find_seeds(ControlContext *ctx) {
-
- double tilt_min;
- double tilt_max;
- double tilt_mid;
- ImageRecord *imagerecord;
- double x_temp, y_temp, z_temp;
- double scale;
- double x, y, z;
- Reflection *centre;
- int i;
- ReflectionList *seeds;
-
- seeds = reflectionlist_new();
-
- /* Locate the 'plane' in the middle of the "wedge".
- * This whole procedure assumes there is just one tilt axis. */
- tilt_min = control_min_tilt(ctx);
- tilt_max = control_max_tilt(ctx);
- tilt_mid = tilt_min + (tilt_max-tilt_min)/2;
- imagerecord = control_image_nearest_tilt(ctx, tilt_mid);
-
- /* Apply the last two steps of the mapping transform to get the direction from the origin
- * towards the middle of the wedge */
- x_temp = 0.0;
- y_temp = cos(deg2rad(imagerecord->tilt));
- z_temp = -sin(deg2rad(imagerecord->tilt));
- x = x_temp*cos(-deg2rad(imagerecord->omega)) + y_temp*sin(-deg2rad(imagerecord->omega));
- y = -x_temp*sin(-deg2rad(imagerecord->omega)) + y_temp*cos(-deg2rad(imagerecord->omega));
- z = z_temp;
-
- /* Find the point in the middle of the "wedge" */
- scale = reflectionlist_largest_g(ctx->reflectionlist)/6;
- x *= scale;
- y *= scale;
- z *= scale;
- reflection_add(ctx->reflectionlist, x, y, z, 1.0, REFLECTION_VECTOR_MARKER_2);
-
- /* Find an "origin" reflection */
- centre = reflectionlist_find_nearest(ctx->reflectionlist, x, y, z);
- if ( !centre ) return NULL;
- centre->found = 1;
- reflection_add(ctx->reflectionlist, centre->x, centre->y, centre->z, 1.0, REFLECTION_GENERATED);
-
- for ( i=1; i<=10; i++ ) {
-
- Reflection *vector;
- Reflection *new_seed;
- int accept, lfom;
- double vx, vy, vz;
-
- do {
-
- Reflection *check;
-
- accept = 1;
-
- /* Find a "candidate vector" reflection */
- vector = reflectionlist_find_nearest_longer_unknown(ctx->reflectionlist, centre->x, centre->y, centre->z, 1e9);
- if ( !vector ) {
- printf("BS: Only found %i seeds\n", i);
- return seeds;
- }
- vector->found = 1;
-
- /* Get vector components (not the coordinates the vector was calculated from!) */
- vx = vector->x - centre->x;
- vy = vector->y - centre->y;
- vz = vector->z - centre->z;
-
- /* Proximity test: don't duplicate seeds */
- check = reflectionlist_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 */
- accept = 0;
- continue;
- }
- }
-
- /* Record lFOM for later analysis */
- lfom = basis_lfom(ctx, vx, vy, vz);
-
- } while ( !accept );
-
- /* 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);
-
- }
-
- return seeds;
-
-}
-
-/* Assemble the most sensible basis from seeds */
-static Basis *basis_assemble_seeds(ReflectionList *seeds) {
-
- Basis *basis;
- Reflection *ref;
- int i, j, b;
- ReflectionList *seeds_sorted;
-
- seeds_sorted = reflectionlist_sort_lfom(seeds);
-
- basis = malloc(10*sizeof(Basis));
-
- for ( b=0; b<10; b++ ) {
-
- 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;
-
- }
- 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);
-
- return basis;
-
-}
-
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-03 05:34:19 +0000