aboutsummaryrefslogtreecommitdiff
path: root/src/ipr.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/ipr.c')
-rw-r--r--src/ipr.c208
1 files changed, 119 insertions, 89 deletions
diff --git a/src/ipr.c b/src/ipr.c
index 47b57e2..dd6ce59 100644
--- a/src/ipr.c
+++ b/src/ipr.c
@@ -27,64 +27,6 @@
#include "ipr.h"
#include "displaywindow.h"
-#if 0
-static int ipr_random_image(ControlContext *ctx) {
- double n;
- n = (double)random()/RAND_MAX;
- n *= ctx->n_images;
- return (int)n;
-}
-#endif
-
-static Basis *ipr_choose_random_basis(ControlContext *ctx) {
-
- Basis *basis;
- 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;
-
- /* 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 = reflection_largest_g(ctx->reflectionctx)/4;
- x *= scale;
- y *= scale;
- z *= scale;
- reflection_add(ctx->reflectionctx, x, y, z, 1.0, REFLECTION_VECTOR_MARKER_2);
-
- centre = reflection_find_nearest(ctx->reflectionctx, x, y, z);
- if ( !centre ) return NULL;
- reflection_add(ctx->reflectionctx, centre->x, centre->y, centre->z, 1.0, REFLECTION_MARKER);
-
- basis = malloc(sizeof(Basis));
-
- basis->a.x = 1.842e9; basis->a.y = 0.0; basis->a.z = 0.0;
- basis->b.x = 0.0; basis->b.y = 1.842e9; basis->b.z = 0.0;
- basis->c.x = 0.0; basis->c.y = 0.0; basis->c.z = 1.842e9;
-
- return basis;
-
-}
-
static double ipr_efom(ReflectionContext *rctx, Basis *basis) {
int n_indexed, n_counted;
@@ -97,7 +39,7 @@ static double ipr_efom(ReflectionContext *rctx, Basis *basis) {
if ( cur->type == REFLECTION_NORMAL ) {
- /* Can this basis approximately account for this reflection? */
+ /* Can this basis "approximately" account for this reflection? */
double det;
double a11, a12, a13, a21, a22, a23, a31, a32, a33;
double h, k, l;
@@ -135,40 +77,126 @@ static double ipr_efom(ReflectionContext *rctx, Basis *basis) {
}
-static Basis *ipr_choose_initial_basis(ControlContext *ctx) {
+static Basis *ipr_find_basis(ControlContext *ctx) {
- Basis *basis;
- Basis *best_basis;
- double best_efom;
- int i;
+ Basis *basis;
+ 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;
- /* Track the best basis throughout the whole procedure */
- best_basis = NULL; best_efom = 0;
+ /* 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);
- printf("IP: Finding initial basis using eFOM...\n");
- for ( i=1; i<=30; i++ ) {
-
- double efom;
+ /* 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 = reflection_largest_g(ctx->reflectionctx)/4;
+ x *= scale;
+ y *= scale;
+ z *= scale;
+ reflection_add(ctx->reflectionctx, x, y, z, 1.0, REFLECTION_VECTOR_MARKER_2);
+
+ /* Find an "origin" reflection */
+ centre = reflection_find_nearest(ctx->reflectionctx, x, y, z);
+ if ( !centre ) return NULL;
+ reflection_add(ctx->reflectionctx, centre->x, centre->y, centre->z, 1.0, REFLECTION_GENERATED);
+
+ basis = malloc(sizeof(Basis));
+ for ( i=1; i<=3; i++ ) {
+
+ Reflection *vector;
+ int accept;
+ double vx, vy, vz;
- basis = ipr_choose_random_basis(ctx);
- efom = ipr_efom(ctx->reflectionctx, basis);
+ do {
+
+ Reflection *check;
+ int orders, j, lfom;
+ double tol;
+
+ accept = 1;
- if ( (efom > best_efom) || !best_basis ) {
- if ( best_basis ) free(best_basis);
- best_efom = efom;
- best_basis = basis;
- printf("IP: %3i: eFOM = %7.3f %%\n", i, efom*100);
+ /* Find a "candidate vector" reflection */
+ vector = reflection_find_nearest_longer(ctx->reflectionctx, centre->x, centre->y, centre->z, 1e9); /* 0.5 A^-1 */
+ if ( !vector ) {
+ printf("IP: Couldn't find enough seeds\n");
+ return NULL;
+ }
+
+ /* 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 */
+ check = reflection_find_nearest_type(ctx->reflectionctx, vx, vy, vz, REFLECTION_MARKER);
+ if ( check ) {
+ if ( distance3d(vx, vy, vz, check->x, check->y, check->z) < 1e9 ) {
+ /* Too close to another seed */
+ accept = 0;
+ continue;
+ }
+ }
+
+ /* Line FoM test */
+ orders = 100e9 / modulus(vx, vy, vz); /* Go out by +/- this number of iterations. Rounding error ignored */
+ lfom = 0;
+ tol = modulus(vx, vy, vz)/10.0;
+ for ( j=-orders; j<=orders; j++ ) {
+ check = reflection_find_nearest(ctx->reflectionctx, centre->x+vx*j, centre->y+vy*j, centre->z+vz*j);
+ if ( check && (distance3d(check->x, check->y, check->z, centre->x+vx*j, centre->y+vy*j, centre->z+vz*j) < tol) ) {
+ lfom++;
+ }
+ }
+ if ( lfom < 1 ) {
+ accept = 0;
+ continue;
+ }
+
+ } while ( !accept );
+ reflection_add(ctx->reflectionctx, vx, vy, vz, 1.0, REFLECTION_MARKER);
+ 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;
+ }
}
}
- printf("IP: Initial basis:\n");
- printf("IP: a = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->a.x, basis->a.y, basis->a.z, basis->a.modulus);
- printf("IP: b = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->b.x, basis->b.y, basis->b.z, basis->b.modulus);
- printf("IP: c = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->c.x, basis->c.y, basis->c.z, basis->c.modulus);
+ printf("IP: eFOM = %7.3f %%\n", ipr_efom(ctx->reflectionctx, basis)*100);
return basis;
+
}
/* Generate list of reflections to analyse */
@@ -191,14 +219,14 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
}
max_res = max_res / 2;
} else {
- max_res = 2e10; //works until I put some more in the reflect.cache header
+ max_res = 1e10; //works until I put some more in the reflect.cache header
}
ordered = reflection_init();
- max_order_a = max_res/basis->a.modulus;
- max_order_b = max_res/basis->b.modulus;
- max_order_c = max_res/basis->c.modulus;
+ 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);
+ max_order_c = max_res/modulus(basis->c.x, basis->c.y, basis->c.z);
for ( h=-max_order_a; h<=max_order_a; h++ ) {
for ( k=-max_order_b; k<=max_order_b; k++ ) {
for ( l=-max_order_c; l<=max_order_c; l++ ) {
@@ -217,6 +245,7 @@ static ReflectionContext *ipr_generate(ControlContext *ctx, Basis *basis) {
}
if ( (h!=0) || (k!=0) || (l!=0) ) {
// reflection_add(ctx->reflectionctx, x, y, z, 1, REFLECTION_GENERATED);
+ reflection_add(ordered, x, y, z, 1, REFLECTION_GENERATED);
}
}
@@ -256,10 +285,11 @@ int ipr_refine(ControlContext *ctx) {
Basis *basis;
int finished;
- srand(time(NULL));
-
- basis = ipr_choose_initial_basis(ctx);
- if ( !basis ) return -1;
+ basis = ipr_find_basis(ctx);
+ if ( !basis ) {
+ printf("IP: Unable to find basis\n");
+ return -1;
+ }
ctx->ipr_basis = basis;
ctx->ipr_lat = ipr_generate(ctx, basis);
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-02 20:58:18 +0000