From 12271165c0948536f9b34603432ade4953c97b4e Mon Sep 17 00:00:00 2001 From: taw27 Date: Fri, 28 Sep 2007 17:11:06 +0000 Subject: 'Satisfactory-ish' basis finding git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@135 bf6ca9ba-c028-0410-8290-897cf20841d1 --- src/displaywindow.c | 4 +- src/ipr.c | 208 ++++++++++++++++++++++++++++++---------------------- src/ipr.h | 2 - src/reflections.c | 57 +++++++++++++- src/reflections.h | 19 +++-- src/reproject.c | 9 +-- src/utils.c | 4 + src/utils.h | 1 + 8 files changed, 195 insertions(+), 109 deletions(-) (limited to 'src') diff --git a/src/displaywindow.c b/src/displaywindow.c index eedad54..6a3adfe 100644 --- a/src/displaywindow.c +++ b/src/displaywindow.c @@ -484,7 +484,6 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) { glEnd(); /* x, y, z pointers */ -#if 0 glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red); DRAW_SHORT_POINTER @@ -499,8 +498,7 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) { glRotatef(-90.0, 0.0, 1.0, 0.0); DRAW_SHORT_POINTER glPopMatrix(); -#endif - + /* Tilt axis */ glPushMatrix(); /* Images rotate clockwise by omega to put tilt axis at +x, 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); diff --git a/src/ipr.h b/src/ipr.h index 3de3724..e511756 100644 --- a/src/ipr.h +++ b/src/ipr.h @@ -23,8 +23,6 @@ typedef struct { double y; double z; - double modulus; /* Convenience */ - } Vector; typedef struct basis_struct { diff --git a/src/reflections.c b/src/reflections.c index 63eaa54..a1f6a2c 100644 --- a/src/reflections.c +++ b/src/reflections.c @@ -102,6 +102,7 @@ void reflection_free(ReflectionContext *reflectionctx) { Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y, double z, double intensity, ReflectionType type) { Reflection *new_reflection; + Reflection *nearest; if ( reflectionctx->list_capped ) return NULL; @@ -110,6 +111,10 @@ Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y, fprintf(stderr, "No further reflections will be stored. Go and fix the peak detection.\n"); reflectionctx->list_capped = 1; } + + nearest = reflection_find_nearest_type(reflectionctx, x, y, z, type); + if ( nearest && distance3d(x, y, z, nearest->x, nearest->y, nearest->z) < 0.1e9 ) return NULL; + reflectionctx->n_reflections++; new_reflection = malloc(sizeof(Reflection)); @@ -119,6 +124,7 @@ Reflection *reflection_add(ReflectionContext *reflectionctx, double x, double y, new_reflection->z = z; new_reflection->intensity = intensity; new_reflection->type = type; + new_reflection->found = 0; reflectionctx->last_reflection->next = new_reflection; reflectionctx->last_reflection = new_reflection; @@ -248,7 +254,56 @@ Reflection *reflection_find_nearest(ReflectionContext *reflectionctx, double x, reflection = reflectionctx->reflections; while ( reflection ) { - if ( reflection->type == REFLECTION_NORMAL ) { + if ( (reflection->type == REFLECTION_NORMAL) && (!reflection->found) ) { + double mod; + mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z); + if ( mod < max ) { + max = mod; + best = reflection; + } + } + reflection = reflection->next; + }; + + if ( best ) best->found = 1; + return best; + + +} + +Reflection *reflection_find_nearest_longer(ReflectionContext *reflectionctx, double x, double y, double z, double min_distance) { + + double max = +INFINITY; + Reflection *reflection; + Reflection *best = NULL; + + reflection = reflectionctx->reflections; + while ( reflection ) { + if ( (reflection->type == REFLECTION_NORMAL) && (!reflection->found) ) { + double mod; + mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z); + if ( (mod < max) && (mod >= min_distance) ) { + max = mod; + best = reflection; + } + } + reflection = reflection->next; + }; + + if ( best ) best->found = 1; + return best; + +} + +Reflection *reflection_find_nearest_type(ReflectionContext *reflectionctx, double x, double y, double z, ReflectionType type) { + + double max = +INFINITY; + Reflection *reflection; + Reflection *best = NULL; + + reflection = reflectionctx->reflections; + while ( reflection ) { + if ( reflection->type == type ) { double mod; mod = modulus(x - reflection->x, y - reflection->y, z - reflection->z); if ( mod < max ) { diff --git a/src/reflections.h b/src/reflections.h index 824fcba..9dd281a 100644 --- a/src/reflections.h +++ b/src/reflections.h @@ -29,16 +29,17 @@ typedef enum { typedef struct reflection_struct { - double x; - double y; - double z; - double intensity; + double x; + double y; + double z; + double intensity; - signed int h; - signed int k; - signed int l; + signed int h; + signed int k; + signed int l; - ReflectionType type; + ReflectionType type; + int found; /* This reflection has been used in the seed-finding process */ struct reflection_struct *next; /* MUST BE LAST in order for caching to work */ @@ -68,6 +69,8 @@ extern void reflection_add_from_reflection(ReflectionContext *rctx, Reflection * extern double reflection_largest_g(ReflectionContext *reflectionctx); extern int reflection_map_to_space(ImageReflection *refl, double *ddx, double *ddy, double *ddz, double *twotheta); extern Reflection *reflection_find_nearest(ReflectionContext *reflectionctx, double x, double y, double z); +extern Reflection *reflection_find_nearest_longer(ReflectionContext *reflectionctx, double x, double y, double z, double min_distance); +extern Reflection *reflection_find_nearest_type(ReflectionContext *reflectionctx, double x, double y, double z, ReflectionType type); #endif /* REFLECTION_H */ diff --git a/src/reproject.c b/src/reproject.c index 0971b25..ed359fc 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -33,12 +33,6 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect tilt = deg2rad(image.tilt); omega = deg2rad(image.omega); - refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection)); - - i = 0; - - reflection = rctx->reflections; - /* Calculate the (normalised) incident electron wavevector * n is rotated "into" the reconstruction, so only one omega step. */ nxt = 0.0; nyt = 0.0; nzt = 1.0; @@ -58,6 +52,9 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect ux = uxt*cos(-omega) + uyt*-sin(omega); uy = -uxt*sin(omega) + uyt*cos(omega); uz = uzt; //reflection_add(ctx->reflectionctx, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2); + refl = malloc(MAX_IMAGE_REFLECTIONS*sizeof(ImageReflection)); + reflection = rctx->reflections; + i = 0; do { double xl, yl, zl; diff --git a/src/utils.c b/src/utils.c index eaaa036..4618b44 100644 --- a/src/utils.c +++ b/src/utils.c @@ -39,6 +39,10 @@ double modulus(double x, double y, double z) { return sqrt(x*x + y*y + z*z); } +double distance3d(double x1, double y1, double z1, double x2, double y2, double z2) { + return modulus(x1-x2, y1-y2, z1-z2); +} + /* Angle between two vectors. Answer in radians */ double angle_between(double x1, double y1, double z1, double x2, double y2, double z2) { diff --git a/src/utils.h b/src/utils.h index 99fc748..1d268e6 100644 --- a/src/utils.h +++ b/src/utils.h @@ -27,6 +27,7 @@ extern double angle_between(double x1, double y1, double z1, double x2, double y extern double angle_between_d(double x1, double y1, double z1, double x2, double y2, double z2); extern double lambda(double voltage); extern void matrix_renormalise(gsl_matrix *m); +extern double distance3d(double x1, double y1, double z1, double x2, double y2, double z2); #define rad2deg(a) ((a)*180/M_PI) #define deg2rad(a) ((a)*M_PI/180) -- cgit v1.2.3