From df1db2ce74a89f7270122e3edcfea73b07485410 Mon Sep 17 00:00:00 2001 From: taw27 Date: Fri, 23 Nov 2007 13:34:46 +0000 Subject: Break off cell delta routine git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@210 bf6ca9ba-c028-0410-8290-897cf20841d1 --- src/basis.c | 12 +++ src/basis.h | 1 + src/refine.c | 314 ++++++++++++++++++++++++++++++++++++++++++----------------- 3 files changed, 236 insertions(+), 91 deletions(-) diff --git a/src/basis.c b/src/basis.c index 08e3edb..f3b9b7a 100644 --- a/src/basis.c +++ b/src/basis.c @@ -68,3 +68,15 @@ double basis_efom(ReflectionList *reflectionlist, Basis *basis) { } +Basis basis_add(Basis u, Basis v) { + + Basis ans; + + ans.a.x = u.a.x + v.a.x; ans.a.y = u.a.y + v.a.y; ans.a.z = u.a.z + v.a.z; + ans.b.x = u.b.x + v.b.x; ans.b.y = u.b.y + v.b.y; ans.b.z = u.b.z + v.b.z; + ans.c.x = u.c.x + v.c.x; ans.c.y = u.c.y + v.c.y; ans.c.z = u.c.z + v.c.z; + + return ans; + +} + diff --git a/src/basis.h b/src/basis.h index 5803cab..f8248c6 100644 --- a/src/basis.h +++ b/src/basis.h @@ -33,6 +33,7 @@ typedef struct basis_struct { } Basis; extern double basis_efom(struct reflectionlist_struct *reflectionlist, Basis *basis); +extern Basis basis_add(Basis u, Basis v); #endif /* BASIS_H */ diff --git a/src/refine.c b/src/refine.c index 12f0851..4ae7ec8 100644 --- a/src/refine.c +++ b/src/refine.c @@ -57,12 +57,113 @@ static const char *refine_decode(RefinementIndex i) { } #endif +/* Return the "cell delta" to add to "cell" in order to make "feature" fit its partner */ +static Basis refine_cell_delta(Basis *cell, ImageFeature *feature, RefinementIndex *sharedr) { + + double dix, diy; + double dx, dy, dz, twotheta; + double old_x, old_y; + RefinementIndex index, delta, shared; + signed int h, k, l; + double a11, a12, a13, a21, a22, a23, a31, a32, a33, det; + double dh, dk, dl; + Basis cd; + + h = feature->reflection->h; + k = feature->reflection->k; + l = feature->reflection->l; + + cd.a.x = 0.0; cd.a.y = 0.0; cd.a.z = 0.0; + cd.b.x = 0.0; cd.b.y = 0.0; cd.b.z = 0.0; + cd.c.x = 0.0; cd.c.y = 0.0; cd.c.z = 0.0; + + /* Determine the difference vector */ + dix = feature->partner->x - feature->x; + diy = feature->partner->y - feature->y; +// printf("RF: Feature %3i: %3i %3i %3i dev = %f %f px\n", i, h, k, l, dix, diy); + + /* Map the difference vector to the relevant tilted plane */ + old_x = feature->partner->x; + old_y = feature->partner->y; + feature->partner->x = dix + feature->partner->parent->x_centre; + feature->partner->y = diy + feature->partner->parent->y_centre; + mapping_map_to_space(feature->partner, &dx, &dy, &dz, &twotheta); + feature->partner->x = old_x; + feature->partner->y = old_y; +// printf("RF: dev=%8e %8e %8e (%5f mrad)\n", dx, dy, dz, twotheta*1e3); + + /* Select the basis vectors which are allowed to be altered */ + index = 0; + if ( h ) index |= INDEX_A; + if ( k ) index |= INDEX_B; + if ( l ) index |= INDEX_C; + assert(index != 0); /* Can't refine using the central beam! */ + + /* Set up the coordinate transform from hkl to xyz */ + a11 = cell->a.x; a12 = cell->a.y; a13 = cell->a.z; + a21 = cell->b.x; a22 = cell->b.y; a23 = cell->b.z; + a31 = cell->c.x; a32 = cell->c.y; a33 = cell->c.z; + + /* Invert the matrix to get dh,dk,dl from dx,dy,dz */ + det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31); + dh = ((a22*a33-a23*a32)*dx + (a23*a31-a21*a33)*dy + (a21*a32-a22*a31)*dz) / det; + dk = ((a13*a32-a12*a33)*dx + (a11*a33-a13*a31)*dy + (a12*a31-a11*a32)*dz) / det; + dl = ((a12*a23-a13*a22)*dx + (a13*a21-a11*a23)*dy + (a11*a22-a12*a21)*dz) / det; +// printf("RF: dev(hkl) = %f %f %f\n", dh, dk, dl); + + delta = 0; + if ( fabs(dh)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_A; + if ( fabs(dk)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_B; + if ( fabs(dl)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_C; + + shared = index & delta; + +// printf("RF: index: %s\nRF: delta: %s\nRF: shared: %s\n", refine_decode(index), refine_decode(delta), refine_decode(shared)); + + if ( shared == 0 ) { + /* No indices left - 'pure shear' (delta is perpendicular (in the abc basis) to index) */ + shared = index; + // printf("RF: Pure shear.\n"); + } + + if ( shared & INDEX_A ) { + double w = (double)abs(h) / (abs(h)+abs(k)+abs(l)); + // printf("RF: w(a) = %f\n", w); + cd.a.x = w*dx / (double)h; + cd.a.y = w*dy / (double)h; + cd.a.z = w*dz / (double)h; + } + if ( shared & INDEX_B ) { + double w = (double)abs(k) / (abs(h)+abs(k)+abs(l)); + // printf("RF: w(b) = %f\n", w); + cd.b.x = w*dx / (double)k; + cd.b.y = w*dy / (double)k; + cd.b.z = w*dz / (double)k; + } + if ( shared & INDEX_C ) { + double w = (double)abs(l) / (abs(h)+abs(k)+abs(l)); + // printf("RF: w(c) = %f\n", w); + cd.c.x = w*dx / (double)l; + cd.c.y = w*dy / (double)l; + cd.c.z = w*dz / (double)l; + } + +// printf("RF: Distortion along x: %+8e = %+8e\n", h*cd.a.x + k*cd.b.x + l*cd.c.x, dx); +// printf("RF: Distortion along y: %+8e = %+8e\n", h*cd.a.y + k*cd.b.y + l*cd.c.y, dy); +// printf("RF: Distortion along z: %+8e = %+8e\n", h*cd.a.z + k*cd.b.z + l*cd.c.z, dz); + + *sharedr = shared; + + return cd; + +} + /* Use the IPR algorithm to make "cell" fit the given image */ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, ReflectionList *cell_lattice) { ImageFeatureList *flist; int i; - Basis cd; /* Cell delta */ + Basis cd; /* Overall cell delta */ int n_a = 0; int n_b = 0; int n_c = 0; @@ -79,98 +180,17 @@ static ImageFeature *refine_fit_image(Basis *cell, ImageRecord *image, Reflectio for ( i=0; irflist->n_features; i++ ) { - double dix, diy; - double dx, dy, dz, twotheta; - double old_x, old_y; - RefinementIndex index, delta, shared; - signed int h, k, l; - double a11, a12, a13, a21, a22, a23, a31, a32, a33, det; - double dh, dk, dl; + Basis this_cd; + RefinementIndex shared; /* Skip if no partner */ if ( !image->rflist->features[i].partner ) continue; - h = image->rflist->features[i].reflection->h; - k = image->rflist->features[i].reflection->k; - l = image->rflist->features[i].reflection->l; - - /* Determine the difference vector */ - dix = image->rflist->features[i].partner->x - image->rflist->features[i].x; - diy = image->rflist->features[i].partner->y - image->rflist->features[i].y; - // printf("RF: Feature %3i: %3i %3i %3i dev = %f %f px\n", i, h, k, l, dix, diy); - - /* Map the difference vector to the relevant tilted plane */ - old_x = image->rflist->features[i].partner->x; - old_y = image->rflist->features[i].partner->y; - image->rflist->features[i].partner->x = dix + image->rflist->features[i].partner->parent->x_centre; - image->rflist->features[i].partner->y = diy + image->rflist->features[i].partner->parent->y_centre; - mapping_map_to_space(image->rflist->features[i].partner, &dx, &dy, &dz, &twotheta); - image->rflist->features[i].partner->x = old_x; - image->rflist->features[i].partner->y = old_y; - // printf("RF: dev=%8e %8e %8e (%5f mrad)\n", dx, dy, dz, twotheta*1e3); - - /* Select the basis vectors which are allowed to be altered */ - index = 0; - if ( h ) index |= INDEX_A; - if ( k ) index |= INDEX_B; - if ( l ) index |= INDEX_C; - assert(index != 0); /* Can't refine using the central beam! */ - - /* Set up the coordinate transform from hkl to xyz */ - a11 = cell->a.x; a12 = cell->a.y; a13 = cell->a.z; - a21 = cell->b.x; a22 = cell->b.y; a23 = cell->b.z; - a31 = cell->c.x; a32 = cell->c.y; a33 = cell->c.z; - - /* Invert the matrix to get dh,dk,dl from dx,dy,dz */ - det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31); - dh = ((a22*a33-a23*a32)*dx + (a23*a31-a21*a33)*dy + (a21*a32-a22*a31)*dz) / det; - dk = ((a13*a32-a12*a33)*dx + (a11*a33-a13*a31)*dy + (a12*a31-a11*a32)*dz) / det; - dl = ((a12*a23-a13*a22)*dx + (a13*a21-a11*a23)*dy + (a11*a22-a12*a21)*dz) / det; - // printf("RF: dev(hkl) = %f %f %f\n", dh, dk, dl); - - delta = 0; - if ( fabs(dh)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_A; - if ( fabs(dk)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_B; - if ( fabs(dl)/(fabs(dh)+fabs(dk)+fabs(dl)) > 0.001 ) delta |= INDEX_C; - - shared = index & delta; - - // printf("RF: index: %s\nRF: delta: %s\nRF: shared: %s\n", refine_decode(index), refine_decode(delta), refine_decode(shared)); - - if ( shared == 0 ) { - /* No indices left - 'pure shear' (delta is perpendicular (in the abc basis) to index) */ - shared = index; - // printf("RF: Pure shear.\n"); - } - - if ( shared & INDEX_A ) { - double w = (double)abs(h) / (abs(h)+abs(k)+abs(l)); - // printf("RF: w(a) = %f\n", w); - cd.a.x += w*dx / (double)h; - cd.a.y += w*dy / (double)h; - cd.a.z += w*dz / (double)h; - n_a++; - } - if ( shared & INDEX_B ) { - double w = (double)abs(k) / (abs(h)+abs(k)+abs(l)); - // printf("RF: w(b) = %f\n", w); - cd.b.x += w*dx / (double)k; - cd.b.y += w*dy / (double)k; - cd.b.z += w*dz / (double)k; - n_b++; - } - if ( shared & INDEX_C ) { - double w = (double)abs(l) / (abs(h)+abs(k)+abs(l)); - // printf("RF: w(c) = %f\n", w); - cd.c.x += w*dx / (double)l; - cd.c.y += w*dy / (double)l; - cd.c.z += w*dz / (double)l; - n_c++; - } - - // printf("RF: Distortion along x: %+8e = %+8e\n", h*cd.a.x + k*cd.b.x + l*cd.c.x, dx); - // printf("RF: Distortion along y: %+8e = %+8e\n", h*cd.a.y + k*cd.b.y + l*cd.c.y, dy); - // printf("RF: Distortion along z: %+8e = %+8e\n", h*cd.a.z + k*cd.b.z + l*cd.c.z, dz); + this_cd = refine_cell_delta(cell, &image->rflist->features[i], &shared); + cd = basis_add(cd, this_cd); + if ( shared & INDEX_A ) n_a++; + if ( shared & INDEX_B ) n_b++; + if ( shared & INDEX_C ) n_c++; done = TRUE; @@ -232,6 +252,112 @@ static gint refine_step(GtkWidget *step_button, ControlContext *ctx) { } +static int refine_sequence_random(ControlContext *ctx, double *fit, double *warp) { + + int i; + + for ( i=0; i<1000; i++ ) { + + int ni; + int nf; + ImageRecord *image; + ImageFeature *feature; + + /* Ensure lattice is up to date */ + reproject_lattice_changed(ctx); + ctx->images->images[i].rflist = NULL; /* Invalidate reprojection for this image - it's wrong */ + + /* Select a random feature from a random image */ + ni = ((float)ctx->images->n_images * random()) / RAND_MAX; + image = &ctx->images->images[ni]; + nf = ((float)image->features->n_features * random()) / RAND_MAX; + feature = &image->features->features[nf]; + + + + } + + *fit = 0.0; + *warp = 0.0; + + return 0; + +} + +static gint refine_random_sequence(GtkWidget *widget, ControlContext *ctx) { + + double omega_offs; + GtkWidget *window_fit; + GtkWidget *graph_fit; + double *fit_vals; + GtkWidget *window_warp; + GtkWidget *graph_warp; + double *warp_vals; + size_t idx; + ImageDisplay *id; + + fit_vals = malloc(401*sizeof(double)); + warp_vals = malloc(401*sizeof(double)); + idx = 0; + + if ( !ctx->cell ) { + displaywindow_error("No reciprocal unit cell has been found.", ctx->dw); + return 0; + } + + /* Temporarily disable ImageDisplay stuff */ + id = ctx->reproject_id; + ctx->reproject_id = NULL; + + for ( omega_offs=-2.0; omega_offs<=2.0; omega_offs+=0.01 ) { + + double fit, warp; + int i; + Basis cell_copy; + + memcpy(&cell_copy, ctx->cell, sizeof(Basis)); + for ( i=0; iimages->n_images; i++ ) { + ctx->images->images[i].omega += omega_offs; + } + + if ( refine_sequence_random(ctx, &fit, &warp) ) { + printf("RF: Sequencer run failed\n"); + return 0; + } + printf("RF: omega_offs=%f, fit=%f, warp=%f\n", omega_offs, fit, warp); + fit_vals[idx] = fit; + warp_vals[idx++] = warp; + + for ( i=0; iimages->n_images; i++ ) { + ctx->images->images[i].omega -= omega_offs; + } + memcpy(ctx->cell, &cell_copy, sizeof(Basis)); + + } + + ctx->reproject_id = id; + reproject_lattice_changed(ctx); + + window_fit = gtk_window_new(GTK_WINDOW_TOPLEVEL); + gtk_window_set_default_size(GTK_WINDOW(window_fit), 640, 256); + gtk_window_set_title(GTK_WINDOW(window_fit), "Omega-Search Graph: Fit"); + graph_fit = gtk_value_graph_new(); + gtk_value_graph_set_data(GTK_VALUE_GRAPH(graph_fit), fit_vals, idx); + gtk_container_add(GTK_CONTAINER(window_fit), graph_fit); + gtk_widget_show_all(window_fit); + + window_warp = gtk_window_new(GTK_WINDOW_TOPLEVEL); + gtk_window_set_default_size(GTK_WINDOW(window_warp), 640, 256); + gtk_window_set_title(GTK_WINDOW(window_warp), "Omega-Search Graph: Warp"); + graph_warp = gtk_value_graph_new(); + gtk_value_graph_set_data(GTK_VALUE_GRAPH(graph_warp), warp_vals, idx); + gtk_container_add(GTK_CONTAINER(window_warp), graph_warp); + gtk_widget_show_all(window_warp); + + return 0; + +} + static int refine_sequence_sweep(ControlContext *ctx, double *fit, double *warp) { int i; @@ -298,7 +424,7 @@ static int refine_sequence_sweep(ControlContext *ctx, double *fit, double *warp) } -static gint refine_sequence(GtkWidget *step_button, ControlContext *ctx) { +static gint refine_sequence(GtkWidget *widget, ControlContext *ctx) { double omega_offs; GtkWidget *window_fit; @@ -389,6 +515,7 @@ void refine_open(ControlContext *ctx) { GtkWidget *label; GtkWidget *step_button; GtkWidget *sequence_button; + GtkWidget *sequence_random_button; if ( ctx->refine_window ) return; @@ -418,10 +545,15 @@ void refine_open(ControlContext *ctx) { gtk_misc_set_alignment(GTK_MISC(label), 0.0, 0.5); gtk_label_set_markup(GTK_LABEL(label), "Sequencing"); gtk_table_attach_defaults(GTK_TABLE(table), label, 1, 2, 4, 5); - sequence_button = gtk_button_new_with_label("Run Sequencer"); + + sequence_button = gtk_button_new_with_label("Run Sweeping Sequencer"); gtk_table_attach_defaults(GTK_TABLE(table), sequence_button, 1, 2, 5, 6); g_signal_connect(G_OBJECT(sequence_button), "clicked", G_CALLBACK(refine_sequence), ctx); + sequence_random_button = gtk_button_new_with_label("Run Random Sequencer"); + gtk_table_attach_defaults(GTK_TABLE(table), sequence_random_button, 1, 2, 6, 7); + g_signal_connect(G_OBJECT(sequence_random_button), "clicked", G_CALLBACK(refine_random_sequence), ctx); + g_signal_connect(G_OBJECT(ctx->refine_window), "response", G_CALLBACK(refine_response), ctx); gtk_widget_show_all(ctx->refine_window); -- cgit v1.2.3