diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-09-28 17:11:06 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-09-28 17:11:06 +0000 |
commit | 12271165c0948536f9b34603432ade4953c97b4e (patch) | |
tree | 8d52c7685ca99b7b836babc05b940ab537c05cf3 /src/ipr.c | |
parent | 1ed23746f6d27f648e3a5f96bf499823069fd171 (diff) |
'Satisfactory-ish' basis finding
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@135 bf6ca9ba-c028-0410-8290-897cf20841d1
Diffstat (limited to 'src/ipr.c')
-rw-r--r-- | src/ipr.c | 208 |
1 files changed, 119 insertions, 89 deletions
@@ -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); |